快速反平方根算法
source link: https://hx-w.github.io/article/39b0/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
什么是反平方根
反平方根即平方根的倒数：
该计算表达式在图形学中的向量正规化中经常用到，对于大规模场景的法线向量计算，如果仅使用
float rsqrt(float number)
{
float y = 1 / sqrt(number);
return y;
}
就显得非常笨拙。
因为在计算机中，一般加法与乘法都是经过硬件优化的，计算速度会很快，求平方根则不同。
为了更快的计算一个浮点数的平方根的倒数，一个更快更奇怪的算法在《雷神之锤3》中被提出（可能更早），该算法也成为了Magic Number的典型案例。
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
//y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
return y;
}
要理解该算法是如何工作的，我们需要掌握一些内容：
- 浮点数在计算机内存中的存储方式
- c语言中的类型转换与重解释
- 牛顿迭代法（求函数的近似零点）
IEEE 754
在IEEE 754标准中，定义了浮点数在计算机中的存储方式：
以32位浮点数为例（在c语言中为float）
其中32位bit被分割成了3部分（科学计数法）：
- Sign 符号位，0为正数，1为负数。
- Exponent 指数位，简写为EE，一共包含8位bit。
- Fraction 或 Mantissa 尾数位，简写为MM，一共包含23bit。
对于求反平方根的算法，符号位可以不进行讨论，因为输入的浮点数只有为正该算法才在实数域中有意义。
对于指数位EE（移码），可以表示-126到127（-127和128被用作特殊值处理），实际表示的2的指数应该是E−127E−127。
对于尾数位MM，表示的是科学计数法中，尾数的小数点后的数字，因为在二进制的科学计数法表示里，尾数的小数点前的数字必为1。
那么，根据IEEE 754的定义，一个浮点数（32位）F就可以表示为：
在本算法中，符号位可以认为恒0，即可化简为：
同时，如果我们将该32位看作是int或long（在c语言中二者可以看做同一类型），那其对应的整形数LL为：
对于同样的32位bit，我们可以将其解释为浮点数FF，也可以解释为整数LL。那么FF和LL之间有没有数量上的关系？
LL与FF的关系
对于FF的表达式，我们可以进行这样的运算：
观察函数f(x)=xf(x)=x与g(x)=log2(x+1)g(x)=log2(x+1)
当0≤x≤10≤x≤1时，f(x)≈g(x)f(x)≈g(x)即可表示为f(x)=g(x)+μf(x)=g(x)+μ，其中μμ为误差。
故log2(1+2−23M)≈2−23M+μlog2(1+2−23M)≈2−23M+μ，即：
又因为L=223E+ML=223E+M，得
至此，我们得到了对于同一32位二进制串的浮点数解释FF与整数解释LL之间的近似数量关系。
c语言中的类型转换与重解释
Type punning (类型双关): c++中可以使用
reinterpret_cast实现双关
在c语言中，如果我们需要将32位浮点数float转为long（或int），我们可以显示得：
float num_float = 3.3;
long num_long = (long)num_float;
得到的结果num_long为3。
但是如果我们更底层的，对于IEEE 754中提到的，我想把32位浮点数，转换成32bit都相同的整形，该如何做呢？
在程序运行中，变量的数据一般是存储在内存中，在c语言里可以通过&来获取变量的地址。而(long
*)
地址即可使用long的指针指向该地址，再使用*获取该指针的内容，即可实现了在c语言里使用指针重解释内存中的数据类型。
即Q_rsqrt函数中的i = * ( long * )
&y;，y本为float类型，使用该方法可以将其bit不变得解释为long类型的i变量。
求解反平方根的近似
在求解之前，我们先回顾上面的结论：
其中FF为32位bit的浮点数解释，LL为32位bit的整形解释。而在c语言中使用指针可以实现FF与LL的转换，更重要的是，转换是极快的。
至此，我们可以明确快速平方根的算法的目的：
- 求解有效32位浮点数的平方根的倒数
- 尽可能的只使用加法、乘法与位运算
- 时间复杂度为常数
- 可以有一些误差（计算机存储浮点数本身就有误差）
值得注意的是，浮点数本身不支持位运算
回到一开始的问题，求解：y=1/x−−√y=1/x
问题可以表示为：
将log2(F)=223L+μ−127log2(F)=223L+μ−127带入上式：
对于上式结果：Ly=32⋅223⋅(127−μ)−12LxLy=32⋅223⋅(127−μ)−12Lx
FxFx表示32位bit的xx的浮点数表示，LxLx表示32位bit的xx的整形表示。在c语言中，二者是可以互相转换的。
前半部分：32⋅223⋅(127−μ)=0x5f3759df32⋅223⋅(127−μ)=0x5f3759df，在μ=0.04505μ=0.04505时成立，更通用的，可以使μ=0.0430μ=0.0430。
后半部分：−12Lx−12Lx在c语言中可以使用位操作实现：-(Lx >>
1)
故使用c语言实现该表达式为：
long Ly = 0x5f3759df - (Lx >> 1);
即原算法实现中的内容。
通过该公式可以求出LyLy，再在c语言中通过Fy = * (float
*) &Ly即可求得FyFy。
在上述推导中，我们用到了log2(x+1)≈x+μlog2(x+1)≈x+μ。
一般情况下μ=0.0430μ=0.0430，该算法实现中取了μ=0.04505μ=0.04505。都是为了使用x+μx+μ逼近logx(x+1)logx(x+1)。
虽然误差总是不可避免，为了减少这种逻辑上带来的误差，我们需要将误差尽可能的减少。
牛顿迭代法是一种用来近似函数零点的方法：
牛顿迭代
对于某复杂函数，我们无法推导出其零点的表达式，或者零点的表达式不利于计算机的计算，我们就可以使用牛顿法来近似得到零点。
牛顿法的核心思想是：切线是曲线的线性逼近
对于某一可导函数，我可以取其定义域内的任意xx：
- 在函数上对应xx的位置求一阶导数，即函数在xx位置的切线。
- 该切线与X轴最多有一个交点，如果有交点的话，我们就以该交点为xx重复求切线，求切线与X轴交点的这个过程。
每次求的切线与X轴的交点都更加逼近函数的真正零点。
根据上图，容易得到：
对于牛顿法求零点，我们需要知道
- 函数的表达式
- 初始迭代值
所以对于本算法，我们需要构造yy的函数，使得其零点为y=1x√y=1x。
经过一些变换，我们就可以写出该函数：
约束y>0y>0，易得：f′(y)=−2⋅y−3f′(y)=−2⋅y−3
使用牛顿法的前提：函数表达式我们有了，初始迭代值我们也在上面求出来了FyFy。
那就可以开始进行牛顿迭代的推导：
我们可以明确上式中：yy是迭代的值，初识迭代值为FyFy，xx为y=1x√y=1x中的xx。
那么就很容易可以转换成c语言中的：
Fy = Fy * (1.5 - 0.5 * x * Fy * Fy);
对应源代码中的：
y = y * ( threehalfs - ( x2 * y * y ) );
在该算法过程中，进行一次迭代即可满足精度的需求。
现在我们就可以把原码进行完整的解释：
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F; // 牛顿迭代中的1.5
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // 将32bit浮点数y重新解释为32bit整形i
i = 0x5f3759df - ( i >> 1 ); // 对i求反平方根的近似
y = * ( float * ) &i; // 把32bit整形i重新解释为32bit浮点数y
y = y * ( threehalfs - ( x2 * y * y ) ); // 第一次牛顿迭代
//y = y * ( threehalfs - ( x2 * y * y ) ); // 第二次迭代，看起来不需要
return y;
}
推广到double类型
以上的算法只针对32位浮点数float，对于64位浮点数double我们也可以应用该思想设计出快速的反平方根算法。
在double类型的64位bit中： 1. 符号位占1bit 2. 指数位EE 占11bit 3. 尾数位MM 占52bit
有了这些数据，我们就可以写出c语言中的函数实现：
double Q_rsqrt(double number)
{
long long i;
double x2, y;
const double threehalfs = 1.5F; // 牛顿迭代中的1.5
x2 = number * 0.5F;
y = number;
i = * ( long long * ) &y; // 将64bit浮点数y重新解释为64bit整形i
i = 0x5fdd3020c49ba400 - ( i >> 1 ); // 对i求反平方根的近似
y = * ( double * ) &i; // 把64bit整形i重新解释为64bit浮点数y
y = y * ( threehalfs - ( x2 * y * y ) ); // 第一次牛顿迭代
//y = y * ( threehalfs - ( x2 * y * y ) ); // 第二次迭代，看起来不需要
return y;
}
其中0x5fdd3020c49ba400根据μμ的实际取值不同而不同。
Recommend
-
96
数据结构与算法—快速排序（Java实现）
-
25
算法原理 下列动图来自 @五分钟学算法 ，演示了快速排序算法的原理和步骤。 步骤: 从数组中选...
-
41
夯实算法-快速排序 2020-03-01 / / 168 点击 / 阅读耗时 7 分钟...
-
40
LeetCode 69. x 的平方根：二分查找法实现自定义的函数：x 的平方根 题目描述 实现 int sqrt(int x) 函数。 计算并返回 x 的平方根，其中 x 是非负整数。 由于返回类型是整数，结果只保留整数的...
-
13
标签：二分 题目：x的平方根 题号：69 题干：实现 int sqrt(int x) 函数。计算并返回 x 的平方根，其中 x 是非负整数。由于返回类型是整数，结果只保留整数的部分，小数部分将被舍去 示...
-
7
求平方根倒数的算法 作者: 张志强 , 发表于 20...
-
3
LeetCode-069-x的平方根发布于 59 分钟前x 的平方根题目描述：实现 int sqrt(int x) 函数。计算并返回 x 的平方根，其中 x 是非负整数...
-
7
迭代法求平方根 作者： physixfan 很久以前看过求平方根的这么一个方法，是用的迭代，按下面这个公式进行迭代即可求出a的平方根： Xn=(Xn-1+a/Xn-1)/2 (n>0) 其中X0可以随便取一个数。当Xn-Xn-1满足精度要...
-
4
#yyds干货盘点# 前端歌谣的刷题之路-第二十八题-平方根 精选 原创 我是歌谣 我有个兄弟 巅峰的时候排名c站总榜19 叫前端小歌谣 曾经...
-
5
旭穹の陋室现代硬件算法[6.5]: 快速平方根倒数发表于2023-08-08|...
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK