平方根倒数速算法

来源：互联网发布：程序员接私活app 编辑：程序博客网时间：2024/04/27 19:28

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

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

 1 import random 2 import timeit 3 from time import clock 4  5  6 def get_random_number(num): 7     '''get random number, no repeated element; use random.sample() method''' 8     return random.sample(range(num), num) 9 10 11 if __name__ == "__main__":12     #use clock() method to calculate time13     start = clock()14     list_a = get_random_number(200)15     finish = clock()16     print(finish - start)17     #check the length of list generated by function18     print(len(list_a))19     print(len(set(list_a))) 20 21     #use timeit.Timer() method22     t1 = timeit.Timer('get_random_number(200)',23                       setup="from __main__ import get_random_number")24     #only excute once25     print(t1.timeit(1))26     #only repeat once, and only excute once27     print(t1.repeat(1, 1))28 29     #use timeit.Timer() and lambda to invoke function30     print(timeit.Timer(lambda: get_random_number（200）).timeit(1))

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

 1 def rqrt_bisection(n): 2     if n < 0: 3         return n 4     low = 0.0 5     up = n 6     mid = (low+up)/2 7     last = 0.0 8     while abs(mid-last)>0.001: 9         if mid*mid>n:10             up = mid11         else:12             low = mid13         last = mid14         mid = (low+up)/215     return 1/mid16 17 print rqrt_bisection(16.0)18 if __name__ == '__main__':19     import timeit20     print timeit.Timer(lambda: rqrt_bisection(16.0)).timeit(10000)

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

 1 def rqrt_newton(n): 2     res = n 3     last = (res + n/res)/2 4     while abs(res-last)>0.001: 5         last = res 6         res = (res + n/res)/2 7     return 1/res 8  9 print rqrt_newton(16.0)10 if __name__ == '__main__':11     import timeit12     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计算平方根倒数的解释告诉我牛逼的算法，

1 def rqrt(n):2     number = np.array([n])3     y = number.astype(np.float32)4     x2 = y*0.55     i=y.view(np.int32)6     i[:]=0x5f3759df-(i>>1)7     y = y*(1.5-x2*y*y)8     return y

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

0 0