«Магический» метод вычисления обратного квадратного корня, по-видимому, восходящий к игре Quake, описан во многих источниках. В Википедии есть хорошая статья об этом: https://en.wikipedia.org/wiki/Fast_inverse_square_root
В частности, я обнаружил, что следующее очень хорошее описание и анализ алгоритма: https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf
Я пытаюсь воспроизвести некоторые из этих результатов в этой статье, но у меня проблемы с точностью. Алгоритм, закодированный на C, имеет следующий вид:
#include <math.h>
#include <stdio.h>
float Q_rsqrt(float number) {
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(long *) &y;
i = 0x5f3759df - (i >> 1);
y = *(float *) &i;
y = y * (threehalfs - (x2 * y * y));
// y = y * (threehalfs - (x2 * y * y));
return y;
}
В документе говорится, что относительная ошибка не превышает 0.0017522874 для всех положительных нормальных чисел с плавающей запятой. (См. Код в Приложении 2 и обсуждение в Разделе 1.4.)
Однако, когда я «подключаю» число 1.4569335e-2F, получаемая мной ошибка превышает этот прогнозируемый допуск:
int main ()
{
float f = 1.4569335e-2F;
double tolerance = 0.0017522874;
double actual = 1.0 / sqrt(f);
float magic = Q_rsqrt(f);
double err = fabs (sqrt(f) * (double) magic - 1);
printf("Input : %a\n", f);
printf("Actual : %a\n", actual);
printf("Magic : %a\n", magic);
printf("Err : %a\n", err);
printf("Tolerance: %a\n", tolerance);
printf("Passes : %d\n", err <= tolerance);
return 0;
}
Результат:
Input : 0x1.dd687p-7
Actual : 0x1.091cc953ea828p+3
Magic : 0x1.08a5dcp+3
Err : 0x1.cb5b716b7b6p-10
Tolerance: 0x1.cb5a044e0581p-10
Passes : 0
Итак, этот конкретный ввод, кажется, нарушает утверждение, сделанное в этой статье.
Мне интересно, проблема ли это в самой статье, или я допустил ошибку в коде. Буду признателен за любые отзывы!
