【算法】非典型牛顿迭代法

何为牛顿迭代法

​ 牛顿迭代法的全名叫做牛顿-拉夫森方法(Newton-Raphson method),这是一种在实数域和复数域上近似求解方程的方法。方法使用函数f(x)泰勒级数的前面几项来寻找方程f(y)=0的根。

​ 对于一般的形如f(x)=0的方程,首先任意估算一个解$x1线x=x_1f(x)线x(x_2,0)x_2f(x)=0x{n+1}$,就能无限接近于方程的精确解。

newton

​ 值得注意的是,该方法在函数不连续或者函数有多个零点的时候会失效。


代码实现

​ 牛顿迭代法多用于以最快速度求一个数字m方根,仅需几次的迭代就能够得到非常精确的结果,下面举例来介绍牛顿迭代法求平方根的数学公式和实现

​ 不妨令xm=n,则方程为f(x)=xmn.

​ 有

f(x)=mxm1
xn+1=xnf(xn)f(xn)=(1m)xn+nxnmxmn

​ 假设我们输入一个数字n,我们现在需要求它的平方根x,则

f(x)=x2n
f(x)=2x
xn+1=xnf(xn)f(xn)=12(xn+nxn)

​ 则我们的迭代方程就是xn+1=xn2+n2xn

​ 现在要确定的是迭代终止条件,假设我们要求的精度是0.00001,则迭代终止的条件为|x2n|<0.00001

​ 由于x=0的时候,n的值就为0,所以我们将迭代的起点设置为x=1,所以我们的程序如下:

1
2
3
4
5
6
7
8
9
10
11
12
#Python3
def mySqrt(self, n: int) -> float:
if n == 0 or n == 1:
return n
#先处理特殊值
x0 = 1.0
ran = abs(x0 * x0 - n) #迭代停止信号
while ran > 0.00001:
x1 = x0 / 2 + n / (x0 * 2)
x0 = x1
ran = abs(x1 * x1 - n)
return x1

​ 当我们输入的值为2147395599时,输出46339.99998,输出耗时为32ms。而Python自带的n**0.5耗时60ms。


非典型算法

​ 这里不得不提到神秘数字0x5f3759df,这个数字首次出现在游戏《雷神之锤III》的源代码中。

​ 在一个名为q_math.c的文件里出现了如下代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//C
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;
}

​ 里面出现了一句让人抓毛的代码:

1
i = 0x5f3759df - ( i >> 1 ); // what the fuck?

​ 这句代码和下面的两行一起

1
2
3
i = 0x5f3759df - ( i >> 1 ); // what the fuck? 
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration

​ 就相当于求平方根的操作。

​ 但其中的原理是什么,答案恐怕连写这个代码的人自己也不知道。(what the f*ck?)


参考

【1】维基百科:牛顿法

【2】力扣题目:求平方根

-------------FIN-------------

本文标题:【算法】非典型牛顿迭代法

文章作者:吃草莓糖葫芦

发布时间:2020年01月31日 - 10:01

最后更新:2020年02月03日 - 15:02

原始链接:https://tsuinterukonsigure.github.io/2020/01/31/%E3%80%90%E7%AE%97%E6%B3%95%E3%80%91%E9%9D%9E%E5%85%B8%E5%9E%8B%E7%89%9B%E9%A1%BF%E8%BF%AD%E4%BB%A3%E6%B3%95/

许可协议: 署名-非商业性使用-禁止演绎 4.0 国际 转载请保留原文链接及作者。

如果文章对你有用,不妨捐助一下作者小哥哥叭~