ガウス＝ルジャンドル法に基づく円周率の計算プログラム(Python版)

Gauss–Legendre.py

# -*- coding:utf-8 -*-
import math
import time 
from decimal import *
N = 3 # 試行回数

a = []
b = []
t = []
p = []

pi = 0.0;

a.insert(0,1.0)
b.insert(0,1/math.sqrt(2))
t.insert(0,0.25)
p.insert(0,1.0)

# N回の試行にかかる時間を計測
start_time = time.clock() 

# N回の試行を開始
for i in range(1, N):
    a.insert(i,(a[i-1]+b[i-1])/2)
    b.insert(i,math.sqrt(a[i-1]*b[i-1]))
    t.insert(i,(t[i-1] - (p[i-1]*(a[i-1]-a[i])*(a[i-1]-a[i]))))
    p.insert(i,2*p[i-1])

#円周率を計算
pi = ((a[N-1]+b[N-1])*(a[N-1]+b[N-1]))/(4*t[N-1])
#実行時間を計算
end_time = time.clock() 
time = end_time - start_time

#　実行結果をUNIX上に表示する
print u"円周率の近似値:%f" % (pi)
print u"実行時間:%f" % (time)