Making sense of numeric literals in log operation in libm

Question

Looking at the implementation of log operation in libm, there are some numeric literals that I have problem understanding.

Download the code from here

Part of the code is shown below. I would like to know the meaning of 0x95f64, 0x6147a and 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;

UPDATE: I am familiar with the hex notation. I am interested in understanding the inner working of the code in relationship with the described algorithm/method in the body header. Why the usage of these specific values, and what is the purpose of its usage?

possible duplicate of [What do numbers using 0x notation mean?](http://stackoverflow.com/questions/8186965/what-do-num-bers-using-0x-notation-mean) — Bernhard Barker, Sep 26 '13 at 10:21
'Magic' with floating-point bit representation (http://en.wikipedia.org/wiki/Floating_point). It is tricky (while quite well-commented) and serves only optimisation purpose. — keltar, Sep 26 '13 at 10:51
I am familiar with the hex notation. I am interested in understanding the inner working of the code in relationship with the described algorithm/method in the body header. Why the usage of these specific values, and what is the purpose of its usage? — user2818825, Sep 27 '13 at 06:12
@keltar please extend this to an answer because I believe it does answer the question — Ivaylo Strandjev, Sep 27 '13 at 08:32

score 3 · Answer 1 · answered Feb 20 '19 at 17:04

The higher 32-bit word of iee754 representation of sqrt(2) is 0x3ff6a09e, where the highest 12 bit (i.e. 0x3ff) stand for exponent, and the lower 20 bit 0x6a09e stand for the first part of mantissa. (1<<20)-0x6a09e is 0x95f62. In the part of algorithm the number 0x95f64 is used, we check if after removing all powers of 2 (which makes x in range 1..2) we still have x>sqrt(2), in which case we divide x by 2. It is not, however, clear to me why 0x95f64 and not 0x95f62 is used.

This part

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) {

Has the following comment in the sources

 /*  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)/

The check if ((hx-0x6147a)|(0x6b851-hx))>0 is actually a check if hx is in range 0x6147a and 0x6b851. The floating point number with higher word 0x3ff6147a is about 1.38, the floating point number with higher bits 0x3ff6b851 is about 1.42, i.e. slightly less than sqrt(2) and slightly more than sqrt(2). Not yet sure if these numbers are significant.

score 0 · Answer 2 · answered Sep 27 '13 at 09:51

Ok, if noone willing to make complete answer - i'll make incomplete.

I don't have so much time to figure out where this exact values came from, so my answer will be generic. It's quite the same float-bit magic as you can see in http://en.wikipedia.org/wiki/Fast_inverse_square_root. Like, say, hx &= 0x000fffff; extracts only mantissa from high-word of double (only 20 bits of highword - higher bits are sign and exponent) - this constant performs integer bit operations under some parts of floating point value (specifically on mantissa, as i see). It takes quite a lot of effort to make this kind of calculations, but in so widely used library as libc even small performance gains could be considered significant.

Why it's done is because integer operations are much faster then floating points one. Maybe performance difference between floats and ints is not this big in current CPUs (especially if you take in account some SSE and other vector instructions - although not every algorithm could get performance boost from SIMD instructions), but it was much higher. So someone have made impressive job in simplifying formula and making as mush as possible calculations in integers instead of floats, - and i assume everyone else just copied this code, as this constants seems to be present in every libc i have access to.

I know it's not an answer you've been looking for - sorry about that. You could also take a look at brilliant http://graphics.stanford.edu/~seander/bithacks.html

Making sense of numeric literals in log operation in libm

2 Answers2