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

Читая Уловки гуру программирования 3D-игр, я наткнулся на эту функцию сортировки, написанную на встроенном ассемблере:

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);
}

x86 optimization assembly square-root sar
person vexe    schedule 21.01.2017    source источник
comment
0x3F800000 — 32-битное представление с плавающей запятой для версии 1.0.   -  person mbil    schedule 22.01.2017
comment
Интересно, вот почему я предполагаю, что получаю неверные результаты, когда меняю значение параметра на int? Похоже, функция работает только для поплавков?   -  person vexe    schedule 22.01.2017
comment
Что еще более важно, это смещение экспоненты. Таким образом, он отменяет смещение, уменьшает показатель степени вдвое, а затем снова добавляет смещение. Это также немного портит мантиссу.   -  person Jester    schedule 22.01.2017
comment
@vexe Что касается общего принципа, вы можете посмотреть этот вопрос   -  person njuffa    schedule 27.01.2017

Ответы (4)


arrow_upward
10
arrow_downward

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

Примените базовую математику,

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%, если входное значение близко к нечетной степени двойки. Для четных степеней двойки результат точен.

person Jester    schedule 22.01.2017

arrow_upward
0
arrow_downward

0x3F800000 в float равно 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).

person D Summy    schedule 21.01.2017
comment
Если вы нарисуете графики для y=(x+1)/2 и y=sqrt(x), вы увидите, что они близки, когда x находится в [1,2]. Так что я предполагаю, что это приближение для значений, которые лежат в этом интервале. - person Roadowl; 22.01.2017
comment
@Roadowl Он не вычисляет (x+1)/2 - person Jester; 22.01.2017
comment
В статье Википедии, на которую вы ссылаетесь, обсуждается совершенно другой метод приближения квадратного корня, чем код, представленный в вопросе. - person Cody Gray; 22.01.2017

arrow_upward
0
arrow_downward

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

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
person mbil    schedule 21.01.2017
comment
Было бы лучше, если бы вы также показывали поля экспоненты/мантиссы для каждого шага, а не только весь двоичный 32-битный шаблон. - person Peter Cordes; 21.04.2018

arrow_upward
0
arrow_downward

Плавающее число f = (1 + m)*[2^(e+127)], где m — часть мантиссы, e — экспоненциальная часть.

таким образом: sqrt (f) = (f) ^ (1/2) = ((1 + m) * [2 ^ (e + 127)] ) ^ (1/2)

-> ((1 + m)* [2^(e+127)] )^(1/2) = (1 + m)^(1/2) * 2^((e + 127)/2)

В части показателя степени 2^((e + 127)/2):

2^((e + 127)/2) = 2^( (e-127/2) + 127)

таким образом, в плавающем представлении это (e - 0x3F800000)/2 + 0x3F800000

В части мантиссы (1 + m)^(1/2):

из формулы биномиального ряда (1 + x)^r = 1 + rx + (r(r - 1)/2)*(x^2) +....

Таким образом, (1 + m) ^ (1/2) равно (1 + m/2 - (m ^ 2)/8 + ...) оно ПРИБЛИЗИТЕЛЬНО равно 1 + m/2 (типичное приближение к первому порядку). Поэтому , часть мантиссы должна быть разделена на 2.

Однако мантисса и показатель степени объединяются как число, сдвиг вправо делит показатель степени и мантисса ОБА.

Чтобы оценить ошибку, вы можете рассмотреть вторые члены биномиального ряда - (m ^ 2)/8.

Поскольку m всегда меньше 1, я подставляю m как 0,9999 (0,5 + 0,25 + 0,125 + ...)

(m^2)/8 = 0,12497500125, это наихудший случай.

person Gaiger Chen    schedule 20.04.2018