Ron的BLOG
0%

解读快速反平方根算法

Posted on

-1. 涉及知识

  1. float单精度的内存表示方式；
  2. 对数和指数；
  3. 牛顿方法/牛顿迭代。

0. 背景

在标准化向量计算的时候，会频繁使用。譬如计算A点到B点的方向，需要计算标准化向量。

1
2
3
4
float len = sqrt(x^2 + y^2 + z^2);
float x_nor = x/len;
float y_nor = y/len;
float z_nor = z/len;

1. 计算下真实的 $1/\sqrt{5}$ = ？答案为：0.44721359

2. 代码计算过程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
// 快速反平方根
float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y = number;
    printf("number为 = %f\n", y);
    i = * ( long * ) &y;         // evil floating point bit hack
    printf("转换为long = %d\n", i);
    i = 0x5f3759df - ( i >> 1 ); // what the fuck?
    printf("long的结果 = %d\n", i);
    y = * ( float * ) &i;
    printf("转换为float = %f\n", y);
    y = y * ( threehalfs - ( x2 * y * y ));  // 1st iteration
    printf("Netown 迭代1次 = %f\n", y);
    y = y * ( threehalfs - ( x2 * y * y ));  // 2nd iteration, can be removed
    printf("Netown 迭代2次 = %f\n", y);
    return y;
}

//number为       = 5.000000
//转换为long     = 1084227584
//long的结果     = 1055349215
//转换为float    = 0.451858
//Netown 迭代1次 = 0.447141
//Netown 迭代2次 = 0.447214

//number = 5.000000, ret = 0.447214

可以发现代码计算的结果 0.447214 与真实结果 0.44721359 非常接近，小数点前五位都是一致的。

3. 步骤说明

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// 快速反平方根
float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

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

    return y;
}

3.1. 第1步：讲number的地址转成long类型地址，并取long类型的值

1
2
3
4
5
6
7
8
i = * ( long * ) &y;
// 让i生成与y一样的内容，但是i是 long类型，y是float
// 不能使用强制类型转换，因为强制类型转换会丢失bit信息
// 此处目的是对float类型的number进行位运算，而float是不能进行位运行的，所以需要转换为long

// &y 获取到y的地址
// ( long * ) &y 然后将float的地址转换为 long的地址
// * ( long * ) &y 然后获取地址上的数值

3.2. 第2步：计算出反平方根的值（long类型下）

此处只有一句代码 i = 0x5f3759df - ( i >> 1 ) 但是这里非常复杂。这里得到的最终结果是 Ly=(3/2)*2^23*(127-μ) - 1/2*Lx (Lx是原来的i，Ly是计算结果，结果重新赋值给了i）。
这里μ为 0.04505 ，所以最终结果为 (3/2)*2^23*(127-μ) - 1/2*Lx 既 Ly = 0x5f3759df - (Lx >>1)

关于此处的解释，后面举例具体来说明。

F与L之间的转换关系为：log2(F)=(2^-23)*L + μ - 127

3.3. 第3步：重新回为浮点数

3.4. 第4步：使用牛顿迭代法提高精度

牛顿迭代是可以提高精度，前面得到了一个靠近的数值，但是由于μ的存在，是存在误差的，这里可以使用牛顿迭代法求到更具体的数值。

牛顿迭代法的公式是：

这里的f(x)是 $f(x) = \frac{1}{x^2}-a$，其导数函数是$f^{‘}(x)=\frac{-2}{x^{-3}}$，把函数和导函数代入牛顿迭代法公式：

1
2
3
y = y * ( threehalfs - ( x2 * y * y ));    // 1st iteration
// 代入值下
y = y * ( 1.5F- ( 0.5F * number * y * y ));   // 1st iteration

4. 举例来说：

number = 5.0F

float的二进制表示为

1
2
3
4
5
6
7
8
9
单精度的float，有32bit，
第1位位符号位，
第2~9位为指数，表示为E
最后23位为尾数，表示为M

1bit     8bit       23bit
符号位S  指数E       尾数M
0       10000001   01000000000000000000000  // 表示为float的时候，值为5.0F
0       10000001   01000000000000000000000  // 表示为long 的时候，值为1084227584

同样的一个32位的地址，可以解读成单精度的float，也可以解读成long类型的整型，下面是公式（F表示float类型的值，L表示long类型的值）：

根据上面两个公式，我们可以得到 float类型的值F，与 long类型的值L之间的关系如下(μ=0.04505）：

我们再把我们想要的内容答案 X = 1/sqrt(number) 代进去。（如下，Fy是我们想要的值）

两边取以2为底的对数，等式是成立的

获取到了上面的公式之后，可以利用 float类型的F与 long类型L的关系，把F换成L

至此，我们得到了代码中第2步的来由。

然后我们将 long 类型的Ly的结果，转成 float 类型，就能得到我们初步的答案。

验证一下：

现在，F = 0.5, L= 1084227584 的情况下，log2(F)=(2^-23)*L + μ - 127 是否成立

log2(F) = 2.3219280949

(2^-23)L + μ - 127 = (2^-23)L + 0.04505 - 127 = 2.29505

此处答案接近。

代入最终计算公式：

Ly=(3/2)2^23(127-μ) - 1/2*Lx

= (3/2)2^23(127-0.04505) - 1/2*1084227584

= 1055349171

1
2
3
4
1bit     8bit       23bit
符号位S  指数E       尾数M
0       10000001   01000000000000000000000
0       01111101   11001110101100110110011

将Ly转为Fy = 0.45185623

牛顿迭代一次：

1
y = y * ( threehalfs - ( x2 * y * y ))

y1 = y ( 1.5- ( 0.55 y y ))
= 0.45185623 ( 1.5- ( 0.55 0.45185623 0.45185623 ))
= 0.44714105083

一些在线二进制转换的工具：

浮点数表示：
https://tooltt.com/ieee/
https://tooltt.com/floatconverter/
进制转换：
https://www.sojson.com/hexconvert.html
LaTeX公式编辑器：
https://www.latexlive.com/