Cython(4) 円周率の計算

コンピュータで延々と計算させる例として、円周率を何万桁、何億桁も求めるのがあるが、ここでは、PythonとCythonの浮動小数点演算の速度比較のために行う。計算方法は色々あるが、桁数の新記録樹立とかではなく、あくまで速度比較なので、非常に簡単で有名なライプニッツの公式（マーダヴァ・ライプニッツの公式）を用いる。この公式は、以下のように表される（Wikipediaより）。

この級数は、収束がとても遅く、円周率を計算するには不適切であるが、収束が遅いために速度比較には最適なので、これを選んだ。さらに、これが一番単純な円周率を計算する公式だという理由もある。実際に計算に利用されるのは、収束が非常に速い公式であり、それらはWikipediaの円周率に関する式を参照されたい。

## [Python] 円周率　マーダヴァ・ライプニッツの公式
def pi_leibniz(n):
    sum = 0.0    
    for k in range(1,n,4):
        sum += 4/k - 4/(k+2)
    return sum

２項ずつまとめて計算するようにしている。こうすると、sumは単調増加で次第に円周率に近づいていく。
nの値により打ち切る項が決まるが、n/2項までの和になる。

テストするプログラムは、n(項数の2倍)を10から10倍ずつ増やしながら表示し、システムの提供する円周率との差と計算時間を示した。

import time
import math
from pi_leibniz import pi_leibniz

for m in range(1,10):
	n = 10**m
	tm_start = time.time()
	pi = pi_leibniz(n)
	tm_end = time.time()
	print("{:12d}\t{}\t{:.16f}\t{:15.10f} sec".format(n,pi, pi-math.pi, tm_end-tm_start))

実行結果は以下のようになった。

$ python3 pi_leibniz_test.py
     	  10	2.9760461760461765  	-0.1655464775436166    	0.0000028610 sec
     	 100	3.1215946525910097  	-0.0199980009987835    	0.0000040531 sec
    	1000	3.1395926555897824  	-0.0019999980000107    	0.0000345707 sec
   	   10000	3.1413926535917893  	-0.0001999999980038    	0.0003449917 sec
  	  100000	3.141572653589808   	-0.0000199999999850    	0.0034744740 sec
 	 1000000	3.1415906535898936  	-0.0000019999998995    	0.0255348682 sec
	10000000	3.14159245358981    	-0.0000001999999832    	0.2365379333 sec
   100000000	3.1415926335405047  	-0.0000000200492885    	2.3944559097 sec
  1000000000	3.1415926445762157  	-0.0000000090135774   24.2318170071 sec

Python版の関数に型宣言を加えてCython版を作った。

## [Cython] 円周率　マーダヴァ・ライプニッツの公式
def pi_leibniz(int n):
    cdef double sum = 0.0
    cdef int   k
    
    for k in range(1,n,4):
        sum += 4/k - 4/(k+2)
    return sum

setup.pyにpi_leibniz.pyxを書き加えてコンパイルし実行すると、こうなった。

$ python3 pi_leibniz_test.py
      	  10    2.9760461760461765    -0.1655464775436166   	0.0000009537 sec
     	 100    3.1215946525910097    -0.0199980009987835   	0.0000004768 sec
    	1000    3.1395926555897824    -0.0019999980000107   	0.0000011921 sec
   	   10000    3.1413926535917893    -0.0001999999980038   	0.0000083447 sec
      100000    3.141572653589808     -0.0000199999999850    	0.0000808239 sec
 	 1000000    3.1415906535898936    -0.0000019999998995   	0.0008053780 sec
	10000000    3.14159245358981      -0.0000001999999832   	0.0070176125 sec
   100000000    3.1415926335405047    -0.0000000200492885   	0.0622570515 sec
  1000000000    3.1415926445762157    -0.0000000090135774   	0.5700106621 sec

計算結果は、PythonとCythonで完全に一致している。これは、Pythonの浮動小数点演算がCのdouble型と同じ精度で行われていることを示していると言えよう。

速度であるが、24.23/0.57 = 42.5倍速になった。まずまずの高速化で、Cython化することで遅くて使えなかったプログラムが実用に耐えるレベルになることがあると考えられる。
時間計測に、time.time()を使っているが、どうやら精度はマイクロ秒が限界のようである。より詳細な時間（短い時間）を計測するためには、time.time_ns()が用意されているが、その紹介はまたの機会に。