fbpx
维基百科

平方根倒数速算法

平方根倒数速算法(英語:Fast Inverse Square Root,亦常以“Fast InvSqrt()”或其使用的十六进制常数0x5f3759df代称)是用于快速计算(即平方根倒数,在此需取符合IEEE 754标准格式的32位浮点数)的一种算法。这一算法的优势在于减少了求平方根倒数时浮点运算操作带来的巨大的运算耗费,而在计算机图形学领域,若要求取照明投影波动角度反射效果,就常需计算平方根倒数。

游戏实现光照和反射效果时以平方根倒数速算法计算波动角度,此图以第一人称射击游戏OpenArena为例。

此算法首先接收一个32位带符浮点数,然后将之作为一个32位整数看待,以将其向右进行一次逻辑移位的方式将之取半,并用在浮點數規格代表近似值的十六进制“魔术数字”0x5f3759df减之,如此即可得对输入的浮点数的平方根倒数的首次近似值;而后重新将其作为浮点数,以牛顿法反复迭代,以求出更精确的近似值,直至求出符合精确度要求的近似值。在计算浮点数的平方根倒数的同一精度的近似值时,此算法比直接使用浮点数除法要快四倍。

此算法最早被认为是由约翰·卡马克于90年代前期在SGI Indigo英语SGI Indigo的开发中使用,后来则于1999年在《雷神之锤III竞技场》的源代码中应用,但直到2002-2003年间才在Usenet一类的公共论坛上出现[1]。后来的调查显示,该算法在这之前就于计算机图形学的硬件与软件领域有所应用,如SGI3dfx就曾在产品中应用此算法,但至今为止仍未能确切知晓算法中所使用的特殊常数的起源。

算法的切入点 编辑

 
法线常在光影效果实现计算时使用,而这就涉及到向量范数的计算。图中所标识的就是与一个面所垂直的一些向量的集合。

浮点数的平方根倒数常用于计算正规化向量[文 1]。3D图形程序需要使用正规化向量来实现光照和投影效果,因此每秒都需做上百万次平方根倒数运算,而在处理坐标转换与光源的专用硬件设备出现前,这些计算都由软件完成,计算速度亦相当之慢;在1990年代这段代码开发出来之时,多数浮点数操作的速度更是远远滞后于整数操作[1],因而针对正规化向量算法的优化就显得尤为重要。下面陈述计算正规化向量的原理:

要将一个向量标准化,就必须计算其欧几里得范数,以求得向量长度,为此便需对向量的各分量的平方和求平方根;而当求取到其长度,并以之除该向量的每个分量后,所得的新向量就是与原向量同向的单位向量,若以公式表示:

 可求得向量v的欧几里得范数,此算法正类如对欧几里得空间的两点求取其欧几里得距离
 求得的就是标准化的向量,若以 代表 ,则有 

可见标准化向量时,对向量分量计算平方根倒数实为必需,所以,对平方根倒数计算算法的优化对计算正规化向量也大有裨益。

为了加速图像处理单元计算,《雷神之锤III竞技场》使用了平方根倒数速算法,而后来采用现场可编程逻辑门阵列顶点着色器也应用了此算法[文 2]

代码概览 编辑

下列代码是《雷神之锤III竞技场》源代码中平方根倒数速算法之實例。示例去除了C预处理器的指令,但附上了原有的注释(括号内为注释的翻译)[2]

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; } 

为计算平方根倒数的值,软件首先要先确定一个近似值,而后则使用某些数值方法不断计算修改近似值,直至达到可接受的精度。在1990年代初(也即该算法发明的大概时间),软件开发时通用的平方根演算法英语Methods of computing square roots多是从查找表中取得近似值[文 3],而这段代码取近似值耗时比之更短,达到精确度要求的速度也比通常使用的浮点除法计算法快四倍[文 4],虽然此算法会损失一些精度,但性能上的巨大优势已足以补偿损失的精度[文 5]。由代码中对原数据的变量类型声明为float可看出,这一算法是针对IEEE 754标准格式的32位浮点数设计的,不过据Chris Lomont和后来的Charles McEniry的研究看,这一算法亦可套用于其他类型的浮点数上。

平方根倒数速算法在速度上的优势源自将浮点数转化为长整型[註 1]以作整数看待,并用特定常数0x5f3759df与之相减。然而对于代码阅读者来说,他们却难以立即领悟出使用这一常数的目的,因此和其它在代码中出现的难以理解的常数一样,这一常数亦被称为“魔术数字[1][文 7][文 8][文 9]。如此将浮点数当作整数先位移后减法,所得的浮点数结果即是对输入数字的平方根倒数的粗略估计值,而后再进行一次牛顿迭代法,以使之更精确后,代码即执行完毕。由于算法所生成的用于输入牛顿法的首次近似值已经相当精确,此算法所得近似值的精度已可接受,而若使用与《雷神之锤III竞技场》同为1999年发布的Pentium III中的SSE指令rsqrtss计算,则计算平方根倒数的收敛速度更慢,精度也更低[4]

将浮点数转化为整数 编辑

 

要理解这段代码,首先需了解浮点数的存储格式。一个浮点数以32个二进制位表示一个有理数,而这32位由其意义分为三段:首先首位为符号位,如若是0则为正数,反之为负数;接下来的8位表示经过偏移处理(这是为了使之能表示-127-128)后的指数;最后23位表示的则是有效数字中除最高位以外的其余数字。将上述结构表示成公式即为 ,其中 表示有效数字的尾数(此处 ,偏移量 [文 10],而指数的值 决定了有效数字(在Lomont和McEniry的论文中称为“尾数”(mantissa))代表的是小数还是整数[文 11]。以上图为例,将描述带入有 ),且 ,则可得其表示的浮点数为 

符号位
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,而后面的各位则用于表示其数值。将浮点数取别名英语Aliasing (computing)存储为整数时,该整数的数值即为 ,其中E表示指数,M表示有效数字;若以上图为例,图中样例若作为浮点数看待有  ,则易知其转化而得的整数型号数值为 。由于平方根倒数函数仅能处理正数,因此浮点数的符号位(即如上的Si)必为0,而这就保证了转换所得的有符号整数也必为正数。以上转换就为后面的计算带来了可行性,之后的第一步操作(逻辑右移一位)即是使该数的长整形式被2所除[文 12]

“魔术数字” 编辑

S(ign,符号) E(xponent,指数) M(antissa,尾数)
1 位 b位 (n-1-b)
n位[文 13]

对猜测平方根倒数速算法的最初构想来说,计算首次近似值所使用的常数0x5f3759df也是重要的线索。为确定程序员最初选此常数以近似求取平方根倒数的方法,Charles McEniry首先检验了在代码中选择任意常数R所求取出的首次近似值的精度。回想上一节关于整数和浮点数表示的比较:对于同样的32位二进制数码,若为浮点数表示时实际数值为 ,而若为整数表示时实际数值则为 [註 2],其中 。以下等式引入了一些由指数和有效数字导出的新元素:

 
 ,其中 

再继续看McEniry 2007里的进一步说明:

 

对等式的两边取二进制对数 ,即函数 反函数),有

 
 

以如上方法,就能将浮点数x和y的相关指数消去,从而将乘方运算化为加法运算。而由于  线性相关,因此在  (即输入值与首次近似值)间就可以线性组合的方式建立方程[文 10]。在此McEniry再度引入新数 描述 与近似值R间的误差[註 3]:由于 ,有 ,则在此可定义 与x的关系为 ,这一定义就能提供二进制对数的首次精度值(此处 ;当R为0x5f3759df时,有 [文 13])。由此将 代入上式,有:

 

参照首段等式代入    ,有:

 

移项整理得:

 

如上所述,对于以浮点规格存储的正浮点数x,若将其作为长整型表示则示值为 ,由此即可根据x的整数表示导出y(在此 ,亦即x的平方根倒数的首次近似值)的整数表示值,也即:

 ,其中 .

最后导出的等式 即与上节代码中i = 0x5f3759df - (i>>1);一行相契合,由此可见,在平方根倒数速算法中,对浮点数进行一次移位操作与整数减法,就可以可靠地输出一个浮点数的对应近似值[文 13]。到此为止,McEniry只证明了,在常数R的辅助下,可近似求取浮点数的平方根倒数,但仍未能确定代码中的R值的选取方法。

关于作一次移位与减法操作以使浮点数的指数被-2除的原理,Chris Lomont的论文中亦有个相对简单的解释:以 为例,将其指数除-2可得 ;而由于浮点表示的指数有进行过偏移处理,所以指数的真实值e应为 ,因此可知除法操作的实际结果为 [文 14],这时用R(在此即为“魔术数字”0x5f3759df)减之即可使指数的最低有效数位转入有效数字域,之后重新转换为浮点数时,就能得到相当精确的平方根倒数近似值。在这里对常数R的选取亦有所讲究,若能选取一个好的R值,便可减少对指数进行除法与对有效数字域进行移位时可能产生的错误。基于这一标准,0xbe即是最合适的R值,而0xbe右移一位即可得到0x5f,这恰是魔术数字R的第一个字节[文 15]

精确度 编辑

 
使用启发式平方根倒数速算法与使用C语言标准库libstdc的函数所计算出的平方根倒数的差值一览,注意这里使用的是双对数坐标系

如上所述,平方根倒数速算法所得的近似值惊人的精确,右图亦展示了以上述代码计算(以平方根倒数速算法计算后再进行一次牛顿法迭代)所得近似值的误差:当输入0.01时,以C语言标准库函数计算可得10.0,而InvSqrt()得值为9.9825822,其间误差为0.017479,相对误差则为0.175%,且当输入更大的数值时,绝对误差不断下降,相对误差也一直控制在一定的范围之内。

牛顿法提高精度 编辑

在进行了如上的整数操作之后,示例程序再度将被转为长整型的浮点数回转为浮点数(对应x = *(float*)&i;),并对其进行一次浮点运算操作(对应x = x*(1.5f - xhalf*x*x);),这里的浮点运算操作就是对其进行一次牛顿法迭代,若以此例说明:

 所求的是y的平方根倒数,以之构造以y为自变量的函数,有 
将其代入牛顿法的通用公式 (其中 为首次近似值),
 ,其中  
整理有 ,对应的代码即为y = y*(threehalfs - xhalf*y*y);

在以上一节的整数操作产生首次近似值后,程序会将首次近似值作为参数送入函数最后两句进行精化处理,代码中的两次迭代(以一次迭代的输出(对应公式中的 )作为二次迭代的输入)正是为了进一步提高结果的精度[文 16],但由于雷神之锤III引擎的图形计算中并不需要太高的精度,所以代码中只进行了一次迭代,二次迭代的代码则被注释[文 9]

历史与考究 编辑

 
id Software的创始人约翰·卡马克。这段代码虽非他所作,但常被认为与他相关。

《雷神之锤III》的代码直到QuakeCon 2005才正式放出,但早在2002年(或2003年)时,平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了[1]。最初人们猜测是卡马克写下了这段代码,但他在回复询问他的邮件时否定了这个观点,并猜测可能是先前曾帮id Software优化雷神之锤的资深汇编程序员Terje Mathisen写下了这段代码;而在Mathisen的邮件裡,他表示,在1990年代初,他只曾作过类似的實作,确切来说这段代码亦非他所作。现在所知的最早實作是由Gary Tarolli在SGI Indigo中實作的,但他亦坦承他仅对常数R的取值做了一定的改进,实际上他也不是作者。在向以发明MATLAB而闻名的Cleve Moler查证后,Rys Sommefeldt则认为原始的算法是Ardent Computer英语Ardent Computer公司的Greg Walsh所发明,但他也没有任何决定性的证据能证明这一点[5]

不仅该算法的原作者不明,人们也仍无法確定当初选择这个“魔术数字”的方法。Chris Lomont曾做了个研究:他推算出了一个函数以討論此速算法的誤差,並找出了使誤差最小的最佳R值0x5f37642f(与代码中使用的0x5f3759df相当接近),但以之代入算法计算并进行一次牛顿迭代后,所得近似值之精度仍略低於代入0x5f3759df的结果[文 17];因此Lomont将目标改为尋找在进行1-2次牛顿迭代后能得到最大精度的R值,在暴力搜尋後得出最优R值为0x5f375a86,以此值代入算法并进行牛顿迭代,所得的结果都比代入原始值(0x5f3759df)更精确[文 17],于是他的论文最后提到“如果可能我想询问原作者,此速算法是以数学推导还是以反复试错的方式求得?”[文 18]。在论文中,Lomont亦指出,64位的IEEE754浮点数(即双精度类型)所对应的魔术数字是0x5fe6ec85e7de30da,但后来的研究表明,代入0x5fe6eb50c7aa19f9的结果精确度更高(McEniry得出的结果则是0x5FE6EB50C7B537AA,精度介于两者之间,Matthew Robertson给出的精度更高的值是0x5FE6EB50C7B537A9[6])。在Charles McEniry的论文中,他使用了一种类似Lomont但更复杂的方法来优化R值:他最开始使用穷举搜索,所得结果与Lomont相同[文 19];而后他尝试用带权二分法寻找最优值,所得结果恰是代码中所使用的魔术数字0x5f3759df,因此,McEniry认为,这一常数最初或许便是以“在可容忍误差范围内使用二分法”的方式求得[文 20]

注释 编辑

  1. ^ 由于现代计算机系统对长整型的定义有所差异,使用长整型会降低此段代码的可移植性。具体来说,由此段浮点转换为长整型的定义可知,如若这段代码正常运行,所在系统的长整型长度应为4字节(32位),否则重新转为浮点数时可能会变成负数;而由于C99标准的广泛应用,在现今多数64位计算机系统(除使用LLP64数据模型的Windows外)中,长整型的长度都是8字节[文 6][3]
  2. ^ 此处“浮点数”所指为标准化浮点数,也即有效数字部分必须满足 ,可参见David Goldberg. What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys. March 1991, 23 (1): 5–48. doi:10.1145/103162.103163. 
  3. ^ Lomont 2003确定R的方式则有所不同,他先将R分解为  ,其中  分别代表R的有效数字域和指数域[文 7]

参考 编辑

脚注 编辑

  1. ^ 1.0 1.1 1.2 1.3 Sommefeldt, Rys. . Beyond3D. 2006-11-29 [2009-02-12]. (原始内容存档于2009-02-09) (英语). 
  2. ^ quake3-1.32b/code/game/q_math.c. Quake III Arena. id Software. [2010-11-15]. (原始内容存档于2011-02-17). 
  3. ^ Meyers, Randy. . drdobbs.com. 2000-12-01 [2010-09-04]. (原始内容存档于2012-05-27). 
  4. ^ Ruskin, Elan. . Some Assembly Required. 2009-10-16 [2010-09-13]. (原始内容存档于2010-08-25). 
  5. ^ Sommefeldt, Rys. . Beyond3D. 2006-12-19 [2008-04-19]. (原始内容存档于2008-05-13). 
  6. ^ Matthew Robertson. A Brief History of InvSqrt (PDF). UNBSJ. 2012-04-24 [2023-11-23]. (原始内容 (PDF)于2023-01-29). 

参考文献 编辑

  • Blinn, Jim. Jim Blinn's Corner: Notation, notation notation. Morgan Kaufmann. 2003. ISBN 1-55860-860-5. 
  • Eberly, David. 3D Game Engine Design. Morgan Kaufmann. 2001. ISBN 978-1-55860-593-0. 
  • Eberly, David. (PDF). Geometric Tools. 2002 [2009-03-22]. (原始内容 (PDF)存档于2009-02-24). 
  • Fog, Agner. (PDF). 2010-02-16 [2010-08-30]. (原始内容 (PDF)存档于2010-08-15). 
  • Hennessey, John; Patterson, David A. Computer Organization and Design 2nd. San Francisco, CA: Morgan Kaufmann Publishers. 1998. ISBN 978-1-55860-491-9. 
  • Lomont, Chris. (PDF). February 2003 [2009-02-13]. (原始内容 (PDF)存档于2009-02-06). 
  • McEniry, Charles. (PDF). August 2007 [2009-02-13]. (原始内容 (PDF)存档于2015-05-11). 
  • Middendorf, Lars; Mühlbauer, Felix; Umlauf, George; Bodba, Christophe. Embedded Vertex Shader in FPGA. Rettberg, Achin (编). Embedded System Design: Topics, Techniques and Trends. IFIP TC10 Working Conference:International Embedded Systems Symposium (IESS). et al. Irvine, California: Springer. 2007-06-01. ISBN 978-0-387-72257-3. 

延伸阅读 编辑

  • Kushner, David. The wizardry of Id. IEEE Spectrum. Aug 2002, 39 (8): 42–47. doi:10.1109/MSPEC.2002.1021943. 
  • 顾森(Matrix67). 旧闻一则:神秘的0x5f3759df 不可思议的Quake III源码. 2007-11-24 [2015-08-01]. (原始内容于2015-08-04). 

外部链接 编辑

  • Origin of Quake3's Fast InvSqrt() (页面存档备份,存于互联网档案馆), Beyond3D.com
  • Quake III Arena source code (页面存档备份,存于互联网档案馆), id Software
  • Margolin, Tomer. . CodeMaestro. The Coding Experience. 2005-08-27. (原始内容存档于2012-04-14). 

平方根倒数速算法, 英語, fast, inverse, square, root, 亦常以, fast, invsqrt, 或其使用的十六进制常数0x5f3759df代称, 是用于快速计算x, displaystyle, textstyle, 即x, displaystyle, textstyle, 的平方根的倒数, 在此x, displaystyle, textstyle, 需取符合ieee, 754标准格式的32位浮点数, 的一种算法, 这一算法的优势在于减少了求平方根倒数时浮点运算操作带来的巨大的运算耗费,. 平方根倒数速算法 英語 Fast Inverse Square Root 亦常以 Fast InvSqrt 或其使用的十六进制常数0x5f3759df代称 是用于快速计算x 1 2 displaystyle textstyle x 1 2 即x displaystyle textstyle x 的平方根的倒数 在此x displaystyle textstyle x 需取符合IEEE 754标准格式的32位浮点数 的一种算法 这一算法的优势在于减少了求平方根倒数时浮点运算操作带来的巨大的运算耗费 而在计算机图形学领域 若要求取照明和投影的波动角度与反射效果 就常需计算平方根倒数 游戏实现光照和反射效果时以平方根倒数速算法计算波动角度 此图以第一人称射击游戏OpenArena为例 此算法首先接收一个32位带符浮点数 然后将之作为一个32位整数看待 以将其向右进行一次逻辑移位的方式将之取半 并用在浮點數規格代表2127 displaystyle textstyle sqrt 2 127 近似值的十六进制 魔术数字 0x5f3759df减之 如此即可得对输入的浮点数的平方根倒数的首次近似值 而后重新将其作为浮点数 以牛顿法反复迭代 以求出更精确的近似值 直至求出符合精确度要求的近似值 在计算浮点数的平方根倒数的同一精度的近似值时 此算法比直接使用浮点数除法要快四倍 此算法最早被认为是由约翰 卡马克于90年代前期在SGI Indigo 英语 SGI Indigo 的开发中使用 后来则于1999年在 雷神之锤III竞技场 的源代码中应用 但直到2002 2003年间才在Usenet一类的公共论坛上出现 1 后来的调查显示 该算法在这之前就于计算机图形学的硬件与软件领域有所应用 如SGI和3dfx就曾在产品中应用此算法 但至今为止仍未能确切知晓算法中所使用的特殊常数的起源 目录 1 算法的切入点 2 代码概览 2 1 将浮点数转化为整数 2 2 魔术数字 2 3 精确度 2 4 牛顿法提高精度 3 历史与考究 4 注释 5 参考 5 1 脚注 5 2 参考文献 6 延伸阅读 7 外部链接算法的切入点 编辑 nbsp 法线常在光影效果实现计算时使用 而这就涉及到向量范数的计算 图中所标识的就是与一个面所垂直的一些向量的集合 浮点数的平方根倒数常用于计算正规化向量 文 1 3D图形程序需要使用正规化向量来实现光照和投影效果 因此每秒都需做上百万次平方根倒数运算 而在处理坐标转换与光源的专用硬件设备出现前 这些计算都由软件完成 计算速度亦相当之慢 在1990年代这段代码开发出来之时 多数浮点数操作的速度更是远远滞后于整数操作 1 因而针对正规化向量算法的优化就显得尤为重要 下面陈述计算正规化向量的原理 要将一个向量标准化 就必须计算其欧几里得范数 以求得向量长度 为此便需对向量的各分量的平方和求平方根 而当求取到其长度 并以之除该向量的每个分量后 所得的新向量就是与原向量同向的单位向量 若以公式表示 v v12 v22 v32 displaystyle boldsymbol v sqrt v 1 2 v 2 2 v 3 2 nbsp 可求得向量v的欧几里得范数 此算法正类如对欧几里得空间的两点求取其欧几里得距离 而v v v displaystyle boldsymbol hat v boldsymbol v boldsymbol v nbsp 求得的就是标准化的向量 若以x displaystyle boldsymbol x nbsp 代表v12 v22 v32 displaystyle v 1 2 v 2 2 v 3 2 nbsp 则有v v x displaystyle boldsymbol hat v boldsymbol v sqrt x nbsp 可见标准化向量时 对向量分量计算平方根倒数实为必需 所以 对平方根倒数计算算法的优化对计算正规化向量也大有裨益 为了加速图像处理单元计算 雷神之锤III竞技场 使用了平方根倒数速算法 而后来采用现场可编程逻辑门阵列的顶点着色器也应用了此算法 文 2 代码概览 编辑下列代码是 雷神之锤III竞技场 源代码中平方根倒数速算法之實例 示例去除了C预处理器的指令 但附上了原有的注释 括号内为注释的翻译 2 float Q rsqrt float number long i float x2 y const float threehalfs 1 5F x2 number 0 5F y number i long amp y evil floating point bit level hacking 邪恶的浮点数位运算黑科技 i 0x5f3759df i gt gt 1 what the fuck 这是什么鬼 y float amp i y y threehalfs x2 y y 1st iteration 第一次迭代 y y threehalfs x2 y y 2nd iteration this can be removed 第二次迭代 可以删除 return y 为计算平方根倒数的值 软件首先要先确定一个近似值 而后则使用某些数值方法不断计算修改近似值 直至达到可接受的精度 在1990年代初 也即该算法发明的大概时间 软件开发时通用的平方根演算法 英语 Methods of computing square roots 多是从查找表中取得近似值 文 3 而这段代码取近似值耗时比之更短 达到精确度要求的速度也比通常使用的浮点除法计算法快四倍 文 4 虽然此算法会损失一些精度 但性能上的巨大优势已足以补偿损失的精度 文 5 由代码中对原数据的变量类型声明为float可看出 这一算法是针对IEEE 754标准格式的32位浮点数设计的 不过据Chris Lomont和后来的Charles McEniry的研究看 这一算法亦可套用于其他类型的浮点数上 平方根倒数速算法在速度上的优势源自将浮点数转化为长整型 註 1 以作整数看待 并用特定常数0x5f3759df与之相减 然而对于代码阅读者来说 他们却难以立即领悟出使用这一常数的目的 因此和其它在代码中出现的难以理解的常数一样 这一常数亦被称为 魔术数字 1 文 7 文 8 文 9 如此将浮点数当作整数先位移后减法 所得的浮点数结果即是对输入数字的平方根倒数的粗略估计值 而后再进行一次牛顿迭代法 以使之更精确后 代码即执行完毕 由于算法所生成的用于输入牛顿法的首次近似值已经相当精确 此算法所得近似值的精度已可接受 而若使用与 雷神之锤III竞技场 同为1999年发布的Pentium III中的SSE指令rsqrtss计算 则计算平方根倒数的收敛速度更慢 精度也更低 4 将浮点数转化为整数 编辑 nbsp 参见 浮点数 要理解这段代码 首先需了解浮点数的存储格式 一个浮点数以32个二进制位表示一个有理数 而这32位由其意义分为三段 首先首位为符号位 如若是0则为正数 反之为负数 接下来的8位表示经过偏移处理 这是为了使之能表示 127 128 后的指数 最后23位表示的则是有效数字中除最高位以外的其余数字 将上述结构表示成公式即为x 1 Si 1 m 2 E B displaystyle textstyle x 1 mathrm Si cdot 1 m cdot 2 E B nbsp 其中m displaystyle textstyle m nbsp 表示有效数字的尾数 此处0 m lt 1 displaystyle textstyle 0 leq m lt 1 nbsp 偏移量B 127 displaystyle textstyle B 127 nbsp 文 10 而指数的值E B displaystyle textstyle E B nbsp 决定了有效数字 在Lomont和McEniry的论文中称为 尾数 mantissa 代表的是小数还是整数 文 11 以上图为例 将描述带入有m 1 2 2 0 250 displaystyle textstyle m 1 times 2 2 0 250 nbsp 且E B 124 127 3 displaystyle textstyle E B 124 127 3 nbsp 则可得其表示的浮点数为x 1 0 250 2 3 0 15625 displaystyle textstyle x 1 0 250 cdot 2 3 0 15625 nbsp 符号位0 1 1 1 1 1 1 1 1270 0 0 0 0 0 1 0 20 0 0 0 0 0 0 1 10 0 0 0 0 0 0 0 01 1 1 1 1 1 1 1 11 1 1 1 1 1 1 0 21 0 0 0 0 0 0 1 1271 0 0 0 0 0 0 0 1288位二进制补码示例 如上所述 一个有符号正整数在二进制补码系统中的表示中首位为0 而后面的各位则用于表示其数值 将浮点数取别名 英语 Aliasing computing 存储为整数时 该整数的数值即为I E 223 M displaystyle textstyle I E times 2 23 M nbsp 其中E表示指数 M表示有效数字 若以上图为例 图中样例若作为浮点数看待有E 124 displaystyle textstyle E 124 nbsp M 1 221 displaystyle M 1 cdot 2 21 nbsp 则易知其转化而得的整数型号数值为I 124 223 221 displaystyle I 124 times 2 23 2 21 nbsp 由于平方根倒数函数仅能处理正数 因此浮点数的符号位 即如上的Si 必为0 而这就保证了转换所得的有符号整数也必为正数 以上转换就为后面的计算带来了可行性 之后的第一步操作 逻辑右移一位 即是使该数的长整形式被2所除 文 12 魔术数字 编辑 S ign 符号 E xponent 指数 M antissa 尾数 1 位 b位 n 1 b 位n位 文 13 对猜测平方根倒数速算法的最初构想来说 计算首次近似值所使用的常数0x5f3759df也是重要的线索 为确定程序员最初选此常数以近似求取平方根倒数的方法 Charles McEniry首先检验了在代码中选择任意常数R所求取出的首次近似值的精度 回想上一节关于整数和浮点数表示的比较 对于同样的32位二进制数码 若为浮点数表示时实际数值为x 1 mx 2ex displaystyle textstyle x 1 m x 2 e x nbsp 而若为整数表示时实际数值则为Ix ExL Mx displaystyle textstyle I x E x L M x nbsp 註 2 其中L 2n 1 b displaystyle textstyle L 2 n 1 b nbsp 以下等式引入了一些由指数和有效数字导出的新元素 mx MxL displaystyle m x frac M x L nbsp ex Ex B displaystyle e x E x B nbsp 其中B 2b 1 1 displaystyle B 2 b 1 1 nbsp 再继续看McEniry 2007里的进一步说明 y 1x displaystyle y frac 1 sqrt x nbsp 对等式的两边取二进制对数 log2 displaystyle textstyle log 2 nbsp 即函数f n 2n displaystyle textstyle f n 2 n nbsp 的反函数 有 log2 y 12log2 x displaystyle log 2 y frac 1 2 log 2 x nbsp log2 1 my ey 12log2 1 mx 12ex displaystyle log 2 1 m y e y frac 1 2 log 2 1 m x frac 1 2 e x nbsp 以如上方法 就能将浮点数x和y的相关指数消去 从而将乘方运算化为加法运算 而由于log2 x displaystyle textstyle log 2 x nbsp 与log2 x 1 2 displaystyle textstyle log 2 x 1 2 nbsp 线性相关 因此在x displaystyle textstyle x nbsp 与y0 displaystyle textstyle y 0 nbsp 即输入值与首次近似值 间就可以线性组合的方式建立方程 文 10 在此McEniry再度引入新数s displaystyle sigma nbsp 描述log2 1 x displaystyle textstyle log 2 1 x nbsp 与近似值R间的误差 註 3 由于0 x lt 1 displaystyle textstyle 0 leq x lt 1 nbsp 有log2 1 x x displaystyle textstyle log 2 1 x approx x nbsp 则在此可定义s displaystyle sigma nbsp 与x的关系为log2 1 x x s displaystyle textstyle log 2 1 x cong x sigma nbsp 这一定义就能提供二进制对数的首次精度值 此处0 s 13 displaystyle 0 leq sigma leq tfrac 1 3 nbsp 当R为0x5f3759df时 有s 0 0450461875791687011756 displaystyle textstyle sigma 0 0450461875791687011756 nbsp 文 13 由此将log2 1 x x s displaystyle textstyle log 2 1 x x sigma nbsp 代入上式 有 my s ey 12mx 12s 12ex displaystyle m y sigma e y frac 1 2 m x frac 1 2 sigma frac 1 2 e x nbsp 参照首段等式代入Mx displaystyle M x nbsp Ex displaystyle E x nbsp B displaystyle B nbsp 与L displaystyle L nbsp 有 My Ey B L 32sL 12Mx 12 Ex B L displaystyle M y E y B L frac 3 2 sigma L frac 1 2 M x frac 1 2 E x B L nbsp 移项整理得 EyL My 32 B s L 12 ExL Mx displaystyle E y L M y frac 3 2 B sigma L frac 1 2 E x L M x nbsp 如上所述 对于以浮点规格存储的正浮点数x 若将其作为长整型表示则示值为Ix ExL Mx displaystyle textstyle I x E x L M x nbsp 由此即可根据x的整数表示导出y 在此y 1x displaystyle textstyle y frac 1 sqrt x nbsp 亦即x的平方根倒数的首次近似值 的整数表示值 也即 Iy EyL My R 12 ExL Mx R 12Ix displaystyle I y E y L M y R frac 1 2 E x L M x R frac 1 2 I x nbsp 其中R 32 B s L displaystyle R frac 3 2 B sigma L nbsp 最后导出的等式Iy R 12Ix displaystyle textstyle I y R frac 1 2 I x nbsp 即与上节代码中i 0x5f3759df i gt gt 1 一行相契合 由此可见 在平方根倒数速算法中 对浮点数进行一次移位操作与整数减法 就可以可靠地输出一个浮点数的对应近似值 文 13 到此为止 McEniry只证明了 在常数R的辅助下 可近似求取浮点数的平方根倒数 但仍未能确定代码中的R值的选取方法 关于作一次移位与减法操作以使浮点数的指数被 2除的原理 Chris Lomont的论文中亦有个相对简单的解释 以10000 104 displaystyle textstyle 10000 10 4 nbsp 为例 将其指数除 2可得10000 1 2 10 2 1 100 displaystyle textstyle 10000 1 2 10 2 1 100 nbsp 而由于浮点表示的指数有进行过偏移处理 所以指数的真实值e应为e E 127 displaystyle textstyle e E 127 nbsp 因此可知除法操作的实际结果为 e 2 127 displaystyle textstyle e 2 127 nbsp 文 14 这时用R 在此即为 魔术数字 0x5f3759df 减之即可使指数的最低有效数位转入有效数字域 之后重新转换为浮点数时 就能得到相当精确的平方根倒数近似值 在这里对常数R的选取亦有所讲究 若能选取一个好的R值 便可减少对指数进行除法与对有效数字域进行移位时可能产生的错误 基于这一标准 0xbe即是最合适的R值 而0xbe右移一位即可得到0x5f 这恰是魔术数字R的第一个字节 文 15 精确度 编辑 nbsp 使用启发式平方根倒数速算法与使用C语言标准库libstdc的函数所计算出的平方根倒数的差值一览 注意这里使用的是双对数坐标系 如上所述 平方根倒数速算法所得的近似值惊人的精确 右图亦展示了以上述代码计算 以平方根倒数速算法计算后再进行一次牛顿法迭代 所得近似值的误差 当输入0 01时 以C语言标准库函数计算可得10 0 而InvSqrt 得值为9 9825822 其间误差为0 017479 相对误差则为0 175 且当输入更大的数值时 绝对误差不断下降 相对误差也一直控制在一定的范围之内 牛顿法提高精度 编辑 主条目 牛顿法 在进行了如上的整数操作之后 示例程序再度将被转为长整型的浮点数回转为浮点数 对应x float amp i 并对其进行一次浮点运算操作 对应x x 1 5f xhalf x x 这里的浮点运算操作就是对其进行一次牛顿法迭代 若以此例说明 y 1x displaystyle y frac 1 sqrt x nbsp 所求的是y的平方根倒数 以之构造以y为自变量的函数 有f y 1y2 x 0 displaystyle f y frac 1 y 2 x 0 nbsp 将其代入牛顿法的通用公式yn 1 yn f yn f yn displaystyle y n 1 y n frac f y n f y n nbsp 其中yn displaystyle y n nbsp 为首次近似值 有yn 1 yn 3 xyn2 2 displaystyle y n 1 frac y n 3 xy n 2 2 nbsp 其中f y 1y2 x displaystyle f y frac 1 y 2 x nbsp f y 2y3 displaystyle f y frac 2 y 3 nbsp 整理有yn 1 yn 3 xyn2 2 yn 1 5 xyn22 displaystyle y n 1 frac y n 3 xy n 2 2 y n 1 5 frac xy n 2 2 nbsp 对应的代码即为y y threehalfs xhalf y y 在以上一节的整数操作产生首次近似值后 程序会将首次近似值作为参数送入函数最后两句进行精化处理 代码中的两次迭代 以一次迭代的输出 对应公式中的yn 1 displaystyle y n 1 nbsp 作为二次迭代的输入 正是为了进一步提高结果的精度 文 16 但由于雷神之锤III引擎的图形计算中并不需要太高的精度 所以代码中只进行了一次迭代 二次迭代的代码则被注释 文 9 历史与考究 编辑 nbsp id Software的创始人约翰 卡马克 这段代码虽非他所作 但常被认为与他相关 雷神之锤III 的代码直到QuakeCon 2005才正式放出 但早在2002年 或2003年 时 平方根倒数速算法的代码就已经出现在Usenet与其他论坛上了 1 最初人们猜测是卡马克写下了这段代码 但他在回复询问他的邮件时否定了这个观点 并猜测可能是先前曾帮id Software优化雷神之锤的资深汇编程序员Terje Mathisen写下了这段代码 而在Mathisen的邮件裡 他表示 在1990年代初 他只曾作过类似的實作 确切来说这段代码亦非他所作 现在所知的最早實作是由Gary Tarolli在SGI Indigo中實作的 但他亦坦承他仅对常数R的取值做了一定的改进 实际上他也不是作者 在向以发明MATLAB而闻名的Cleve Moler查证后 Rys Sommefeldt则认为原始的算法是Ardent Computer 英语 Ardent Computer 公司的Greg Walsh所发明 但他也没有任何决定性的证据能证明这一点 5 不仅该算法的原作者不明 人们也仍无法確定当初选择这个 魔术数字 的方法 Chris Lomont曾做了个研究 他推算出了一个函数以討論此速算法的誤差 並找出了使誤差最小的最佳R值0x5f37642f 与代码中使用的0x5f3759df相当接近 但以之代入算法计算并进行一次牛顿迭代后 所得近似值之精度仍略低於代入0x5f3759df的结果 文 17 因此Lomont将目标改为尋找在进行1 2次牛顿迭代后能得到最大精度的R值 在暴力搜尋後得出最优R值为0x5f375a86 以此值代入算法并进行牛顿迭代 所得的结果都比代入原始值 0x5f3759df 更精确 文 17 于是他的论文最后提到 如果可能我想询问原作者 此速算法是以数学推导还是以反复试错的方式求得 文 18 在论文中 Lomont亦指出 64位的IEEE754浮点数 即双精度类型 所对应的魔术数字是0x5fe6ec85e7de30da 但后来的研究表明 代入0x5fe6eb50c7aa19f9的结果精确度更高 McEniry得出的结果则是0x5FE6EB50C7B537AA 精度介于两者之间 Matthew Robertson给出的精度更高的值是0x5FE6EB50C7B537A9 6 在Charles McEniry的论文中 他使用了一种类似Lomont但更复杂的方法来优化R值 他最开始使用穷举搜索 所得结果与Lomont相同 文 19 而后他尝试用带权二分法寻找最优值 所得结果恰是代码中所使用的魔术数字0x5f3759df 因此 McEniry认为 这一常数最初或许便是以 在可容忍误差范围内使用二分法 的方式求得 文 20 注释 编辑 由于现代计算机系统对长整型的定义有所差异 使用长整型会降低此段代码的可移植性 具体来说 由此段浮点转换为长整型的定义可知 如若这段代码正常运行 所在系统的长整型长度应为4字节 32位 否则重新转为浮点数时可能会变成负数 而由于C99标准的广泛应用 在现今多数64位计算机系统 除使用LLP64数据模型的Windows外 中 长整型的长度都是8字节 文 6 3 此处 浮点数 所指为标准化浮点数 也即有效数字部分必须满足0 mx lt 1 displaystyle textstyle 0 leq m x lt 1 nbsp 可参见David Goldberg What Every Computer Scientist Should Know About Floating Point Arithmetic ACM Computing Surveys March 1991 23 1 5 48 doi 10 1145 103162 103163 Lomont 2003确定R的方式则有所不同 他先将R分解为R1 displaystyle textstyle R 1 nbsp 与R2 displaystyle textstyle R 2 nbsp 其中R1 displaystyle textstyle R 1 nbsp 与R2 displaystyle textstyle R 2 nbsp 分别代表R的有效数字域和指数域 文 7 参考 编辑脚注 编辑 1 0 1 1 1 2 1 3 Sommefeldt Rys Origin of Quake3 s Fast InvSqrt Beyond3D 2006 11 29 2009 02 12 原始内容存档于2009 02 09 英语 quake3 1 32b code game q math c Quake III Arena id Software 2010 11 15 原始内容存档于2011 02 17 Meyers Randy The New C Integers in C99 Part 1 drdobbs com 2000 12 01 2010 09 04 原始内容存档于2012 05 27 Ruskin Elan Timing square root Some Assembly Required 2009 10 16 2010 09 13 原始内容存档于2010 08 25 Sommefeldt Rys Origin of Quake3 s Fast InvSqrt Part Two Beyond3D 2006 12 19 2008 04 19 原始内容存档于2008 05 13 Matthew Robertson A Brief History of InvSqrt PDF UNBSJ 2012 04 24 2023 11 23 原始内容存档 PDF 于2023 01 29 参考文献 编辑 Blinn 2003年 第130页 Middendorf 2007年 第155 164页 Eberly 2001年 第504页 Lomont 2003年 第1页 McEniry 2007年 第1页 Fog 2010年 第6页 7 0 7 1 Lomont 2003年 第3页 McEniry 2007年 第2 amp p 16页 9 0 9 1 Eberly 2002年 第2页 10 0 10 1 McEniry 2007年 第2页 Hennessey 1998年 第276页 Hennessey 1998年 第305页 13 0 13 1 13 2 McEniry 2007年 第3页 Hennessey 1998年 第278页 Lomont 2003年 第5页 McEniry 2007年 第6页 17 0 17 1 Lomont 2003年 第10页 Lomont 2003年 第10 11页 McEniry 2007年 第11 12页 McEniry 2007年 第16页 Blinn Jim Jim Blinn s Corner Notation notation notation Morgan Kaufmann 2003 ISBN 1 55860 860 5 Eberly David 3D Game Engine Design Morgan Kaufmann 2001 ISBN 978 1 55860 593 0 Eberly David Fast Inverse Square Root PDF Geometric Tools 2002 2009 03 22 原始内容 PDF 存档于2009 02 24 Fog Agner Calling conventions for different C compilers and operating systems Chapter 3 Data Representation PDF 2010 02 16 2010 08 30 原始内容 PDF 存档于2010 08 15 Hennessey John Patterson David A Computer Organization and Design 2nd San Francisco CA Morgan Kaufmann Publishers 1998 ISBN 978 1 55860 491 9 引文使用过时参数coauthors 帮助 Lomont Chris Fast Inverse Square Root PDF February 2003 2009 02 13 原始内容 PDF 存档于2009 02 06 McEniry Charles The Mathematics Behind the Fast Inverse Square Root Function Code PDF August 2007 2009 02 13 原始内容 PDF 存档于2015 05 11 Middendorf Lars Muhlbauer Felix Umlauf George Bodba Christophe Embedded Vertex Shader in FPGA Rettberg Achin 编 Embedded System Design Topics Techniques and Trends IFIP TC10 Working Conference International Embedded Systems Symposium IESS et al Irvine California Springer 2007 06 01 ISBN 978 0 387 72257 3 引文使用过时参数coauthors 帮助 延伸阅读 编辑Kushner David The wizardry of Id IEEE Spectrum Aug 2002 39 8 42 47 doi 10 1109 MSPEC 2002 1021943 顾森 Matrix67 旧闻一则 神秘的0x5f3759df 不可思议的Quake III源码 2007 11 24 2015 08 01 原始内容存档于2015 08 04 外部链接 编辑Origin of Quake3 s Fast InvSqrt 页面存档备份 存于互联网档案馆 Beyond3D com Quake III Arena source code 页面存档备份 存于互联网档案馆 id Software Margolin Tomer Magical Square Root Implementation In Quake III CodeMaestro The Coding Experience 2005 08 27 原始内容存档于2012 04 14 取自 https zh wikipedia org w index php title 平方根倒数速算法 amp oldid 80833370, 维基百科,wiki,书籍,书籍,图书馆,

文章

,阅读,下载,免费,免费下载,mp3,视频,mp4,3gp, jpg,jpeg,gif,png,图片,音乐,歌曲,电影,书籍,游戏,游戏。