平方根倒数速算法

今天看《python科学计算》学习到的平方根倒数速算法这里实践并记录：

由于刚刚接触py不久，首先查了一下自定义函数用timeit模块计算时间性能的用法如下，

import random
import timeit
from time import clock


def get_random_number(num):
    '''get random number, no repeated element; use random.sample() method'''
    return random.sample(range(num), num)


if __name__ == "__main__":
    #use clock() method to calculate time
    start = clock()
    list_a = get_random_number(200)
    finish = clock()
    print(finish - start)
    #check the length of list generated by function
    print(len(list_a))
    print(len(set(list_a))) 

    #use timeit.Timer() method
    t1 = timeit.Timer('get_random_number(200)',
                      setup="from __main__ import get_random_number")
    #only excute once
    print(t1.timeit(1))
    #only repeat once, and only excute once
    print(t1.repeat(1, 1))

    #use timeit.Timer() and lambda to invoke function
    print(timeit.Timer(lambda: get_random_number（200）).timeit(1))

再开始做普通二分法实现的rsqrt()，及其时间效率，

def rqrt_bisection(n):
    if n < 0:
        return n
    low = 0.0
    up = n
    mid = (low+up)/2
    last = 0.0
    while abs(mid-last)>0.001:
        if mid*mid>n:
            up = mid
        else:
            low = mid
        last = mid
        mid = (low+up)/2
    return 1/mid

print rqrt_bisection(16.0)
if __name__ == '__main__':
    import timeit
    print timeit.Timer(lambda: rqrt_bisection(16.0)).timeit(10000)

精度在0.001时的10000次计算耗时26.5ms。接下来尝试牛顿迭代的方式，

def rqrt_newton(n):
    res = n
    last = (res + n/res)/2
    while abs(res-last)>0.001:
        last = res
        res = (res + n/res)/2
    return 1/res

print rqrt_newton(16.0)
if __name__ == '__main__':
    import timeit
    print timeit.Timer(lambda: rqrt_newton(16.0)).timeit(10000)

耗时12.4ms有了较大的提升。这种方法的原理很简单，我们仅仅是不断用(x,f(x))的切线来逼近方程x^2-a=0的根。根号a实际上就是x^2-a=0的一个正实根，这个函数的导数是2x。也就是说，函数上任一点(x,f(x))处的切线斜率是2x。那么，x-f(x)/(2x)就是一个比x更接近的近似值。代入 f(x)=x^2-a得到x-(x^2-a)/(2x)，也就是(x+a/x)/2。但是依旧与pythonmath带的sqrt()差距很大，

def rqrt_py(n):
    return 1/math.sqrt(n)

print rqrt_py(16.0)
if __name__ == '__main__':
    import timeit
    print timeit.Timer(lambda: rqrt_py(16.0)).timeit(10000)

耗时只有2.3ms。差距在哪里呢？维基百科关于《雷神之锤》中使用0x5f3759df计算平方根倒数的解释告诉我牛逼的算法，

def rqrt(n):
    number = np.array([n])
    y = number.astype(np.float32)
    x2 = y*0.5
    i=y.view(np.int32)
    i[:]=0x5f3759df-(i>>1)
    y = y*(1.5-x2*y*y)
    return y

此算法首先接收一个32位带符浮点数，然后将之作为一个32位整数看待，以将其向右进行一次逻辑移位的方式将之取半，并用十六进制“魔术数字”0x5f3759df减之，如此即可得对输入的浮点数的平方根倒数的首次近似值；而后重新将其作为浮点数，以牛顿法反复迭代，以求出更精确的近似值，直至求出符合精确度要求的近似值。在计算浮点数的平方根倒数的同一精度的近似值时，此算法比直接使用浮点数除法要快四倍。实验显示耗时2.53ms略比pythonmath的1/sqrt()慢一点点，应该是写法上用numpy的array导致的吧。