今天看《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导致的吧。