0x5f3759df这个快速开方中的常数的数学依据和原理

2024-08-23 11:32

本文主要是介绍0x5f3759df这个快速开方中的常数的数学依据和原理,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

 

Quake-III Arena (雷神之锤3)是90年代的经典游戏之一。

 

该系列的游戏不但画面和内容不错,而且即使计算机配置低,也能极其流畅地运行。这要归功于它3D引擎的开发者约翰-卡马克(John Carmack)。事实上早在90年代初DOS时代,只要能在PC上搞个小动画都能让人惊叹一番的时候,John Carmack就推出了石破天惊的Castle Wolfstein, 然后再接再励,doom, doomII, Quake...每次都把3-D技术推到极致。他的3D引擎代码资极度高效,几乎是在压榨PC机的每条运算指令。当初MS的Direct3D也得听取他的意见,修改了不少API。


最近,QUAKE的开发商ID SOFTWARE遵守GPL协议,公开了QUAKE-III的原代码,让世人有幸目睹Carmack传奇的3D引擎的原码。


我们知道,越底层的函数,调用越频繁。3D引擎归根到底还是数学运算。

那么找到最底层的数学运算函数(game/code/q_math.c),必然是精心编写的。里面有很多有趣的函数,很多都令人惊奇,估计我们几年时间都学不完。


在/code/game/q_math.c里发现了这样一段代码。

它的作用是将一个数开平方并取倒,经测试这段代码比(float)(1.0/sqrt(x))快4倍:

[cpp]  view plain copy
在CODE上查看代码片 派生到我的代码片
  1. float Q_rsqrt( float number )  
  2. {  
  3.     long i;  
  4.     float x2, y;  
  5.     const float threehalfs = 1.5F;  
  6.   
  7.     x2 = number * 0.5F;  
  8.     y  = number;  
  9.     i  = * ( long * ) &y;                       // evil floating point bit level hacking  
  10.     i  = 0x5f3759df - ( i >> 1 );               // what the fuck?  
  11.     y  = * ( float * ) &i;  
  12.     y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration  
  13. //  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed  
  14.   
  15. #ifndef Q3_VM  
  16. #ifdef __linux__  
  17.     assert( !isnan(y) ); // bk010122 - FPE?  
  18. #endif  
  19. #endif  
  20.     return y;  
  21. }  


函数返回1/sqrt(x),这个函数在图像处理中比sqrt(x)更有用。

注意到这个函数只用了一次叠代!(其实就是根本没用叠代,直接运算)。编译,实验,这个函数不仅工作的很好,而且比标准的sqrt()函数快4倍!要知道,编译器自带的函数,可是经过严格仔细的汇编优化的啊!


这个简洁的函数,最核心,也是最让人费解的,就是标注了“what the fuck?”的一句:i = 0x5f3759df - ( i >> 1 );


再加上y = y * ( threehalfs - ( x2 * y * y ) );

两句话就完成了开方运算!而且注意到,核心那句是定点移位运算,速度极快!特别在很多没有乘法指令的RISC结构CPU上,这样做是极其高效的。


算法的原理其实不复杂,就是牛顿迭代法,用x-f(x)/f'(x)来不断的逼近f(x)=a的根。

简单来说比如求平方根,f(x)=x^2=a ,f'(x)= 2*x,f(x)/f'(x)=x/2,把f(x)代入x-f(x)/f'(x)后有(x+a/x)/2,现在
我们选a=5,选一个猜测值比如2,那么我们可以这么算

 

5/2 = 2.5;

(2.5+2)/2 = 2.25;

5/2.25 = xxx;

(2.25+xxx)/2 = xxxx

...


这样反复迭代下去,结果必定收敛于sqrt(5),没错,一般的求平方根都是这么算的,但是卡马克(quake3作者)真正牛B的地
方是他选择了一个神秘的常数0x5f3759df 来计算那个猜测值,就是我们加注释的那一行,那一行算出的值非常接近1/sqrt(n),这样我们只需要2次牛顿迭代就可以达到我们所需要的精度.好吧 如果这个还不算NB,接着看:



普渡大学的数学家Chris Lomont看了以后觉得有趣,决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。Lomont也是个
牛人,在精心研究之后从理论上也推导出一个最佳猜测值,和卡马克的数字非常接近, 0x5f37642f。卡马克真牛,他是外星人吗?


传奇并没有在这里结束。Lomont计算出结果以后非常满意,于是拿自己计算出的起始值和卡马克的神秘数字做比赛,看看谁的
数字能够更快更精确的求得平方根。结果是卡马克赢了... 谁也不知道卡马克是怎么找到这个数字的。


最后Lomont怒了,采用暴力方法一个数字一个数字试过来,终于找到一个比卡马克数字要好上那么一丁点的数字,虽然实际上
这两个数字所产生的结果非常近似,这个暴力得出的数字是0x5f375a86。


Lomont为此写下一篇论文,"Fast Inverse Square Root"。

最后,给出最精简的1/sqrt()函数:

[cpp]  view plain copy
在CODE上查看代码片 派生到我的代码片
  1. float InvSqrt(float x)  
  2. {  
  3.     float xhalf = 0.5f * x;  
  4.     int i = *(int *)&x;  
  5.     i = 0x5f3759df - (i>>1);  
  6.     x = *(float *)&i;  
  7.     x = x * (1.5f - xhalf * x * x);  
  8.     return x;  
  9. }  

同时:


单精度的浮点数结构是这样:
那么现在,每个正浮点数 y 可以用尾数和指数的形式写成  2^e(1+m),其中 m 是尾数部分,取值范围是  [0, 1);e 是指数部分,一个整数。每个浮点数所对应的「整数形式」则可以用 EL+M表示,其中 L 是指数部分需要的位移次数(用 2 的幂表示),E 和 M 是指数部分和小数部分的整数版本。两者之间的关系是
\left\{ \begin{array}{ll}E  =  e + B \\M  =  mL\end{array} \right.
对单精度浮点而言, L=2^{23},B=127

考虑对数  \log_2 y=e+\log_2(1+m),由于 m\in\left[0, 1\right)\log_2(1+m)\in\left[0,1\right),取近似 \log_2(1+m)\approx m+\delta,可以算出整体「偏差」最小的 \delta=0.0430357 ,此时两者基本相当。因此我们可以说 

\log_2 y\approx e+m+\delta……………… (1)

那么,对于 y 的整数形式 Y 而言,展开并带代入 (1) 有:
Y & = & EL+M\\&= &L(e+B+m) \\& \approx & L(\log_2 y + B - \delta)


\log_2 y\approx {Y \over L}-(B-\delta)………………(2)

那么对于  y'={1 \over \sqrt{y}}来说, \log_2{y'}\approx -{1\over 2}\log_2 y,代入 (2) 得:
{Y' \over L}-(B-\delta)={1 \over 2}\left[(B-\delta)-{Y\over L}\right]
解得
Y'={3 \over 2}L(B-\delta)-{1 \over 2}Y
这个就是代码
i  = 0x5F3759DF - ( i >> 1 );
的秘密所在。0x5F3759DF 的具体取值是根据 δ 变的,至于卡马克为啥用 0x5F3759DF 而不是其他相近的值,估计是做实验测的吧……

同时:

假设原浮点数为0EM,

E为8位指数部分,M为23位尾数部分。

位操作的主要目的是为了让新的浮点数的e = -(e0)/2

我们知道E = e0+127

而我们想让位操作后的E'=-e0/2+127

右移一位可以得到:

E_s = E/2 = e0/2 + 127/2

如果我们用190去减E_s, 就可以得到

E_1 = 190-E_s=190-e0/2-63=-e0/2+127

而190正是这个魔法数的第31-23位。


因此这个魔法数帮助我们得到了正确的新浮点数的指数部分。小数部分则比较复杂,因为在原指数部分为奇数的情况下右移一位会无端端给尾数增加1/2,且相减后尾数可能会借用指数部分导致指数变小,因此通过调整尾数来修正。如需详细了解可参考这篇paper: lomont.org/Math/Papers/


需要指出现在的Intel CPU已经提供快速rsqrt指令,所以无需自己去实现这个算法。


这篇关于0x5f3759df这个快速开方中的常数的数学依据和原理的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1099252

相关文章

深入探索协同过滤:从原理到推荐模块案例

文章目录 前言一、协同过滤1. 基于用户的协同过滤(UserCF)2. 基于物品的协同过滤(ItemCF)3. 相似度计算方法 二、相似度计算方法1. 欧氏距离2. 皮尔逊相关系数3. 杰卡德相似系数4. 余弦相似度 三、推荐模块案例1.基于文章的协同过滤推荐功能2.基于用户的协同过滤推荐功能 前言     在信息过载的时代,推荐系统成为连接用户与内容的桥梁。本文聚焦于

电脑桌面文件删除了怎么找回来?别急,快速恢复攻略在此

在日常使用电脑的过程中,我们经常会遇到这样的情况:一不小心,桌面上的某个重要文件被删除了。这时,大多数人可能会感到惊慌失措,不知所措。 其实,不必过于担心,因为有很多方法可以帮助我们找回被删除的桌面文件。下面,就让我们一起来了解一下这些恢复桌面文件的方法吧。 一、使用撤销操作 如果我们刚刚删除了桌面上的文件,并且还没有进行其他操作,那么可以尝试使用撤销操作来恢复文件。在键盘上同时按下“C

hdu4407(容斥原理)

题意:给一串数字1,2,......n,两个操作:1、修改第k个数字,2、查询区间[l,r]中与n互质的数之和。 解题思路:咱一看,像线段树,但是如果用线段树做,那么每个区间一定要记录所有的素因子,这样会超内存。然后我就做不来了。后来看了题解,原来是用容斥原理来做的。还记得这道题目吗?求区间[1,r]中与p互质的数的个数,如果不会的话就先去做那题吧。现在这题是求区间[l,r]中与n互质的数的和

uva 10014 Simple calculations(数学推导)

直接按照题意来推导最后的结果就行了。 开始的时候只做到了第一个推导,第二次没有继续下去。 代码: #include<stdio.h>int main(){int T, n, i;double a, aa, sum, temp, ans;scanf("%d", &T);while(T--){scanf("%d", &n);scanf("%lf", &first);scanf

uva 10025 The ? 1 ? 2 ? ... ? n = k problem(数学)

题意是    ?  1  ?  2  ?  ...  ?  n = k 式子中给k,? 处可以填 + 也可以填 - ,问最小满足条件的n。 e.g k = 12  - 1 + 2 + 3 + 4 + 5 + 6 - 7 = 12 with n = 7。 先给证明,令 S(n) = 1 + 2 + 3 + 4 + 5 + .... + n 暴搜n,搜出当 S(n) >=

uva 11044 Searching for Nessy(小学数学)

题意是给出一个n*m的格子,求出里面有多少个不重合的九宫格。 (rows / 3) * (columns / 3) K.o 代码: #include <stdio.h>int main(){int ncase;scanf("%d", &ncase);while (ncase--){int rows, columns;scanf("%d%d", &rows, &col

【生成模型系列(初级)】嵌入(Embedding)方程——自然语言处理的数学灵魂【通俗理解】

【通俗理解】嵌入(Embedding)方程——自然语言处理的数学灵魂 关键词提炼 #嵌入方程 #自然语言处理 #词向量 #机器学习 #神经网络 #向量空间模型 #Siri #Google翻译 #AlexNet 第一节:嵌入方程的类比与核心概念【尽可能通俗】 嵌入方程可以被看作是自然语言处理中的“翻译机”,它将文本中的单词或短语转换成计算机能够理解的数学形式,即向量。 正如翻译机将一种语言

hdu 4565 推倒公式+矩阵快速幂

题意 求下式的值: Sn=⌈ (a+b√)n⌉%m S_n = \lceil\ (a + \sqrt{b}) ^ n \rceil\% m 其中: 0<a,m<215 0< a, m < 2^{15} 0<b,n<231 0 < b, n < 2^{31} (a−1)2<b<a2 (a-1)^2< b < a^2 解析 令: An=(a+b√)n A_n = (a +

v0.dev快速开发

探索v0.dev:次世代开发者之利器 今之技艺日新月异,开发者之工具亦随之进步不辍。v0.dev者,新兴之开发者利器也,迅速引起众多开发者之瞩目。本文将引汝探究v0.dev之基本功能与优势,助汝速速上手,提升开发之效率。 何谓v0.dev? v0.dev者,现代化之开发者工具也,旨在简化并加速软件开发之过程。其集多种功能于一体,助开发者高效编写、测试及部署代码。无论汝为前端开发者、后端开发者

hdu4407容斥原理

题意: 有一个元素为 1~n 的数列{An},有2种操作(1000次): 1、求某段区间 [a,b] 中与 p 互质的数的和。 2、将数列中某个位置元素的值改变。 import java.io.BufferedInputStream;import java.io.BufferedReader;import java.io.IOException;import java.io.Inpu