c - 在 libm 的日志操作中理解数字文字

Question

查看libm中日志操作的实现，有一些数字文字我理解有问题。

从here下载代码

部分代码如下所示。我想知道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;

更新:我熟悉十六进制符号。我有兴趣了解代码的内部工作与正文 header 中描述的算法/方法的关系。为什么使用这些特定值，使用目的是什么？

Answer 1

sqrt(2)的iee754表示的高32位字是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的 float 约为1.38，高位0x3ff6b851的 float 约为1.42，即略小于sqrt(2)，略大于sqrt(2)。尚不确定这些数字是否重要。