计算机与数学 —— 雷神之锤3源码中的快速逆平方根算法
2016-06-01 18:39
507 查看
这篇博客介绍了在雷神之锤3源代码中快速求逆平方根的算法。
函数接受一个浮点数作为输入,输出的结果是平方根的倒数。
在这32位中,最高位为符号位,后面的8位为整数ex,代表浮点数的指数,而最后的23位表示的是小数部分mx,小数部分第一位表示2−12^{-1},第二位表示2−22^{-2}……
所以如果x是一个正浮点数,则有:
x=2ex(1+mx) x = 2^{e_{x}} (1+m_{x})
另外,如果需要把一个浮点数转化为整数形式,则需要做如下的运算:
Ix=EXL+Mx=L(ex+B+mx)Ix = E_XL + M_x=L(e_x+B+m_x)
在上面的公式中,L指的是指数部分需要的唯一次数,这里是2232^{23},B是127,而M是小数部分对应的整数版本。
假设对于如下的函数,我们想要求这个函数的根:
该如何求得这个根呢?首先,我们先猜一个解x0x_{0},并且认为它是函数的解。
但是由于它其实并不是函数的解,那么现在,我们需要将这个解进行迭代,从而让其逼近真正的解。
在此处,我们可以在(x0,y0)(x_{0}, y_{0})处作其切线,求得该直线的方程:
y−f(x0)=f′(x0)(x−x0)y - f(x_0) = f'(x_0)(x - x_0)
并且求直线的根,此时会发现已经对于真正的解逼近了一步。
推广到n,继续迭代,就可以足够逼近真正的解了:
xi+1=xi−f(xi)f′(xi) x_{i+1} = x_{i} - \dfrac{f(x_i)}{f'(x_i)}
此时发现f(xi)f(x_i)与f′(xi)f'(x_i)可以被一个统一的函数g(xi)g(x_i)来表示:
g(x)=f(x)f′(x)g(x) = \dfrac{f(x)}{f'(x)}
令ε为当前的解与真正解r的距离:
ϵi=xi−r\epsilon_i = x_i - r
综合上面三个方程,可得:
ϵi+1=ϵi−g(xi)\epsilon_{i+1} = \epsilon_i - g(x_i)
因此只要ε值小于某个特定的值,我们可以认为此时的x和方程的解已经很接近了。
y=1x√y = \dfrac{1}{\sqrt{x}}
转化为关于y的方程,有:
f(y)=1y2−x=0f(y) = \dfrac{1}{y^2} - x = 0
转化为牛顿发使用的方程,有:
yn+1=yn(3−xy2n)2y_{n+1} = \dfrac{y_n(3 - xy_n^2)}{2}
此时,对原本的方程等号两边同时取2的对数,就有:
log2y=−12log2(1+mx)\log_{2}{y} = -\dfrac{1}{2}\log_{2}{(1+m_x)}
因为mx≥0m_x\ge0并且mx<1m_x<1,那么在这个区间内,可以取近似为:
log2(1+mx)≈mx+σ\log_2(1+m_x)\approx m_x + \sigma
根据方差的计算,当σ=0.0430357\sigma=0.0430357时,整体的偏差是最小的,此时上面的等号两边应该相当。
因此,把上面的完全整合起来,最终的IxI_x可以写成:
Llog2x+L(B−σ)L\log_2{x}+L(B-\sigma)
则Iy≈32L(B−σ)−12IxI_y\approx \dfrac{3}{2}L(B-\sigma)-\dfrac{1}{2}I_x
最终,写成代码就是:
32L(B−σ)=0x5f3759df\dfrac{3}{2}L(B-\sigma) = 0x5f3759df
<全文完>
源码
雷神之锤3中的逆平方根算法如下: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 ; }
函数接受一个浮点数作为输入,输出的结果是平方根的倒数。
单精度浮点数
在计算机中,单精度浮点数使用32位来储存表示:在这32位中,最高位为符号位,后面的8位为整数ex,代表浮点数的指数,而最后的23位表示的是小数部分mx,小数部分第一位表示2−12^{-1},第二位表示2−22^{-2}……
所以如果x是一个正浮点数,则有:
x=2ex(1+mx) x = 2^{e_{x}} (1+m_{x})
另外,如果需要把一个浮点数转化为整数形式,则需要做如下的运算:
Ix=EXL+Mx=L(ex+B+mx)Ix = E_XL + M_x=L(e_x+B+m_x)
在上面的公式中,L指的是指数部分需要的唯一次数,这里是2232^{23},B是127,而M是小数部分对应的整数版本。
牛顿法
牛顿迭代法(简称牛顿法),是用于计算机求解任意连续函数的根值的一种方法。假设对于如下的函数,我们想要求这个函数的根:
该如何求得这个根呢?首先,我们先猜一个解x0x_{0},并且认为它是函数的解。
但是由于它其实并不是函数的解,那么现在,我们需要将这个解进行迭代,从而让其逼近真正的解。
在此处,我们可以在(x0,y0)(x_{0}, y_{0})处作其切线,求得该直线的方程:
y−f(x0)=f′(x0)(x−x0)y - f(x_0) = f'(x_0)(x - x_0)
并且求直线的根,此时会发现已经对于真正的解逼近了一步。
推广到n,继续迭代,就可以足够逼近真正的解了:
xi+1=xi−f(xi)f′(xi) x_{i+1} = x_{i} - \dfrac{f(x_i)}{f'(x_i)}
此时发现f(xi)f(x_i)与f′(xi)f'(x_i)可以被一个统一的函数g(xi)g(x_i)来表示:
g(x)=f(x)f′(x)g(x) = \dfrac{f(x)}{f'(x)}
令ε为当前的解与真正解r的距离:
ϵi=xi−r\epsilon_i = x_i - r
综合上面三个方程,可得:
ϵi+1=ϵi−g(xi)\epsilon_{i+1} = \epsilon_i - g(x_i)
因此只要ε值小于某个特定的值,我们可以认为此时的x和方程的解已经很接近了。
算法分析
如果需要求得一个浮点数的平方根倒数,方程如下:y=1x√y = \dfrac{1}{\sqrt{x}}
转化为关于y的方程,有:
f(y)=1y2−x=0f(y) = \dfrac{1}{y^2} - x = 0
转化为牛顿发使用的方程,有:
yn+1=yn(3−xy2n)2y_{n+1} = \dfrac{y_n(3 - xy_n^2)}{2}
此时,对原本的方程等号两边同时取2的对数,就有:
log2y=−12log2(1+mx)\log_{2}{y} = -\dfrac{1}{2}\log_{2}{(1+m_x)}
因为mx≥0m_x\ge0并且mx<1m_x<1,那么在这个区间内,可以取近似为:
log2(1+mx)≈mx+σ\log_2(1+m_x)\approx m_x + \sigma
根据方差的计算,当σ=0.0430357\sigma=0.0430357时,整体的偏差是最小的,此时上面的等号两边应该相当。
因此,把上面的完全整合起来,最终的IxI_x可以写成:
Llog2x+L(B−σ)L\log_2{x}+L(B-\sigma)
则Iy≈32L(B−σ)−12IxI_y\approx \dfrac{3}{2}L(B-\sigma)-\dfrac{1}{2}I_x
最终,写成代码就是:
i = 0x5f3759df - ( i >> 1 );
32L(B−σ)=0x5f3759df\dfrac{3}{2}L(B-\sigma) = 0x5f3759df
<全文完>
相关文章推荐
- cogs_14_搭配飞行员_(二分图匹配+最大流,网络流24题#01)
- Broadcast Receiver监听网络状态
- 83款 网络爬虫开源软件
- 目前网络上开源的网络爬虫以及一些简介和比较
- 数据结构 - 用递归算法解决实际问题
- Longest Palindromic Substring --leetcode 数据结构第五题
- Spring Boot: HttpMediaTypeNotAcceptableException: Could not find acceptable representation原因及解决方法
- 数据结构 - 串的基本运算实现
- qt ,使用tcp/ip协议网络传输数据时,字节序转换方法
- 数据结构 - 栈和队列的基本运算实现
- 数据结构—线索化二叉树(中序)
- 聊聊HTTPS和SSL/TLS协议
- ios学习--网络流量统计
- HTTP的请求方法OPTIONS
- 后台发送http请求 类
- Java正则、Mysql、网络编程、多线程总结
- B+树 习题解
- 通过HttpURLConnection模拟post表单提交
- HTTP状态码
- ES6学习笔记(六)--set,map数据结构和for...of遍历