快速平方根倒数算法深度理解

快速平方根倒数算法是什么？
简单来说这个算法避开了开方和除法运算快速实现了
y = 1 x y= \frac{1}{\sqrt x} y=x 1

快速平方根倒数算法首次出现是在著名的雷神之锤3的图形计算优化代码中
话不多说先上c代码
一.代码

小白的我以为这么写的：
float y = 1/sqrt(x);开发游戏的大佬是这么写的：
float Q_rsqrt(float number)
{long i; float x2,y;const float threehalfs = 1.5F;x2= number * 0.5F;y= number;i= * (long * ) &y;i= 0x5f3759df-(i>>1);y= * (float * ) &i;//上述代码给牛顿迭代法提供了一个较好的初始值y=y * (threehalfs - (x2*y*y)); //牛顿迭代一次   //y=y * (threehalfs - (x2*y*y))    可持续迭代
}

经大佬测试代码比(float)(1/sqrt(x))快4倍。

初步观看代码：
①定义了32位整型变量i
②定义了两个32位浮点数
③将1.5存入一个常量
④将形参number除2存入x2
⑤将形参number存入y

之后的四行代码乍一看？？？莫慌下面开始正文

二.理解算法前的知识铺垫
1.IEEE 754 （二进制浮点数算术标准）

32位分为三个部分：符号位指数部分

①符号位：0为正数 1为负数

②8位指数部分：8位可以表示0到255所有的数字，但我们还需要得到负指数。
因此我们对0-255减去127就可以得到-127到128之间所有的整数。
举个例子如果指数部分我想要4
那这八位必须表示4+127=131 即10000011
这样通过移码我们可以将正负指数都变成0-255之间的整数
③23位尾数：利用23位我们可以表示0到2^23-1之间的整数
对于十进制科学计数法，尾数是1到10
对于二进制科学计数法，尾数是1-2
11000 = 1.1 ∗ 2 4 11000=1.1*{2}^{4} 11000=1.1∗24
0.0101 = 1.01 ∗ 2 − 3 0.0101=1.01*{2}^{-3} 0.0101=1.01∗2−3
在科学计数法中，第一位按照定义不可以是0，对于二进制这个数只能是1。因此我们不存储这一位。

说了一堆理论的我们举一个例子看一下浮点数是怎么存储在计算机中的：
我们拿 8.25 举例子

Ⅰ.写成二进制：1000.01
Ⅱ.变成科学计数法： 1.00001 ∗ 2 3 = ( 1 + 0.00001 ) ∗ 2 3 1.00001 *{2^3}=(1 + 0.00001) * {2^3} 1.00001∗23=(1+0.00001)∗23
Ⅲ.解释
指数部分我们用 8bit 表示 3+127
3对应的就是2^3中的3
(1 + 0.00001)：加号前面的1我们丢掉不存
对于0.00001中小数点后的数我们用23位尾数去存

罗嗦了这么久现在我们可以轻松理解以下内容：
对于一个浮点数x(十进制或者二进制)有：
x = ( − 1 ) s ∗ ( 1 + m ) ∗ 2 e x=(-1)^s*(1+m)*{2}^{e} x=(−1)s∗(1+m)∗2e

对于下式 ( 1 + 0.00001 ) ∗ 2 3 (1 + 0.00001) * {2^3} (1+0.00001)∗23
m : 00001 m:00001 m:00001
e : 3 e:3 e:3
我们将存取之后的8位指数部分用E表示
将存取之后的23位尾数部分用M表示

E : 3 + 127 E:3+127 E:3+127
M : 00001 ⋯ （补 18 个 0 填满 23 位） M:00001\cdots（补18个0填满23位） M:00001⋯（补18个0填满23位）
我们可以得出E e和M m的关系式
E = e + 127 E=e+127 E=e+127
M = 2 23 ∗ m M={2}^{23}*m M=223∗m
三.数学推导
y = 1 x y= \frac{1}{\sqrt x} \quad y=x 1
取对数：
① log ⁡ 2 y = − 0.5 ∗ log ⁡ 2 ①\log_2 y=-0.5* \log_2 ①log2y=−0.5∗log2
将y和x分别用二进制科学计数法表示：
② y = ( 1 + m y ) ∗ 2 e y ②y=(1+m_y)*{2}^{e_y} ②y=(1+my)∗2ey
③ x = ( 1 + m x ) ∗ 2 e x ③x=(1+m_x)*{2}^{e_x} ③x=(1+mx)∗2ex

联合上述三个式子化简：
l o g 2 ( 1 + m y ) + e y = − 0.5 ∗ ( l o g 2 ( 1 + m x ) + e x ) log_2 (1+m_y)+e_y=-0.5*(log_2(1+m_x)+e_x) log2(1+my)+ey=−0.5∗(log2(1+mx)+ex)
通过一次线性近似：
m y + b + e y = − 0.5 ∗ ( m x + b + e x ) m_y+b+e_y=-0.5*(m_x+b+e_x) my+b+ey=−0.5∗(mx+b+ex)
m y = M y 2 23 m_y= \frac{M_y}{{2}^{23}} \quad my=223My
m x = M x 2 23 m_x= \frac{M_x}{{2}^{23}} \quad mx=223Mx
E y = e y + 127 E_y=e_y+127 Ey=ey+127
E x = e x + 127 E_x=e_x+127 Ex=ex+127
上述五个式化简：
M y + 2 23 ∗ E y = 1.5 ∗ 2 23 ∗ ( 127 − b ) − 0.5 ∗ ( M x + 2 23 ∗ E x ) M_y+{2}^{23}*E_y=1.5*{2}^{23}*(127-b)-0.5*(M_x+{2}^{23}*E_x) My+223∗Ey=1.5∗223∗(127−b)−0.5∗(Mx+223∗Ex)
令
I y = M y + 2 23 ∗ E y I_y=M_y+{2}^{23}*E_y Iy=My+223∗Ey
I x = M x + 2 23 ∗ E x I_x=M_x+{2}^{23}*E_x Ix=Mx+223∗Ex
对于 1.5 ∗ 2 23 ∗ ( 127 − b ) 1.5*{2}^{23}*(127-b) 1.5∗223∗(127−b)我们令其为K，可以知道通过调节b的大小就可以改变K的值

前辈们给出了当
b = 0.0450465 b=0.0450465 b=0.0450465时有K=
K = 0 x 5 f 3759 d f K=0x5f3759df K=0x5f3759df
这刚好对应代码中：

i= 0x5f3759df-(i>>1);

通过一系列不怎么复杂的数学推导我们究竟得到了什么?
我们知道了一个输出结果的整形和输入整形之间的关系，即：
I y = K − 0.5 ∗ I x I_y=K-0.5*I_x Iy=K−0.5∗Ix
四.代码解释
①

long i;
float x2,y;
const float threehalfs = 1.5F;
x2= number * 0.5F;
y= number;
i= * (long * ) &y;

我们将数字number存入y中准备对其进行位操作，但是众所周知浮点数没有位操作的语句。
但是长整型数可以进行位操作！

如果执行下面这一行代码：

i=long(y);
y=4.456时i=4

但是我们想要的是将y逐位转化成长整型，换句话说按照上面的代码我们改变了二进制逐位的值。
因此我们采取下面代码：

i= (long * ) &y;

这句话干了什么呢？

简单来说它转换了y对应的内存地址而不是y这个数字本身。
或者说这句话告诉cpu把y指向的内存数据当成长整形来读取

&y              取到浮点数y的地址
(long * ) &y    将该地址转化成长整型
*(long * ) &y    再取其值

6666指针确实是个神奇的东西，我们继续向下
通过这句话浮点数y的二进制表示存入到了i中。
②

i= 0x5f3759df-(i>>1);

这句的理解在第三大部分可以轻松找到。
简单来说我们通过放缩系数和平移操作代替了除法和开方的过程，这大大提高了代码运行效率。
③

y = * (float *) &i;

和①类似我们再次通过这种操作将二进制位转换成浮点数。

通过上述三步我们得到了一个近似值。因此我们还需要进行牛顿迭代法继续使结果精确。
④

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

下面进行牛顿迭代法的数学推导：
y = 1 x y= \frac{1}{\sqrt x} y=x 1
f ( y ) = 1 y 2 − x f(y)=\frac{1}{y^2}-x f(y)=y21−x
y n + 1 = y n − 1 y n 2 − x − 2 y n 3 = y n ∗ ( 3 2 − 1 2 ∗ x ∗ y n 2 ) y_{n+1}=y_n-\frac{\frac{1}{y_n^2}-x}{\frac{-2}{y_n^3}}=y_n*(\frac{3}{2}-\frac{1}{2}*x*y_n^2) yn+1=yn−yn3−2yn21−x=yn∗(23−21∗x∗yn2)

五.总结
短短的几行代码背后是强大的数学支撑，万万想不到之前学过的高数会在算法中这样体现。
第一次试写涉及到数学公式推导的文章，如有疑问请大家多多留言！

快速平方根倒数算法深度理解相关推荐

MATLAB可视化实战系列（二十八）-贪心算法求快速平方根倒数算法中的“魔术数字”【含matlab源代码】
前言快速平方根倒数算法(Fast InvSqrt)是一种快速计算平方根的倒数的算法,常用于向量标准化运算,在光照渲染中有重要应用.此算法最早可能是于90年代前期由SGI所发明,后来于1999年在&l ...
均方根误差不超过_快速平方根倒数算法
论文地址戳这里www.lomont.org 一. 介绍快速平方根倒数算法也称为平方根倒数速算法(Fast Inverse Square Root)是用于快速计算的一种算法.此算法由于出现在< ...
雷神之锤3快速计算算术平方根倒数算法中魔法数字的另一种求法（1）
对于雷神之锤3中快速计算算术平方根倒数算法中魔法数字的真正来源一直是个悬案. 本人对此进行了一番研究,有幸参悟其中奥秘,特分享给大家. 本人不爱码字,直接上图上图中代码出自雷神之锤3,本人对其中魔法 ...
《雷神之锤III》平方根倒数算法学习笔记
今天刷到个算法视频,觉得很有意思,所以打算把它记录下来. 平方根倒数算法源代码二进制浮点数运算对数技巧平方根倒数牛顿迭代源代码 float Q_rsqrt(float number) {l ...
【转】卡马克快速平方根——平方根倒数算法
-------------------------------------------------------------------------------- 快速平方根(平方根倒数)算法日前在书 ...
卡马克快速平方根倒数
在1999 年<雷神之锤III竞技场>源代码中出现的这个代码而被用户所熟知牛顿迭代法(Newton-Raphson Method,简称NR) 而牛顿迭代法的基础则是泰勒级数(Taylor ...
快速开平方根倒数算法(Fast inverse square root)的一点探究
文章目录一.写在前面 1. 提示 2. 背景与前情二.正文 1. 需求分析 2. 必备工具之IEEE-754浮点数表示方法 3. 同一储存单元32bits的两种不同意义 4. 公式推导 4. 本文 ...
《雷神之锤 Ⅲ》平方根倒数速算法魔术数字的另一种求法(2)
上图代码是出自雷神之锤3计算一个正实数的算术平方根的倒数的算法,一般我们去计算一个数的平方根采用最好的方法是牛顿迭代法不会牛顿迭代法,可以去看下这位大哥的帖子 www.cnblogs.com/hou ...
【华为云技术分享】深度理解AI概念、算法及如何进行AI项目开发
莫衷一是的AI 做了多年的业务工作,一直希望能够用机器代替人力,把人从繁琐的具体工作中解放出来,从技术发展看AI或许可以支撑实现这个愿景. 但最近关于AI的讨论和争论比较多,学术上,纽约大学的Gary ...

快速平方根倒数算法深度理解

快速平方根倒数算法深度理解

快速平方根倒数算法深度理解相关推荐

最新文章

热门文章