猫史档案馆

【圆周率】又一个神奇的计算圆周率的程序 能计算到小数点后数十万位

用户:SunshineCodeSunshineCode查看:171 回复:8 评论:13 创建时间:2019-06-05T20:20:23

最近搜到一个C语言程序，是这样的：

#include<stdio.h>
#include<iostream>


using namespace std;

int main(void)

{

    printf("本程序每四位数输出结果，如果请求计算的位数不是4的整数倍，最后输出可能会少1~3位");

    long a[2]={956,80},b[2]={57121,25},i=0,j,k,p,q,r,s=2,t,u,v,N,M=10000;
    printf("计算公式：Machin%4cpi=16arctan(1/5)-4arctan(1/239)\n计算到的位数？\n",32,32);
    cin>>N,N=N/4+3;
    long *pi=new long[N],*e=new long[N];
    while(i<N)pi[i++]=0;
    while(--s+1)
    {
        for(*e=a[k=s],i=N;--i;)e[i]=0;
        for(q=1;j=i-1,i<N;e[i]?0:++i,q+=2,k=!k)
        for(r=v=0;++j<N;pi[j]+=k?u:-u)u=(t=v*M+(e[j]=(p=r*M+e[j])/b[s]))/q,r=p%b[s],v=t%q;
    }
    while(--i)(pi[i]=(t=pi[i]+s)%M)<0?pi[i]+=M,s=t/M-1:s=t/M;
    for(cout<<"3.";++i<N-2;)printf("%04ld",pi[i]);
    delete []pi,delete []e,cin.ignore(),cin.ignore();
    return 0;
}

 

这个程序可以计算制定位数的圆周率值，

你想计算到多少位，只要你有耐心，电脑有足够的运行内存，都是可以算出来的。

迭代速度很快，计算到小数点后20000位，只用了两三秒（xp32位机）

所用的计算公式是马青公式：

pi = 16 arctan 1/5 + 4 arctan 1/239

*其中arctan又可写作arctg

arctan是反正切，计算公式为：

arctan x = x - 1/3 x^3 + 1/5 x^5 - 1/7 x^7 + 1/9 x^9 - ...

 

这个公式计算圆周率收敛速度很快，每迭代一次就能计算出1.39位的圆周率值。

于是，我最近就在据这个程序研究如何在编程猫中实现。

回复

上一页第 1 页 / 共 1 页下一页
SunshineCodeSunshineCode

注："^"表示次方

点赞0

评论

SunshineCodeSunshineCode

补：原C程序改写Python：

def show(n,plt_type='pie'):

result = f(n)

a = pd.Series([int(i) for i in result])

b = a.value_counts()

if plt_type=='pie':

plt.figure(figsize=(10,8))

plt.pie(b,labels=b.index,explode=[0.05]*10,shadow=True,autopct='%1.1f%%')

plt.show()

elif plt_type=='bar':

plt.figure(figsize=(10,8))

plt.bar(b.index,b.values)

plt.show()

else:

print('Type Wrong')

 

点赞0

评论

FATAL_ERRORFATAL_ERROR

这么好的东西

没人顶顶可惜了emotion_编程猫_点赞

点赞0

评论

急_开_锁_办_证810864急_开_锁_办_证810864

少侠，这是你要的编程猫实现。编程猫没有函数返回值，要实现高级一点的功能复用的话，都要靠全局变量传递结果，很难受。center_image

其实对于演示圆周率的近似计算，Nilakantha级数和著名的Basel问题(Leonhard Euler最早给出证明的那个)都更适合用于在Scratch类的软件中练习。

点赞0

评论

输出更多位数的圆周率的值，算法确实很复杂，数学不好的话一不小心就被绕晕了。不过这种算法，在很多时候都是有价值的，因此我认为这个算法值得研究。

圆周率是常数，不需要在运行的时候临时用算法去算；提前算到一定的精度，然后定义为常数直接使用就行了。必须进行即时演算的事情才需要进行算法设计和优化。有关圆周率计算方法的研究，就像总有人去挑战圆周率的背诵一样，更像是一种智力游戏。

不过你可以把这个单纯计算圆周率的程序变得更有趣一些：请用户随意输入1个长度有限的数字序列，然后通过程序搜寻，“证明”圆周率中总能找出这样的连续数字序列。比起圆周率能算多快，数学上还有很多更重要的关于圆周率的谜没有解开，这个就是其中之一：目前为止还没有任何人能证明圆周率是一个合取数，但你能编写程序验证例举并验证足够多的个例。

这个已经证明了，圆周率是一个正规数，但是这个寻找数串的前提是你有圆周率的值

圆周率的数字规律和质数分布一样都是世界难题，历史上进展极其缓慢，应该没那么容易证明。计算的话，即使不用级数，也可以用蒙特卡洛方法去随机撒足够多的点，看看有多少能落在单位圆内。当然这些算法很慢而已，不是只能依赖级数。17年高教社出版的《Python语言程序设计基础》(第2版)第4章第6小节里有提到“迄今为止求解圆周率最好的方法是利用BBP公式，该公式如下...”，不过谁知道呢。

啊不嘟啊不嘟

####################导入时间模块
import time
###############计算当前时间
time1=time.time()
################算法根据马青公式计算圆周率####################


number = int(input('请输入想要计算到小数点后的位数:'))
print('计算中......')

# 多计算10位，防止尾数取舍的影响
number1 = number+10

# 算到小数点后number1位
b = 10**number1
# 求b*4/5的首项
x1 = b*4//5
# 求b/-239的首项
x2 = b//-239

# 求第一大项
he = x1+x2
#设置下面循环的终点，即共计算n项
number *= 2

#循环初值=3，末值2n,步长=2
for i in range(3,number,2):
    # 求每个含1/5的项及符号
    x1 //= -25
    # 求每个含1/239的项及符号
    x2 //= -57121
    # 求两项之和
    x = (x1+x2) // i
    # 求总和
    he += x

#求出π
    pi = he*4
#舍掉后十位
    pi //= 10**10


#输出圆周率π的值
f = open('pi.txt','w')
paistring=str(pi)
result=paistring[0]+str('.')+paistring[1:len(paistring)]
f.write(result)
f.close()
time2=time.time()
print('π的值已存储到当前目录的pi.txt文件中')
print ('总共耗时：'+ str((time2 - time1)//1) + '秒')

python版本

点赞1

评论

啊不嘟啊不嘟

源码编辑器版

点赞0

评论

由于kitten的bug只能算到16位

其实用kitten也是可以算到更多位数的，只不过算法有点复杂。我准备暑假研究能算出更多位数的，比如算个200位。 输出更多位数的值，关键在于算法。 C语言也是只能保存16个小数位的，但是这个程序能计算出数十万位的圆周率的值，就在于它的算法。算法是怎样的，我正在研究。

kitten会出现E数

我python算了十万位，还可以继续算，花了40秒

存在了一个txt文件里

我易语言用很复杂的算法，丘德诺夫斯基公式，算圆周率算10万位0.8秒，10亿位80分钟

SunshineCodeSunshineCode

分享一下：现在我们研究计算圆周率，都是用的马青公式，因为这个公式比较简单，容易实现，不过这个公式收敛速度较慢，每迭代一次只能算出1.39位的圆周率值。如果大家想要超高精度超高速度的圆周率计算，不妨使用丘德诺夫斯基公式。这个公式每迭代一次就能计算出14位的圆周率值。

代码如下（python）：

# -*- 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

点赞0

评论

没有缩进好难受

啊不嘟啊不嘟

# -*- 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))
    with open('pi.txt','w') as f:
        f.write('3.'+str(p)[1:n-1])
    print('pi的值已保存到当前目录：pi.txt')
# end

 

点赞0

评论

加个缩进顺便优化下