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

idealguy 2018-10-02 10:50:39

本文来自 idealguy 的CSDN 博客，全文地址请点击
https://blog.csdn.net/idealguy/article/details/82929032?utm_source=copy

【原创】

本文利用Python自带的高精度整型计算功能，实现丘德诺夫斯基（Chudnovsky）法計算高精度圓周率。计算速度比马青公式要快。

经过艰苦的调试，实现该计算，10万位的计算在秒级完成。以下是计算截图：

以下是代码：



# -*- 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=time.clock()

    # 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.clock()-t))

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

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

# end