二分法和牛顿迭代法

先说一个面试题：问 1.2 - 0.2  ==  1 ？

　　答案是False！ 为什么？

其原因在于十进制和二进制的转换上，计算机先要把十进制的数转化为二进制，然后再计算。
但是，在转化中，浮点数转化为二进制，就出问题了，例如：
十进制的 0.1，转化为二进制是：0.0001100110011001100110011001100110011001100110011…（不能精确）
也就是说，转化为二进制后，不会精确等于十进制的 0.1。
同时，计算机存储的位数是有限制的，所以，就出现上述现象了。

这种问题不仅仅是 Python 中有，所有支持浮点数运算的编程语言都会遇到，它不是 Python 的 bug 

 

求一个数的平方根函数sqrt(int num)，在大多数语言中都提供实现。那么要求一个数的平方根，是怎么实现的哪？

　　实际上求平方根的算法方法主要有两种：二分法( binary search)和牛顿迭代法( Newton iteration):

说到这两种方法，心里就感慨重重，记得大学开了一门叫《计算方法》的课程，后来就忘了，直到有一天老师讲到

这里的时候才想起这种计算方法自己曾几何时学过！

1、二分法：（夹逼法）　　

　　以求根号5为例：

　　a:折半：       5/2=2.5

　　b:平方校验:  2.5*2.5=6.25>5，并且得到当前上限2.5

　　c:再次向下折半:2.5/2=1.25

　　d:平方校验：1.25*1.25=1.5625<5,得到当前下限1.25

　　e:再次折半:2.5-(2.5-1.25)/2=1.875

　　f:平方校验：1.875*1.875=3.515625<5,得到当前下限1.875

每次得到当前和5进行比较，并且记下下限和上线，依次迭代，逐渐逼近我们设好的误差范围：

import math
from math import sqrt
 
def sqrt_binary(num):
    x=sqrt(num)    # 求出系统给我们的准确值
    y=num/2.0　　　　#中分
    low=0.0
    up=num*1.0
    count=1
    while abs(y-x)>0.00000001:
        print count,y          
        count+=1           
        if (y*y>num):
            up=y           #二分后中间值比预期大了
            y=low+(y-low)/2
        else:
            low=y
            y=up-(up-y)/2   #二分后中间值比预期小了
    return y
 
print(sqrt_binary(5))
print(sqrt(5))

运行结果:

　　1 2.5
　　2 1.25
　　3 1.875
　　4 2.1875
　　5 2.34375
　　6 2.265625
　　7 2.2265625
　　8 2.24609375
　　9 2.236328125
　　10 2.2314453125
　　11 2.23388671875
　　12 2.23510742188
　　13 2.23571777344
　　14 2.23602294922
　　15 2.23617553711
　　16 2.23609924316
　　17 2.23606109619
　　18 2.23608016968
　　19 2.23607063293
　　20 2.23606586456
　　21 2.23606824875
　　22 2.23606705666
　　23 2.2360676527
　　24 2.23606795073
　　25 2.23606809974
　　26 2.23606802523
　　27 2.23606798798
2.23606796935
2.2360679775

经过27次二分迭代，得到的值和系统sqrt()差别在0.0000001，

2、牛顿迭代法：

仔细思考一下就能发现，我们需要解决的问题可以简单化理解。

从函数意义上理解：我们是要求函数f(x) = x²，使f(x) = num的近似解，即x² - num = 0的近似解。

从几何意义上理解：我们是要求抛物线g(x) = x² - num与x轴交点（g(x) = 0）最接近的点。

 

我们假设g(x0)=0，即x0是正解，那么我们要做的就是让近似解x不断逼近x0，这是函数导数的定义

 可以由此得到： (可以写出f'(x)的导数形式）   

 

从几何图形上看，因为导数是切线，通过不断迭代，导数与x轴的交点会不断逼近x0。

所以对于一般情况：

将m=2带入得：

代码如下：

def sqrt_newton(num):
    x=sqrt(num)
    y=num/2.0
    count=1
    while abs(y-x)>0.00000001:
        print count,y
        count+=1
        y=((y*1.0)+(1.0*num)/y)/2.0000
    return y
 
print(sqrt_newton(5))
print(sqrt(5))

运算结果： 

 1 2.5
 2 2.25
 3 2.23611111111

对于牛顿迭代法来说：还可以求立方，甚至等多（只需要修改m得值即可）

def cube_newton(num):
    x=num/3.0
    y=0
    count=1
    while abs(x-y)>0.00000001:
        print count,x
        count+=1
        y=x
        x=(2.0/3.0)*x+(num*1.0)/(x*x*3.0)
    return x
 
print(cube_newton(27))