平方根倒数速算法是适用于快速计算(积的平方根的倒数,在此需取符合IEEE 754标準格式的32位浮点数)的一种算法。
平方根倒数速算法(英语:Fast Inverse Square Root,亦常以“Fast InvSqrt()”或其使用的十六进制常数0x5f3759df代称)是用于快速计算(积的平方根的倒数,在此需取符合IEEE 754标準格式的32位浮点数)的一种算法。此算法最早可能是于90年代前期由SGI所发明,后来则于1999年在《雷神之锤III竞技场》的原始码中套用,但直到2002-2003年间才在Usenet一类的公共论坛上出现。这一算法的优势在于减少了求平方根倒数时浮点运算操作带来的巨大的运算耗费,而在计算机图形学领域,若要求取照明和投影的波动角度与反射效果,就常需计算平方根倒数。
此算法首先接收一个32位带符浮点数,然后将之作为一个32位整数看待,以将其向右进行一次逻辑移位的方式将之取半,并用十六进制“魔术数字”0x5f3759df减之,如此即可得对输入的浮点数的平方根倒数的首次近似值;而后重新将其作为浮点数,以牛顿法反覆叠代,以求出更精确的近似值,直至求出符合精确度要求的近似值。在计算浮点数的平方根倒数的同一精度的近似值时,此算法比直接使用浮点数除法要快四倍。
基本介绍
- 中文名:平方根倒数速算法
- 外文名:Fast Inverse Square Root
- 别名:Fast InvSqrt(),0x5f3759df
- 定义:适用于快速计算的一种算法
- 备注:最早被认为由约翰·卡马克所发明
算法定义
平方根倒数速算法最早被认为是由约翰·卡马克所发明,但后来的调查显示,该算法在这之前就于计算机图形学的硬体与软体领域有所套用,如SGI和3dfx就曾在产品中套用此算法。而就现在所知,此算法最早由Gary Tarolli在SGI Indigo的开发中使用。虽说在随后的相关研究中也提出了一些可能的来源,但至今为止仍未能确切知晓此常数的起源。
算法切入
浮点数的平方根倒数常用于计算正规化矢量。3D图形程式需要使用正规化矢量来实现光照和投影效果,因此每秒都需做上百万次平方根倒数运算,而在处理坐标转换与光源的专用硬体设备出现前,这些计算都由软体完成,计算速度亦相当之慢;在1990年代这段代码开发出来之时,多数浮点数操作的速度更是远远滞后于整数操作,因而针对正规化矢量算法的最佳化就显得尤为重要。下面陈述计算正规化矢量的原理:
要将一个矢量标準化,就必须计算其欧几里得範数以求得矢量长度,而这时就需对矢量的各分量的平方和求平方根;而当求取到其长度并以之除该矢量的每个分量后,所得的新矢量就是与原矢量同向的单位矢量,若以公式表示:
- 可求得矢量v的欧几里得範数,此算法正类如对欧几里得空间的两点求取其欧几里得距离,
- 而求得的就是标準化的矢量,若以代表,则有,
可见标準化矢量时需要用到对矢量分量的平方根倒数计算,所以对平方根倒数计算算法的最佳化对计算正规化矢量也大有裨益。
为了加速图像处理单元计算,《雷神之锤III竞技场》使用了平方根倒数速算法,而后来採用现场可程式逻辑门阵列的顶点着色器也套用了此算法。
浮点化整
要理解这段代码,首先需了解浮点数的存储格式。一个浮点数以32个二进制位表示一个有理数,而这32位由其意义分为三段:首先首位为符号位,如若是0则为正数,反之为负数;接下来的8位表示经过偏移处理(这是为了使之能表示-127-128)后的指数;最后23位表示的则是有效数字中除最高位以外的其余数字。将上述结构表示成公式即为,其中表示有效数字的尾数(此处,偏移量,而指数的值决定了有效数字(在Lomont和McEniry的论文中称为“尾数”(mantissa))代表的是小数还是整数。以上图为例,将描述带入有),且,则可得其表示的浮点数为。
符号位 | |||||||||
0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | = | 127 |
0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | = | 2 |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | = | 1 |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | = | 0 |
1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | = | −1 |
1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | = | −2 |
1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | = | −127 |
1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | = | −128 |
8位二进制整数补码示例
如上所述,一个有符号正整数在二进制补码系统中的表示中首位为0,而后面的各位则用于表示其数值。将浮点数取别名存储为整数时,该整数的数值即为,其中E表示指数,M表示有效数字;若以上图为例,图中样例若作为浮点数看待有,,则易知其转化而得的整数型号数值为。由于平方根倒数函式仅能处理正数,因此浮点数的符号位(即如上的Si)必为0,而这就保证了转换所得的有符号整数也必为正数。以上转换就为后面的计算带来了可行性,之后的第一步操作(逻辑右移一位)即是使该数的长整形式被2所除。
历史考究

《雷神之锤III》的代码直到QuakeCon 2005才正式放出,但早在2002年(或2003年)时平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了。最初人们猜测是卡马克写下了这段代码,但他在询问邮件的回覆中否定了这个观点,并猜测可能是先前曾帮id Software最佳化雷神之锤的资深彙编程式设计师Terje Mathisen写下了这段代码;而在Mathisen的邮件里他表示在1990年代初他只曾作过类似的实现,确切来说这段代码亦非他所作。现在所知的最早实现是由Gary Tarilli在SGI Indigo中实现的,但他亦坦承他仅对常数R的取值做了一定的改进,实际上他也不是作者。Rys Sommefeldt则在向以发明MATLAB而闻名的Cleve Moler查证后认为原始的算法是Ardent Computer公司的Greg Walsh所发明,但他也没有任何决定性的证据能证明这一点。
目前不仅该算法的原作者不明,人们也仍无法明确当初选择这个“魔术数字”的方法。Chris Lomont在研究中曾做了个试验:他编写了一个函式,以在一个範围内遍历选取R值的方式将逼近误差降到最小,以此方法他计算出了线性近似的最优R值0x5f37642f(与代码中使用的0x5f3759df相当接近),但以之代入算法计算并进行一次牛顿叠代后,所得近似值与代入0x5f3759df的结果相比精度却仍略微更低;而后Lomont将目标改为遍历选取在进行1-2次牛顿叠代后能得到最大精度的R值,并由此算出最优R值为0x5f375a86,以此值代入算法并进行牛顿叠代后所得的结果都比代入原始值(0x5f3759df)更精确,于是他的论文最后以“原始常数是以数学推导还是以反覆试错的方式求得”的问题作结。在论文中Lomont亦指出64位的IEEE754浮点数(即双精度类型)所对应的魔术数字是0x5fe6ec85e7de30da,但后来的研究表明代入0x5fe6eb50c7aa19f9的结果精确度更高(McEniry得出的结果则是0x5FE6EB50C7B537AA,精度介于两者之间)。在Charles McEniry的论文中,他使用了一种类似Lomont但更複杂的方法来最佳化R值:他最开始使用穷举搜寻法,所得结果与Lomont相同;而后他尝试用带权二分法寻找最优值,所得结果恰是代码中所使用的魔术数字0x5f3759df,因此McEniry确信这一常数或许最初便是以“在可容忍误差範围内使用二分法”的方式求得。
注释信息
- 由于现代计算机系统对长整型的定义有所差异,使用长整型会降低此段代码的可移植性。具体来说,由此段浮点转换为长整型的定义可知,如若这段代码正常运行,所在系统的长整型长度应为4位元组(32位),否则重新转为浮点数时可能会变成负数;而由于C99标準的广泛套用,在现今多数64位计算机系统(除使用LLP64数据模型的Windows外)中,长整型的长度都是8位元组。
- 此处“浮点数”所指为标準化浮点数,也即有效数字部分必须满足,可参见David Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys. 1991.March, 23 (1): 5–48. doi:10.1145/103162.103163.
- Lomont 2003确定R的方式则有所不同,他先将R分解为与,其中与分别代表R的有效数字域和指数域。