开发者社区> 蒲烧> 正文
阿里云
为了无法计算的价值
打开APP
阿里云APP内打开

开方的几种写法

简介: 1.牛顿迭代 牛顿迭代其实不仅仅是用来做开平方的,而是快速逼近求得方程的逼近解。以下内容来自果壳文章《求牛顿开方法的算法及其原理,此算法能开任意次方吗?》 假设方程 在  附近有一个根,那么根据f(x)在x0附近的值和斜率,就能估计f(x)和x轴的交点用。
+关注继续查看

1.牛顿迭代

牛顿迭代其实不仅仅是用来做开平方的,而是快速逼近求得方程的逼近解。以下内容来自果壳文章《求牛顿开方法的算法及其原理,此算法能开任意次方吗?》

假设方程 c0da7cf18442d200c0963291fbe1e2f8ffdc4ff5在 efbda784ad565c1c5201fdc948a570d0426bc6e6 附近有一个根,那么根据f(x)在x0附近的值和斜率,就能估计f(x)和x轴的交点用。迭代式子:aefc36d26e2b81c7381b838b5cbc1f3c5dad173a 依次计算593f4cff5d4210d46e140db57bafc4f692493f76a8728ff397f08f1999170f64ff5838333f75538049510ae6a807501d0723acc06f628c6bbd2f68e2、……,那么序列将无限逼近方程的根。
用牛顿迭代法开平方】令:
1848f3b7365287b7962f2dcaed9ca482438d1c60 
所以f(x)的一次导是:
44df64cc9dc5767483c43d3c46fbd72cb78e7f97
牛顿迭代式:
768b411ac2ba3b8e23a2e85f2dfc565ecf2d13b9随便一个迭代的初始值,例如0c1d7f319728a07a57d000f2379b5215e4130147,代入上面的式子迭代。例如计算bfe16f27ebc966df6f10ba356a1547b6e7242dd7,即a=2。
0c1d7f319728a07a57d000f2379b5215e4130147
c1838533d54bb7d1bf09520443906c90c861cba2
ade97b4a977688e3486015b4ae430d6c902139b0
……
计算器上可给出477c57e54790c6a773dd634f07f2e8e5486d983e

 

根据果壳文章以上阐述,由此我们可以写出如下代码:


    /**
     * 牛顿迭代开平方
     */

public static double newtonInvertSqrt(double target) {
    double point = 1;
    double precision = 1e-2;
    while (Math.abs(target - point * point) > precision) {
        point = (point + target / point) / 2;    

    }

    return point;

}

求9的平方根,尝试 : 3次 , point : 3 , point的平方 : 9 , point平方与target的差值 : 0

到此为此,似乎问题已经得到了圆满的答案,但紧接着我就发现了更有趣的解法

2.魔法数开方倒数

暴雪是一家神奇的公司,它的代码总是能给人惊喜,例如 One-way Hash 被业内赞为最优秀的hash表改良算法(timer33算法应该才是最快的,One-way Hash 在Timer33的基础上改用两次Hash计算来比对Key是否一致,在数据结构没有根本变化的情况下,理论上来说应该更慢)。

下面是暴雪的开平方算法源码,它对上述的牛顿迭代做了一些改进:


#include <math.h>  
float InvSqrt(float x)  
{  
 float xhalf = 0.5f*x;  
 int i = *(int*)&x; // get bits for floating VALUE  
 i = 0x5f375a86- (i>>1); // gives initial guess y0  
 x = *(float*)&i; // convert bits BACK to float  
 x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy  
 return x;  
}  
  
int main()  
{  
  printf("%lf",1/InvSqrt(3));  
  
   return 0;  
}  

这个算法的牛逼之处在于,经典的牛顿迭代里猜测值是随便取的一个值,而这段代码里却用了一个魔法数0x5f3759df来作为猜测值。这个算法只用了一次位运算,根本没有多次迭代就得到了结果,比直接用sqrt(n)快了大约四倍。

但是,我要说但是了,这个算法并不是暴雪原创的。它更早的时候出现在 雷神之锤3 这个游戏的3D引擎代码里,作者是约翰-卡马克(John Carmack),江湖人称卡神。

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

#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}

 注意看代码
i = 0x5f3759df - ( i >> 1 ); // what the fuck?

很显然,卡神也不知道这个魔法数是哪来的。根据gamedev一群大神的挖掘,发现早在70年代NASA的代码里就有了这个魔法数。

这个魔法数 0x5f3759df 。它是从哪来的,又起到了什么作用?
普渡大学的数学家Chris Lomont在gamedev上看到这个讨论以后觉得有趣,决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。为此写下一篇论文,Fast Inverse Square Root 并通过暴力穷举的方式,得出一个更好的开方猜测值:0x5f375a86

为什么这个值是0x5f3759df,Chris Lomont的论文Fast Inverse Square Root给出了答案(点我看WIKI)知乎上有同学给出了汉化版答案(点我看知乎原文):

单精度的浮点数结构是这样:

那么现在,每个正浮点数 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 );  
的秘密所在。

程序员多学点数学,总是会有用的,与君共勉
 
附录:java版算法实现

public static double fastInvertSqrt(double number) {
    long i;
    double x2, y;
    final double threehalfs = 1.5F;
    x2 = number * 0.5F;
    y = number;

    // evil floating point bit level hacking
    i = Double.doubleToLongBits(y);

    // gives initial guess y0
    i = 0x5fe6eb50c7b537a9L - (i >> 1);

    y = Double.longBitsToDouble(i);

    // 1st iteration
    y = y * (threehalfs - (x2 * y * y));

    // 2nd iteration
    y = y * (threehalfs - (x2 * y * y));

    return 1 / y;
}


引用文献:

[1].果壳《求牛顿开方法的算法及其原理,此算法能开任意次方吗? 》http://www.guokr.com/question/461510/

[2].维基百科《Fast inverse square root》https://en.wikipedia.org/wiki/Fast_inverse_square_root

[3].知乎《0x5f3759df这个快速开方中的常数的数学依据是什么?》https://www.zhihu.com/question/20690553

版权声明:本文内容由阿里云实名注册用户自发贡献,版权归原作者所有,阿里云开发者社区不拥有其著作权,亦不承担相应法律责任。具体规则请查看《阿里云开发者社区用户服务协议》和《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容,填写侵权投诉表单进行举报,一经查实,本社区将立刻删除涉嫌侵权内容。

相关文章
lodash两数相除
lodash两数相除
0 0
lodash两数相减
lodash两数相减
0 0
考点:函数参数传参、求和、奇数、偶数、输入输出、range步长灵活使用【Python习题04】
考点:函数参数传参、求和、奇数、偶数、输入输出、range步长灵活使用【Python习题04】
0 0
复习C部分:1.看代码求值题 2.写三个整数代码从大到小输出 3.打印1~100中所有3的倍数 4.给定两个数,求最大公约数(递减法,辗转相除法)
复习C部分:1.看代码求值题 2.写三个整数代码从大到小输出 3.打印1~100中所有3的倍数 4.给定两个数,求最大公约数(递减法,辗转相除法)
0 0
【初识C语言】简单计算n的阶乘以及多个数字的阶乘求和
【初识C语言】简单计算n的阶乘以及多个数字的阶乘求和
0 0
力扣题 两数相除:画图解析 采用递归计算除法(不使用乘法、除法和 mod 运算符)
这是力扣上的一道题目,难度为中等,两数相除:给定两个整数,被除数 dividend 和除数 divisor。将两数相除,要求不使用乘法、除法和 mod 运算符。
0 0
二分法的两种写法
二分法的两种写法
0 0
F#表达式化简
如果使用过Matlab或者Maple等软件,应该知道这类数学软件的符号计算引擎非常强大,可以进行数学公式的推导,比如可以对数学公式进行化简。当然,实现一个功能完备的化简引擎还是不容易的。这里用F# 实现一个简单的化简函数: ( a + x )+( a - x ) => 2 * a
0 0
最大公约数的几种写法
    特点及意义 最大公约数指某几个整数共有因子中最大的一个。 GCD即Greatest Common Divisor. 例如,12和30的公约数有:1、2、3、6,其中6就是12和30的最大公约数。
785 0
+关注
蒲烧
nerd and proud
文章
问答
文章排行榜
最热
最新
相关电子书
更多
低代码开发师(初级)实战教程
立即下载
阿里巴巴DevOps 最佳实践手册
立即下载
冬季实战营第三期:MySQL数据库进阶实战
立即下载