转: 关于黑魔法常数

转自： http://www.guokr.com/post/91389/

看了之前求0x5f3759df的数学原理之后花了两天时间把@Sheldon 给的那篇论文看完了，感谢@Uroboros 的讲解，作为一个没什么计算机背景的人，能看完（好吧其实没全看）这篇文章我还是非常有成就感的，所以写下这个帖子，基本上全是翻译了，了解计算机的人可以抱着挑错的心态看。

基础知识1：i>>1
操作i>>1表示将二进制数i向右移动一位，也就是将最后一位删掉并在最前一位添加0
注意到我们将一个n位的十进制数M删掉最后一位之后就变成了n-1位，可以看做将这个十进制数除以10之后向下取整floor(M/10)
同样的，讲一个二进制数删掉最后一位之后相当于将这个数除以2并向下取整floor(M/2)

这样看上去i>>1就像是floor(i/2)，因为函数f(x)=1/sqrt(x)的一阶导数是-1/2*(x)^-3/2，正好有个-1/2在前面，不禁让人感觉 0x5f3759df - ( i >> 1 )是函数1/sqrt(x)在某一个点的一阶泰勒展开。在Fast Inverse Square Root 里面有这样一段：

Eberly’s explanation was that this produced linear approximation, but is incorrect; we’ll see the guess is piecewise linear, and the function being approximated is not the same in all cases.

“Eberly的解释是说这相当于做了线性近似，但是这个解释是不对的。我们会看到这个估计值是分段线性的，并且这个近似函数在各种情况下也并不相同。”

为什么这么说呢？这里需要用到基础知识2：浮点数存储方式
各种类型浮点数的存储方式可以通过查看IEEE745（完全不知道是什么东西）了解
这里用到的是32位单精度浮点数，并且总是正数，表示方式为：
0|E|M
其中0代表正数
E是指数，E=0相当于2^-127
M表示一个绝对值小于1的数，但是注意到这里省略了一位。也就是说，当转化位十进制的时候应该表示为（1+M）
那么换算之后的数就应该是：（1+M）2^（E-127）
这样我们发现其实i>>1并不完全是floor(i/2)而是将一个数变成(floor(M/2)+1)*2^floor((E-127)/2)
而且根据E的奇偶性尾数可能还需要加上1/2

不妨令e=E-127，注意到1/sqrt(x)是让原数的指数变为-e/2，这么说来卡马克可能仅仅是希望产生一个e/2而用上了位移，接下来就是要找到一个相减之后产生-e/2并且让尾数尽量和(1+M)^-1/2相近

由于这个数必然为正数假设这个数值为：
0|R1|R2
接下来便是要分情况讨论：
假设E为偶数，这时候指数的右移并不会影响到尾数的数值：
这时候e是奇数，令e=2d+1
那么相减之后指数部分变为：

注意这里的相减其实是将浮点数转化为整型（也就是正则化）之后再相减，而不是普通的浮点数加减。
由于初始值一定要是正数，所以我们需要上式一定为正，因为0=<E<=256，所以R1>=128
因为我们讨论的E为偶数，也就是末位数一定是0，所以不用考虑他右移后对M的影响，所以两数相减之后的结果是：

注意这里用M/2而不是floor(M/2)因为这一点点的误差相较于其他误差来说太小了
当然，还存在一种情况，那就是R2<M/2，其实二进制的加减和十进制差不多，如果尾数小了，那么就要像更高位数“借位”，也就是说这种情况下相减之后的结果是：

我们定义：

那么我们可以将这两种情况合并为一个函数：

这个函数就是我们对函数1/sqrt(x)的近似了，那么我们的目标就是让这个函数的相对误差|(y-y0)/y|尽量小：

这样我们得到一个误差方程：

之后：

注意这里的1/8其实是凑出来的，具体凑法可以先假设一个小于一的正数a，由于0<R2<1，0<M<1，我们可以通过展开得到R1关于a的一个区间。让a尽量小，使得这个区间范围小于一。根据R1一定是整数的特性，我们可以确定使得误差最小的R1。这里得出R1=190，带入十六进制里面并右移（注意开头有一个表示符号的0）就是0x5f，正好是黑魔法常数的头几位。