因为吹水的能力不佳,所以要先打个草稿,今天的吹水过程大概是:
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?(这他妈的是怎么回事?)
关于这个值后来也有人专门去推理过,还写成论文,但是,我肯定是看不懂的!!!
因为这一步,取到一个非常合适的初始值,只需要一次牛顿迭代就能知道目标值啦,一次!!!!
最后
没有最后,今天吹水结束!!!!!!!!