左耳听风00013-1

你好，我是陈皓，网名左耳朵耗子。

下列代码是在《雷神之锤III竞技场》源代码中的一个函数（已经剥离了C语言预处理器的指令）。其实，最早在2002年（或2003年）时，这段平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了，并且也在程序员圈子里引起了热烈的讨论。

我先把这段代码贴出来，具体如下：

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y; // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );  // what the fuck? 
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );  // 1st iteration 
    // 2nd iteration, this can be removed
    // y  = y * ( threehalfs - ( x2 * y * y ) ); 

    return y;
}

这段代码初读起来，我是完全不知所云，尤其是那个魔数0x5f3759df，根本不知道它是什么意思，所以，注释里也是 What the fuck。今天这节课，我主要就是想带你来了解一下这个函数中的代码究竟是怎样出来的。

其实，这个函数的作用是求平方根倒数，即 $x^{-1/2}x−1/2，也就是下面这个算式：$

$\frac{1}{\sqrt{x}}x1$

当然，它算的是近似值。只不过这个近似值的精度很高，而且计算成本比传统的浮点数运算平方根的算法低太多。在以前那个计算资源还不充分的年代，在一些3D游戏场景的计算机图形学中，要求取照明和投影的光照与反射效果，就经常需要计算平方根倒数，而且是大量的计算——对一个曲面上很多的点做平方根倒数的计算。也就是需要用到下面的这个算式，其中的x,y,z是3D坐标上的一个点的三个坐标值。

$\frac{1}{\sqrt{x^{2}+y^{2}+z^{2}}}x2+y2+z21$

基本上来说，在一个3D游戏中，我们每秒钟都需要做上百万次平方根倒数运算。而在计算硬件还不成熟的时代，这些计算都需要软件来完成，计算速度非常慢。

我们要知道，在上世纪90年代，多数浮点数操作的速度更是远远滞后于整数操作。所以，这段代码所带来的作用是非常大的。

计算机的浮点数表示

为了讲清楚这段代码，我们需要先了解一下计算机的浮点数表示法。在C语言中，计算机的浮点数表示用的是IEEE 754 标准，这个标准的表现形式其实就是把一个32bits分成三段。

第一段占1bit，表示符号位。代称为S（sign）。
第二段占8bits，表示指数。代称为E（Exponent）。
第三段占23bits，表示尾数。代称为M（Mantissa）。

如下图所示：

然后呢，一个小数的计算方式是下面这个算式：

$(-1)^{S}\ast(1+\frac{M}{2^{23}})\ast 2^{(E-127)}(−1)S∗(1+223M)∗2(E−127)$

但是，这个算式基本上来说，完全就是让人一头雾水，摸不着门路。对于浮点数的解释基本上就是下面这张漫画里表现的样子。

下面，让我来试着解释一下浮点数的那三段表示什么意思。

第一段符号位。对于这一段，我相信应该没有人不能理解。
第二段指数位。什么叫指数？也就是说，对于任何数x，其都可以找到一个 $nn，使得2^{n}2n<=x<=2^{n+1}2n+1。比如：对于3来说，因为 2 < 3 < 4，所以 n=1。而浮点数的这个指数为了要表示0.00x的小数，所以需要有负数，这8个bits本来可以表示0-255。为了表示负的，取值要放在 [-127,128] 这个区间中。这就是为什么我们在上面的公式中看到的 2^{(E-127)}2(E−127)这一项了。也就是说，n = E-127n=E−127，如果n=1n=1，那么EE就是128了。$
第三段尾数位。也就是小数位，但是这里叫偏移量可能好一些。这里的取值是在[ 0 - $2^{23}223]中。你可以认为，我们把一条线分成2^{23}223个线段，也就是8388608个线段。也就是说，把2^{n}2n到2^{n+1}2n+1分成了8388608个线段。而存储的M值，就是从2^n2n到 x 要经过多少个段。这要计算一下，2^{n}2n到x的长度占2^{n}2n到2^{n+1}2n+1长度的比例是多少。$

我估计你对第三段还是有点不懂，那么我们来举一个例子。比如说，对3.14这个小数。

是正数。所以，S = 0。
$2^121 < 3.14 <2^222。所以，n=1， n+127 = 128。所以，E=128。$
(3.14 - 2) / (4 - 2) = 0.57，而 $0.57*2^{23} = 4781506.560.57∗223=4781506.56，四舍五入，得到M = 4781507。因为有四舍五入，所以，产生了浮点数据的精度问题。$

把S、E、M转成二进制，得到 3.14的二进制表示。

我们再用IEEE 754的那个算式来算一下：

${(-1)}^0*({1+\frac{4781507}{2^{23}}})*2^{(128-127)}(−1)0∗(1+2234781507)∗2(128−127)$

你看，浮点数的精度问题出现了。

我们再来看一个示例，小数 0.015。

是正数。所以，S = 0。
$2^{-7}< 0.015 < 2^{-6}2−7<0.015<2−6 。所以，n=-7， n+127 = 120。所以，E=120。$
$2^{-7}) / (2^{-6} - 2^{-7})(0.015−2−7)/(2−6−2−7) = 0.0071875/0.0078125=0.920.0071875/0.0078125=0.92。而0.92 * 2^{23} = 7717519.360.92∗223=7717519.36，四舍五入，得到 M = 7717519。$

于是，我们得到0.015的二进制编码：

其中：

120 的二进制是01111000
7717519的二进制是11101011100001010001111

返回过来算一下：

$(-1)^{0}\ast (1+\frac{7717519}{2^{23}})\ast 2^{(120-127)}(−1)0∗(1+2237717519)∗2(120−127)$

你看，浮点数的精度问题又出现了。

我们来用C语言验证一下：

int main() {
    float x = 3.14;
    float y = 0.015;
    return 0;
}

在我的Mac上用lldb 工具 Debug 一下。

(lldb) frame variable
(float) x = 3.1400001
(float) y = 0.0149999997

(lldb) frame variable -f b
(float) x = 0b01000000010010001111010111000011
(float) y = 0b00111100011101011100001010001111

从结果上，完全验证了我们的方法。

好了，不知道你看懂了没有？我相信你应该看懂了。

posted @ 2023-01-12 11:31 易先讯阅读(12) 评论(0) 收藏举报

刷新页面返回顶部

易先讯

关注程序员个人发展

左耳听风00013-1

计算机的浮点数表示