C语言 在 libm 中的日志操作中理解数字文字



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

从这里下载代码

部分代码如下所示。我想知道0x95f640x6147a0x6b851的含义。

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;

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

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

检查是否 ((hx-0x6147a(|(0x6b851-hx((>0 实际上是检查 hx 是否在范围内0x6147a和0x6b851。字0x3ff6147a较高的浮点数约为 1.38,位0x3ff6b851较高的浮点数约为 1.42,即略小于 sqrt(2( 和略高于 sqrt(2(。尚不确定这些数字是否重要。

好吧,如果没有人愿意做出完整的答案 - 我会做不完整的。

我没有太多时间弄清楚这个确切的值来自哪里,所以我的答案将是通用的。这与你在 http://en.wikipedia.org/wiki/Fast_inverse_square_root 中看到的浮点位魔术完全相同。比如说,hx &= 0x000fffff;只从双精度的高字中提取尾数(只有 20 位高字 - 高位是符号和指数(——这个常量在浮点值的某些部分下执行整数位运算(特别是在尾数上,如我所见(。进行这种计算需要花费相当多的努力,但是在libc这样广泛使用的库中,即使是很小的性能提升也可以被认为是显着的。

为什么这样做是因为整数运算比浮点运算快得多。在当前的 CPU 中,浮点数和整数之间的性能差异可能并不大(特别是如果您考虑到一些 SSE 和其他矢量指令 - 尽管并非每种算法都可以从 SIMD 指令中获得性能提升(,但它要高得多。因此,有人在简化公式和在整数而不是浮点数中进行尽可能模糊的计算方面做得令人印象深刻,- 我假设其他人只是复制了这段代码,因为这个常量似乎存在于我可以访问的每个 libc 中。

我知道这不是你一直在寻找的答案 - 对此感到抱歉。您还可以看看精彩的 http://graphics.stanford.edu/~seander/bithacks.html

最新更新