算π

 

# 丘德诺夫斯基法計算高精度圓周率程序

import time 

def Sqrt10005(): # Sqrt(10005L) by Imitate-Manual Method
    n1=0
    c=10002499687 #100.02499687
    mc=8; m=mc
    f1=10**mc
    f2=f1*f1
    a=10005*f2-c*c
    while mc<n:
        a*=f2
        b=c*2*f1
        d=a//b
        c=c*f1+d
        a-=d*(b+d)
        mc+=m
        if mc*2>n: m=n-mc
        else: m=mc
        f1=10**m
        f2=f1*f1
        n1+=1
    return c

# Main Program
#
print ("Chudnovsky法計算高精度圓周率程序")
while 1:
    n=int(input('計算位數[1..50000],0:退出：'))
    if n<=10: break
    n+=2
    base=10**n
    t=time.perf_counter ()
# Start Calculating
    A=13591409*base; B=A
    c3=13591409
    i=1
    while abs(A)>5:
        c1=((108-72*i)*i-46)*i+5
        c2=10939058860032000*i**3
        c4=c3; c3+=545140134
        i+=1
        A=A*c1*c3//(c2*c4) # Must in form: A=A*...
        B+=A
    p=426880*base*Sqrt10005()//B//100
# Post access
    print ("用時= %8.3f 秒" % (time.perf_counter ()-t))
    print ("PI="+str(p))

# end