Осмысление числовых литералов в работе журнала в libm

Question

Глядя на реализацию операции журналирования в libm, есть некоторые числовые литералы, которые у меня проблемы с пониманием.

Загрузите код отсюда

Часть кода показана ниже. Я хотел бы знать значение 0x95f64, 0x6147a а также 0x6b851,

if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */
k += (i>>20);
f = x-1.0;
if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
if(f==zero) { if(k==0) return zero;  else {dk=(double)k;
                           return dk*ln2_hi+dk*ln2_lo;}}
    R = f*f*(0.5-0.33333333333333333*f);
    if(k==0) return f-R; else {dk=(double)k;
             return dk*ln2_hi-((R-dk*ln2_lo)-f);}
}
s = f/(2.0+f); 
dk = (double)k;
z = s*s;
i = hx-0x6147a;
w = z*z;
j = 0x6b851-hx;
t1= w*(Lg2+w*(Lg4+w*Lg6)); 
t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 
i |= j;
R = t2+t1;

ОБНОВЛЕНИЕ: я знаком с шестнадцатеричной нотацией. Я заинтересован в понимании внутренней работы кода в связи с описанным алгоритмом / методом в заголовке тела. Почему использование этих конкретных значений, и какова цель его использования?

Answer 1

Старшее 32-битное слово представления iee754 в sqrt(2) равно 0x3ff6a09e, где старшие 12 бит (то есть 0x3ff) обозначают показатель степени, а младшие 20 бит 0x6a09e обозначают первую часть мантиссы. (1<<20)-0x6a09e равно 0x95f62. В части алгоритма используется число 0x95f64, мы проверяем, если после удаления всех степеней 2 (что делает x в диапазоне 1..2), мы все еще имеем x>sqrt(2), и в этом случае мы делим x на 2. Однако мне не ясно, почему используется 0x95f64, а не 0x95f62.

Эта часть

i = hx-0x6147a;
w = z*z;
j = 0x6b851-hx;
t1= w*(Lg2+w*(Lg4+w*Lg6));
t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
i |= j;
R = t2+t1;
if(i>0) {

Имеет следующий комментарий в источниках

 /*  In order to guarantee error in log below 1ulp, we compute log
  *  by
  *      log(1+f) = f - s*(f - R)    (if f is not too large)
  *      log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)/

Проверка if ((hx-0x6147a)|(0x6b851-hx))>0 на самом деле является проверкой, находится ли hx в диапазоне 0x6147a и 0x6b851. Число с плавающей запятой с более высоким словом 0x3ff6147a составляет около 1,38, число с плавающей запятой с старшими битами 0x3ff6b851 составляет около 1,42, то есть немного меньше, чем sqrt(2) и немного больше, чем sqrt(2). Еще не уверен, значат ли эти цифры.

Answer 2

Хорошо, если никто не хочет дать полный ответ - я сделаю неполный.

У меня не так много времени, чтобы выяснить, откуда пришли эти точные значения, поэтому мой ответ будет общим. Это та же магия с плавающей точкой, которую вы можете увидеть в http://en.wikipedia.org/wiki/Fast_inverse_square_root. Как, скажем, hx &= 0x000fffff; извлекает только мантиссу из старшего слова double (только 20 битов старшего слова - старшие биты являются знаком и показателем степени) - эта константа выполняет операции с целочисленными битами в некоторых частях значения с плавающей запятой (в частности, для мантиссы, как я вижу). Для такого рода вычислений требуется немало усилий, но в столь широко используемой библиотеке, как libc, даже небольшой прирост производительности можно считать значительным.

Это делается потому, что целочисленные операции выполняются намного быстрее, чем операции с плавающей запятой. Возможно, разница в производительности между числами с плавающей запятой и целочисленными значениями не так велика в современных процессорах (особенно если принять во внимание некоторые SSE и другие векторные инструкции - хотя не каждый алгоритм может получить повышение производительности от инструкций SIMD), но она была намного выше. Таким образом, кто-то проделал впечатляющую работу по упрощению формул и выполнению как можно большего количества вычислений в целых числах, а не в числах с плавающей точкой, - и я предполагаю, что все остальные просто скопировали этот код, поскольку эти константы, по-видимому, присутствуют в каждой библиотеке libc, к которой у меня есть доступ.

Я знаю, что это не тот ответ, который вы искали - извините за это. Вы также можете взглянуть на блестящий http://graphics.stanford.edu/~seander/bithacks.html