牛顿迭代算法 sqrt函数

因为吹水的能力不佳，所以要先打个草稿，今天的吹水过程大概是：
1、牛顿迭代法的演绎过程
2、牛顿迭代法求n次方根
3、牛顿迭代法求n次方根改进版
4、牛逼哄哄的invsqrt求平方根倒数

1、牛顿迭代法的演绎过程

乍一听，好像很高大上，其实理解上没什么难度。
牛顿迭代法就是在不断迭代的过程中逼近曲线的根，可以看做，它就是用来求曲线的根的。
所以它是怎么个迭代的过程，先上个动图：

牛顿迭代法的过程

我用吹水的姿势描述一遍动图的过程，先是有一个曲线，然后要找到曲线和x轴的交点，这个时候，脑袋一拍，选x1作为第一点，在这个点上垂直于x轴的线与曲线的交点，在这个交点做与曲线的切线，而切线与x轴的交点就是下一个迭代的点，不断重复这个过程，直到找到曲线和x轴交点这个目标值。值得注意的是，这个方法是不断逼近，所以得到的值可能会与目标值存在误差，而误差小于一定范围就可以认为是目标值了。一个更加图文并茂更加通俗易懂的解释可以点击 这个这个这个这个。

2、牛顿迭代法求n次方根

那为什么说牛顿迭代法能用来求n次根呢，其实在1中的目标值，即曲线与x轴相交点的x值，假设你的方程为 f(x) = x² - 64，当f(x) = 0时，x就是64的2次方根，以此类推，把2换成你想要求的n次方，就可以求出64的n次方根了。

明白了演绎过程，那我们就来探索它的代数实现。在上代码前，先上个图

 

牛顿迭代法某一点的求值图解

每次迭代要求的点就是右绿色箭头指向的这个点，在图中可以看到黑色箭头标识的线段长度为-f(x)，f'(x)是对f(x)的求导结果，也就是切线的斜率，所以绿色箭头线段的长度为-f(x)/f'(x)，而左绿色箭头的点为x，所有右绿色箭头指向的值为x加上绿色线段的长度，即x-(f(x)/f'(x))。

上面的代数实现完成后，就可以上代码了。嘻嘻嘻，最近在看Python版的SICP，所以贴上Python的代码，这篇文章也是因为看了这本书的1.6章节写的。def就是定义一个方法的意思。

 

def newton_nth_root_of_a(n, a, x = 1, e = 1e-13):
    def f(x):
        return x ** n - a
    def df(x):
        return n * x ** (n - 1)
    
    gap = abs(f(x))
    while gap > e:
        x = x - f(x)/df(x)
        gap = abs(f(x))
    return x
    
print(newton_nth_root_of_a(1, 64))  #64.0
print(newton_nth_root_of_a(2, 64))  #8.0
print(newton_nth_root_of_a(3, 64))  #4.0
print(newton_nth_root_of_a(4, 64))  #2.82842712474619
print(newton_nth_root_of_a(5, 64))  #2.29739670999407
print(newton_nth_root_of_a(6, 64))  #2.0

这里是在线可运行版本 ，不多不少，刚好十行，newton_nth_root_of_a方法定义中，n代表幂，a代表对a求根，x代表起始值，e代表误差范围，x和e已经设置了缺省值。f(x)是曲线，df(x)是f(x)的求导结果，在while循环中，x - f(x)/df(x)就是我们上面代数实现中的结果，直到f(x)的值小于误差，x就是最后的结果。

3、牛顿迭代法求n次方根改进版

明白了1、2其实整个过程就是这样了，那为啥有3，因为要硬着头皮接着吹水😂。
2中的10行代码很简洁，但是假如我要不断的求某些值的2次根，就要不断调用newton_nth_root_of_a(2, a)。那么我们来改进下代码，保持上面的代码不变，增加下面的代码：

 

def curry_nth_root(n):
    def curry_nth_root_of_a(a):
        return newton_nth_root_of_a(n, a)    
    return curry_nth_root_of_a

_2th_root = curry_nth_root(2)

print(_2th_root(64)) #8.0
print(_2th_root(49)) #7.000000000000001
print(_2th_root(36)) #6.0
print(_2th_root(25)) #5.0
print(_2th_root(16)) #4.0
print(_2th_root(9))  #3.0

这里是在线可运行版本 ，哇咔咔，在curry化之后，就可以方便的求n次方根了，不用每次都输入n。

curry化其实也很简单，那么我们就来个极端点的抽象代码：

 

def improve_guess(update, close, guess=1):
    '''返回猜想值，update为更新方法，close为逼近方法，guess为初始猜想值'''
    while not close(guess):
        guess = update(guess)
    return guess
    
def newton_update(f, df):
    '''返回牛顿迭代的计算函数, f为原函数，df是对其求导'''
    def update(x):
        return x - f(x) / df(x)
    return update
    
def newton_find_zero_of_f(f, df, tolerance=1e-13):
    '''找出f(x)=0时，x的值，tolerance为与目标真实值的误差'''
    def near_zero(x):
        return is_eq(f(x), 0, tolerance)
    return improve_guess(newton_update(f, df), near_zero)
    
def is_eq(x, y, tolerance):
    '''x与y的差值是否在误差范围内'''
    return abs(x - y) < tolerance
    
def nth_root_of_a(n, a):
    def f(x):
        return x ** n - a
    def df(x):
        return n * x ** (n - 1)
    return newton_find_zero_of_f(f, df)
    
print(nth_root_of_a(1, 64))  #64.0
print(nth_root_of_a(2, 64))  #8.0
print(nth_root_of_a(3, 64))  #4.0
print(nth_root_of_a(4, 64))  #2.82842712474619
print(nth_root_of_a(5, 64))  #2.29739670999407
print(nth_root_of_a(6, 64))  #2.0

这里是在线可运行版本 ，把每一步都抽象出来，哈哈哈哈哈，这看上去很帅，虽然实际上没啥必要，最后也可以再curry化。这是个极端的例子，但是根据业务的需要，我们可以把需要复用的步骤抽象出来，而抽象的程度也要根据业务开展。

4、牛逼哄哄的invsqrt求平方根倒数

吹水时间又到了，究竟这个牛顿迭代法有个毛线用。那你就要google下invsqrt这个改变游戏史进程的函数，可以看下知乎上的这个回答 有哪些算法惊艳到了你？。这个函数到底是干嘛的，就是求某个数的2次方根的倒数。因为这个函数在3d引擎的中经常使用，所以提高这个函数的效率是至关重要的，这个函数比系统的1/sqrt(x)快了一半。虽然是倒数，但是也是能用牛顿迭代法的啊。为毛啊，因为这也是曲线啊，过程跟1中的图相识。

但是！！！！请看代码

 

float InvSqrt (float x){
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f - xhalf*x*x);
    return x;
}

第六行就是利用牛顿迭代法的值x - f(x)/f'(x)演算后的结果，在这里f(x) = 1/x²，则f'(x)=2/x³。
第三行是把一个float指针强转成int指针，从而进行右移运算，因为float类型不能进行位运算，从而得到x值指数部分/2的结果，额，这里我不打算展开说啊，因为我说不清楚啊啊啊啊，我引用一段别人的解释，详情请移步这里：

所以，32位的浮点数用十进制实数表示就是：M2E。开根然后倒数就是：M(-1/2)2^(-E/2)。现在就十分清晰了。语句i> >1其工作就是将指数除以2，实现2(E/2)的部分。而前面用一个常数减去它，目的就是得到M(1/2)同时反转所有指数的符号。

至于0x5f3759df，这就相当无解了。wiki的介绍中，这一句的注释是这样的：

 

    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?（这他妈的是怎么回事？）

关于这个值后来也有人专门去推理过，还写成论文，但是，我肯定是看不懂的！！！
因为这一步，取到一个非常合适的初始值，只需要一次牛顿迭代就能知道目标值啦，一次！！！！

最后

没有最后，今天吹水结束！！！！！！！！