看一下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;
更新:我熟悉十六进制表示法。我有兴趣了解代码的内部工作与body头中描述的算法/方法的关系。为什么要使用这些特定值,以及它的用途是什么?
好吧,如果没有人愿意做出完整的答案 - 我会做不完整的。
我没有那么多时间来弄清楚这些确切的值来自何处,所以我的答案将是通用的。这与你在http://en.wikipedia.org/wiki/Fast_inverse_square_root中看到的浮点魔法完全相同。比方说,hx &= 0x000fffff;
只从双字的高字中提取尾数(只有20位的高位字 - 高位是符号和指数) - 这个常量在浮点值的某些部分下执行整数位操作(特别是在尾数上,如我所见) )。进行这种计算需要花费很多精力,但在像libc这样广泛使用的库中,即使是小的性能提升也可能被认为是重要的。
为什么这样做是因为整数运算比浮点运算快得多。浮点数和整数之间的性能差异在当前的CPU中并不是那么大(特别是如果你考虑一些SSE和其他向量指令 - 虽然不是每个算法都能从SIMD指令中获得性能提升),但它要高得多。所以有人在简化公式方面做了令人印象深刻的工作,并在整数而不是浮点数中尽可能地进行计算,而且我假设其他人都只是复制了这些代码,因为这个常量似乎存在于我有权访问的每个libc中。
我知道这不是你一直在寻找的答案 - 抱歉。你也可以看看辉煌的http://graphics.stanford.edu/~seander/bithacks.html
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)。尚不确定这些数字是否重要。