Как работает эта встроенная функция сборки в sqrt?

Question

Чтение через уловок 3D Game Programming Gurus, я наткнулся на эту функцию сортировки, написанной в инлайн сборки:

inline float FastSqrt(float Value) 
{ 
    float Result; 

    _asm 
    { 
     mov eax, Value 
     sub eax, 0x3F800000 
     sar eax, 1 
     add eax, 0x3F800000 
     mov Result, eax 
    } 

    return(Result); 
} 

Это приближение фактического квадратного корня, но точность достаточно для мои потребности.

Как это работает? Что это за волшебство 0x3F800000 значение? Как мы достигаем квадратного корня путем вычитания, вращения и добавления?

Вот как это выглядит в коде C/C++:

inline float FastSqrt_C(float Value) 
{ 
    float Result; 

    long Magic = *((long *)&Value); 
    Magic -= 0x3F800000; 
    Magic >>= 1; 
    Magic += 0x3F800000; 
    Result = *((float *)&Magic); 

    return(Result); 
} 

Answer 1

Многие люди отмечают, что 0x3f800000 является представление 1.0. Хотя это верно, это не имеет никакого отношения к тому, как работает расчет. Чтобы понять это, вам нужно знать, как хранятся неотрицательные поплавки. f = (1+m)*2^x, с 0 <= m < 1 и m являясь мантиссой, x экспонентом. Также обратите внимание, что x хранится с предубеждением, поэтому то, что на самом деле находится в двоичном формате, равно x+127. 32-битное значение составлено из знакового бита (который равен нулю в нашем случае), за которым следуют 8 бит хранителя экспоненты x+127 и, наконец, 23 бита мантиссы, m. (См. wikipedia article).

Применить основы математики,

sqrt(f) = sqrt((1+m)*2^x) 
     = sqrt(1+m)*sqrt(2^x) 
     = sqrt(1+m)*2^(x/2) 

Итак, как грубое приближение, мы должны вдвое сократить показатель, но из-за смещения мы не можем просто сделать x/2 нам нужно (x-127)/2 + 127. Этот 127 переместился в соответствующее положение бит - это волшебство 0x3f800000.

Разделение на 2 достигается с помощью сдвига вправо на один бит. Так как это действует на весь плавающий, это также оказывает побочное влияние на мантису.

Во-первых, предположим, что оригинальный показатель был ровным. Тогда младший значащий бит, который сдвигается, равен нулю. Таким образом, мантисса также сокращается вдвое, поэтому в итоге мы получаем: sqrt(f) = (1+m/2)*2^(x/2). Мы получили показатель правильно, но мантисса (1+m/2) вместо sqrt(1+m). Максимальная относительная погрешность для этого равна (1.5 - sqrt(2))/sqrt(2) ~ 6%, которая возникает, если m составляет почти 1, что означает, что значение f близко, но меньше нечетной степени 2. Возьмем, например, f=7.99. Формула дает нам около 2.998 вместо 2.827, который действительно имеет ошибку 6%.

Теперь, если показатель экспоненциальности был нечетным, то младшим значащим битом будет 1, и это при смещении в мантиссу приведет к увеличению вдвое. Таким образом, мы получаем sqrt(f) = (1.5+m/2)*2^((x-1)/2). Максимальная погрешность для этого на самом деле равна m=0, и это будет (1.5/sqrt(2)-sqrt(1))/sqrt(1), что снова составляет 6%. Это происходит для чисел, близких к нечетной степени из двух сверху.

Эти два случая означают, что худшая погрешность составляет около 6%, если входное значение оказывается вблизи нечетной мощности двух. Для четных степеней двух результат является точным.

Answer 2

0x3F800000 в поплавке 1. Это происходит потому, что пути поплавков хранятся. Вы можете увидеть визуальное представление на https://gregstoll.dyndns.org/~gregstoll/floattohex/.

Это хорошее приближение, я верю в sqrt. Происхождение этого происходит из игры Quake для обратного sqrt (https://en.wikipedia.org/wiki/Fast_inverse_square_root#Aliasing_from_floating_point_to_integer_and_back).

Answer 3

Вот пример механики это в действии:

FastSqrt (4,0) == 2,0

4.0 to hex -> 0x40800000 
0x40800000 - 0x3f800000 = 0x1000000 
0x1000000 to binary -> 00000001 00000000 00000000 00000000 
shift toward the lsb (sar) -> 00000000 10000000 00000000 00000000 
00000000 10000000 00000000 00000000 back to hex -> 0x00800000 
0x00800000 + 0x3f800000 = 0x40000000 
0x40000000 to dec -> 2.0