Python 实现丘德诺夫斯基（Chudnovsky）法計算高精度圓周率

# -*- coding: UTF-8 -*-

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

# Calculating PI with Chudnovsky-Series

# Author: Idealguy，2018, Shanghai

#

import time

# In following functions, High-Prec Nums are both amplified 10**n

# pre-defined: Base=10**n

#

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=()

# 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 秒" % (()-t))

s=input('是否显示结果(Y/N):')

if (s=='Y')|(s=="y"): print ("PI="+str(p))

# end