1.牛顿迭代
牛顿迭代其实不仅仅是用来做开平方的,而是快速逼近求得方程的逼近解。以下内容来自果壳文章《求牛顿开方法的算法及其原理,此算法能开任意次方吗?》
假设方程 在 附近有一个根,那么根据f(x)在x0附近的值和斜率,就能估计f(x)和x轴的交点用。迭代式子: 依次计算 、 、 、……,那么序列将无限逼近方程的根。
用牛顿迭代法开平方】令:
所以f(x)的一次导是:
牛顿迭代式:
随便一个迭代的初始值,例如 ,代入上面的式子迭代。例如计算 ,即a=2。
……
计算器上可给出
根据果壳文章以上阐述,由此我们可以写出如下代码:
/**
* 牛顿迭代开平方
*/
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 可以用尾数和指数的形式写成 ,其中 m 是尾数部分,取值范围是 ;e 是指数部分,一个整数。每个浮点数所对应的「整数形式」则可以用表示,其中 L 是指数部分需要的位移次数(用 2 的幂表示),E 和 M 是指数部分和小数部分的整数版本。两者之间的关系是
对单精度浮点而言,。
考虑对数 ,由于,,取近似,可以算出整体「偏差」最小的,此时两者基本相当。因此我们可以说
……………… (1)
那么,对于 y 的整数形式 Y 而言,展开并带代入 (1) 有:
即
………………(2)
那么对于 来说,,代入 (2) 得:
解得
。
这个就是代码
i = 0x5f3759df - ( i >> 1 );
程序员多学点数学,总是会有用的,与君共勉
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