计算圆周率PI的算法

时间：2021-11-05 18:34:53
现代的计算机可以把圆周率计算到数亿位,请问有没有人知道这是如何实现的?
类似的还有自然对数底e

期待高手,解决立马给分!

44 个解决方案

#1

关注

#2

你学过高等数学吗？
得到pi 的公式多的很，同让e也一样。
现在关键是要找到一个收敛速度快的函数来计算。

#3

http://sq.k12.com.cn:9000/tg/tgq163/sxyzl.htm

看看这篇文章，他的内容不是最新的了，但是值得参考。

#4

数学是很重要的!

#5

利用arctg45=1来计算

把arctgx展开成无穷级数，所求结果是这个无穷级数的4倍

#6

这问题有点意思，楼主你怎么想到问这样的问题的？

#7

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    Label1: TLabel;
    Label2: TLabel;
    Edit1: TEdit;
    Edit2: TEdit;
    procedure Edit1KeyPress(Sender: TObject; var Key: Char);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Edit1KeyPress(Sender: TObject; var Key: Char);
var
  k:integer;
  p,t,s,e:real;
  m:string;
begin
  if key=#13 then
    begin
      k:=strtoint(edit1.Text);
      if k>50 then
        begin
          showmessage('数值太大可能会使程序失去响应');
          edit1.Text:='18';
          k:=18;
        end;
      p:=0; s:=2;
      e:=exp(k*ln(0.1));
      repeat
        t:=s;
        p:=sqrt(2+p);
        s:=s*2/p;
      until abs(t-s)<e;
      m:=format('%20.18f',[s]);
      edit2.Text:=copy(m,0,k+2);
      edit1.SelStart:=0;
      edit1.SelLength:=length(edit1.Text);
      edit1.SetFocus;
    end;
end;

end.
说明一下：我是菜鸟，刚开始学delphi。别骂我啊！这是我书上的一道例题。

#8


我听说有人按 PI 的数位编出乐谱，再用计算机演奏出来很好听的

大家有兴趣可以做一个啊
做出来别忘了给我一份  ljlblue@163.com

#9

long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(;b-c;)
  f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
  for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

#10

你自己改成Pascal的
Okay？

#11

up

#12

这个程序我上学时做过习题,现在源码没有了,不过可以提两点建议.
1.直接按arctg45的级数展开是不可取的,因为它收敛太慢,正确的级数是两个反正切的差,具体我记不清了,但是你可以在<初中代数的教学参考书中找到>,我当时是在那里找的.
2.必须定义一种可以表示任意精度(也可说是任意长度整数)的数据结构,计算机本身是不提供的,需要你去设计.还要能过进行加\减\乘\除.实现时你可以研究一下汇编,看计算机是怎么实现的,然后将这种计算扩展到任意长.我当时做的就在Dephi中嵌的汇编.

完成这两件事,你的题就完成了
估计要几个小时的时间才行.

#13

我说的几小时是你编程的时间,不是计算机计算的时间
我算到2000位也用不了几秒种

自然对数e简单些,直接展开收敛也很快

#14


to： mingyeh(插死陈文登) 

   太历害了！佩服！
   短短两行代码可以算出那么多位PI！

   可是不知道你是依据什么公式计算的，什么原理，能达到多少位的精度
   望赐教！

#15


to: wan_ming(wan ming) 
*********************************************************

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    Label1: TLabel;
    Label2: TLabel;
    Edit1: TEdit;
    Edit2: TEdit;
    procedure Edit1KeyPress(Sender: TObject; var Key: Char);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Edit1KeyPress(Sender: TObject; var Key: Char);
var
  k:integer;
  p,t,s,e:real;
  m:string;
begin
  if key=#13 then
    begin
      k:=strtoint(edit1.Text);
      if k>50 then
        begin
          showmessage('数值太大可能会使程序失去响应');
          edit1.Text:='18';
          k:=18;
        end;
      p:=0; s:=2;
      e:=exp(k*ln(0.1));
      repeat
        t:=s;
        p:=sqrt(2+p);
        s:=s*2/p;
      until abs(t-s)<e;
      m:=format('%20.18f',[s]);
      edit2.Text:=copy(m,0,k+2);
      edit1.SelStart:=0;
      edit1.SelLength:=length(edit1.Text);
      edit1.SetFocus;
    end;
end;

end.

这个算法的依据是:
   置初值  p=0,s=2
      重复计算：
          p=sqrt(2+p);    s=s*2/p;
   每重复一次 s 越接近于 PI

但这个算法最多只能达到高级语言本身数据类型(如 c 的 double )的精度，达不到我提出的要求

#16

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char * argv[])
{
  long * pi, * t, m, n, r, s;
  int t0[][3] = {48, 32, 20, 24, 8, 4}, k0[][3] = {1, 1, 0, 1, 1, 1};
  int n0[][3] = {18, 57, 239, 8, 57, 239}, d, i, j, k, p, q;

  d = (argc > 1) ? (((i = atoi(argv[1])) < 0) ? 0 : i) : 1000;
  q = (argc > 2) ? 1 : 0;
  printf("pi= %s%d * arctg(1/%d) %s %d * arctg(1/%d) %s %d * arctg(1/%d) [%s]\n",
    k0[q][0] ? "" : "-", t0[q][0], n0[q][0], k0[q][1] ? "+" : "-", t0[q][1],
    n0[q][1], k0[q][2] ? "+" : "-", t0[q][2], n0[q][2], q ? "Stomer" : "Gauss");
  if ((t = (long *)calloc((d += 5) + 1, sizeof(long))) == NULL) return 1;
  if ((pi = (long *)calloc(d + 1, sizeof(long))) == NULL) return 2;
  for (i = d; i >= 0; i--) pi[i] = 0;
  for (p = 0; p < 3; p++) {
    for (k=k0[q][p], n=n0[q][p], t[i=j=d]=t0[q][p], i--; i >= 0; i--) t[i] = 0;
    for (r = 0, i = j; i >= 0; i--) {
      r = (m = 10 * r + t[i]) % n;
      t[i] = m / n;
      k ? (pi[i] += t[i]) : (pi[i] -= t[i]);
    }
    while (j > 0 && t[j] == 0) j--;
    for (k = !k, s = 3, n *= n; j > 0; k = !k, s += 2) {
      for (r = 0, i = j; i >= 0; i--) {
        r = (m = 10 * r + t[i]) % n;
        t[i] = m / n;
      }
      while (j > 0 && t[j] == 0) j--;
      for (r = 0, i = j; i >= 0; i--) {
        r = (m = 10 * r + t[i]) % s;
        m /= s;
        k ? (pi[i] += m) : (pi[i] -= m);
      }
    }
  }
  for (n = i = 0; i <= d; pi[i++] = r) {
    n = (m = pi[i] + n) / 10;
    if ((r = m % 10) < 0) r += 10, n--;
  }
  printf("pi= %ld.", pi[d]);
  for (i = d - 1; i >= 5; i--)
    printf("%ld%s", pi[i], ((m = d - i + 5) % 65) ? ((m % 5) ? "" : " ") : "\n");
  printf("%sDIGITS: %d\n", (m % 65) ? "\n" : "", d - 5);
  return 0;
}

#17

!!!!!!!!!!!!!!!!
long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(;b-c;)
  f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
  for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

#18

上面的代码我也觉得奇怪，b没附初值这段这么来的，for(;b-c;)
  f[b++]=a/5;
这也能编译通过？？

#19

#define N 2800
long a=10000,b,c=N,d,e,f[N+1],g;
main()
{
while(b<c)
  f[b++]=a/5;

while(c>0){
    d=0; g=c*2; b=c;
    do{
       d+=f[b]*a; f[b]=d%--g;
       d/=g--;    b--; 
       if(b) d*=b; else break;
    }while(1);
    printf("%.4d",e+d/a);
    e=d%a;
    c-=14;
}
}

从 mingyeh(插死陈文登) 改进的（可读性好一些），VC 6 编译通过，没有问题，int 型默认初值为 0

#20

我试了，dcyu(Dd）和 mingyeh(插死陈文登) 的算法在 10000 位以内都是对的，10000以后就不知道了

#21

拉马努金的改进公式：
设f(n)=(13591409+545140134n)/640320^(3*n)*(-1)^n*(6n)!/(n!)^3/(3n)!
则1/pi=12/640320^(3/2)*(f(0)+f(1)+.......)
这个级数每增加一项能提高14位小数的精度

迭代方法：
y[0]=Sqrt(2)-1
z[n]=(1-y[n-1]^4)^(1/4)
y[n]=(1-z[n])/(1+z[n])
a[0]=6-4*Sqrt(2)
a[n]=(1+y[n])^4*a[n-1]-2^(2n+1)*y[n]*(1+y[n]+y[n]^2)
8次迭代能精确到小数点178814位

#22

计算e非常简单，而且收敛很快，用公式：
e=1+1/1!+1/2!+1/3!+1/4!+...

#23

我在书店里看到一本小册子专门讲如何算pi的

#24

long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(;b-c;)
  f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
  for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

这种算法的数据理论依据, 请参见这篇文章:
Stanley Rabinowitz, Stan wagon
A Spigot Algorithm for the Digits of pi
American Mathematical Monthly, 
vol. 102, Iss. 3 (Mar. 1995), pp195-203.

该期刊电子版网址为:
http://uk.jstor.org/journals/maa.html

==================================================
有关pi值计算的一些算法, 也可参见:
http://numbers.computation.free.fr/Constants/TinyPrograms/tinycodes.html
http://hades.ph.tn.tudelft.nl/Internal/PHServices/Documentation/MathWorld/math/math/p/p303.htm
http://www1.physik.tu-muenchen.de/~gammel/matpack/html/Mathematics/Pi.html

#25

呵呵，我有原代码，可是只能算几万的，算不了几亿的，算几亿靠我的
塞杨二566恐怕不行，！！！
如果只是几百万的话，还没问题，到几亿，而你又不想等的话，就必须靠
fft（算是FFT的变种吧！）具体贴起来可就多了。

#26

有二次收敛算法，就是如果这次精度是1000位，下一次计算可以达到2000位精度，就是比较复杂，要是几亿位，恐怕编程的时间要远远超过最后算出PI的时间.
sqrt()表示平方根运算
x0=sqrt(2)
pi0=2+sqrt(2)
y0=sqrt(sqrt(2))

x(i+1)=1/2*(sqrt(x(i)) + 1/sqrt(x(i)))
pi(i+1)=pi(i)*(x(i+1)+1)/(y(i)+1)
y(i+1)=(y(i)*sqrt(x(i+1)) + 1/sqrt(x(i+1))) / (y(i) + 1)

希望没有打错。
算法用到了平方根运算，乘除法，必须运用FFT算法。

#27

e就好算多了，直接求1/n!的和就可以了

#28

还有一种用概率来求的。用园和矩形，然后随机生成点。根据点落的在矩形内的多少。然后依据几何概型来求！

#29

恐怕编程的时间要远远超过最后算出PI的时间.
呵呵，这倒不至于，不过确实非常的麻烦！
而且很多都没办法，不得不用汇编写（用C的话效率就要低一倍以上，没办法。还的稍微利用一点硬件的特性！
哎算几亿，太麻烦了，
算个2万吧，又太容易。哎。

#30

去年大一的时候，因为老师悬赏，我写了一个能计算pi的程序。
用turbo pascal 7.0写的，当时整整花了一周的时间，什么事都没干。
程序能计算到8000位，1000位的时间大概在15分钟。当然由于当时急着交作业，也就没对算法进行优化。总的思路如下：
1.要定义4个高精度的计算，加、减、乘、除。
2.然后找一个收敛比较快的公式，（现在我忘了）
3.按照公式计算，只不过所有的运算都必须是高精度运算。
由于在dos下，内存的空间有限，因此数组最多只能到8000个元素，如果在windows中用delphi或者是c++builder可以计算得更多。
我写了200行的源代码，不过我觉得这挺有意思的，至少能够说明自己编程的水平。有兴趣的话，或者需要源程序，我们可以联系.
我的电邮是supercctv@sina.com

#31

http://www.jason314.com/
这个网站把PI介绍的很详细
我用其中的Bailey-Borwein-Plouffe算法的改进算法经过优化后计算PI小数点后100,000位仅用1秒钟。
还有更优的计算公式，不过我不大会实现。

#32

以100000000为例，我说的方法内存占用大于pi的内存占用的8倍，而pi占用
100000000/8*0.3010，所以总的内存占用为330MB左右，要是没有大的内存，就免了。还有，总共需要26步迭代，每次要是用到1分钟的时间，你的程序就很快了。记住，计算乘法使用的FFT算法，最大的运算位数只是几千位，要是计算一亿位，就会需要特殊的算法.

#33

vgrtgfg

#34

参考了手边的几本书, 比较的结果是在众多的展开式中, 以下式收敛较快

pi = 24 * arctan(1/8) + 8 * arctan(1/57) + 4 * arctan(1/239)

其中, 在计算时将 arctan 函数用级数展开.

不知大家有何看法?

#35

大家不用忙着贴代码，先听我慢慢说。下面的公式一般我们叫做马庭公式：
pi = 24 * arctan(1/8) + 8 * arctan(1/57) + 4 * arctan(1/239)
如果不是要求太高，用上面的公式就行了，原因就是因为编程比较简单。去年，我曾经花了几个月去研究这个问题。我找了不少资料，参照了这个网站的信息：http://www.jason314.com/ ，后来改进了站长的汇编程序，如果没记错，速度可以在几小时内算到一百万位。但这已经是马庭公式这类型的算法的极限。我国内下载的程序也没有超越这个速度的。这类的公式还有更快的，如下，但是我想速度也不会有什么突破。
Pi/4 = 22*arctan(1/28)+2*arctan(1/443)-5*arctan(1/1393)-10*arctan(1/11018)
Pi/4 = 44*arctan(1/57)+7*arctan(1/239)-12*arctan(1/682)+24*arctan(1/12943)

#36

现在回答楼主的问题，真正能算PI的程序我只见过两个。
第一个是流行测试 CPU 的 SuperPI 程序。这个程序是由日本人用 F2C 把   Fortran 的源程序改成 C 的，改得适合在PC机上使用，使用的就是 Gauss-Legendre algorithm （高斯－勒让德）算法。而那个Fortran程序就是由近年来一直刷新世界纪录的那两个日本人  Yasumasa KANADA (Associate Professor)  和  D.Takahashi  写的，其实是刷新自己的记录而已。他们使用的是Gauss-Legendre algorithm算法和Borwein's 4-th order convergent algorithm算法，也许Fortran更适合科学计算。具体的算法我就不贴出来了，如果大家有兴趣可以找一下，简单的说就是用迭代法。关键问题已经不是寻找算法，而是怎样在更好的在计算机里实现算法，如何存储结果等。

#37

Introduction of Kanada Lab.

Field of Research:
 High Performance Computing, Virtual Workstation System, Research on Research, Computer Algebra

Name of Laboratory:
 Laboratory for Computational tools

Location:
 University of Tokyo. Information science.

Name of staff:
 Yasumasa KANADA (Associate Professor)
 Akira YOSHIOKA  (Research Associate)

网站在东京大学里，要找还是可以的。

#38

连续的回复不能超过三次！

#39

第二个是PIFast，现在最新的版本是4.1，是一个法国人写的。日本人的程序是在大型机上运行的，而PIFast完完全全是在PC机运行的，主要是用了FFT，所以速度很快，快得你不相信，比PC版的SuperPI快很多。但是具体算法就不是很清楚，大家可以到他的网站看看。
http://numbers.computation.free.fr/Constants/constants.html

#40

网址我快不记得了，不过我下载的那个程序我特意也superPI比了一下，
算的少的话，看不出差距，越大差距越大，
那个程序比super pi 快多了！
不过他用了那个Apf233，里面一堆汇编文件，我好不容易才在VC6。0下编译通
过，运行。

#41

忘了告诉大家，目前PC机上流行的最快的圆周率计算程序是PiFast，PC机上的最高计算记录12,884,901,372位。使用Chudnovsky公式计算，速度最快。还可以使用Ramanujan公式，速度慢大约40%。PiFast可以利用磁盘缓存，突破物理内存的限制进行超高精度的计算，最高计算位数可达240亿位，并提供基于Fabrice Bellard公式的验算功能。它除了计算圆周率，还可以计算e和sqrt(2)。

#42

还可以用
多特魔德法！！！！

#43

说 错了!和和 !嘿
摩特卡罗法!

#44

up

#1

关注

#2

你学过高等数学吗？
得到pi 的公式多的很，同让e也一样。
现在关键是要找到一个收敛速度快的函数来计算。

#3

http://sq.k12.com.cn:9000/tg/tgq163/sxyzl.htm

看看这篇文章，他的内容不是最新的了，但是值得参考。

#4

数学是很重要的!

#5

利用arctg45=1来计算

把arctgx展开成无穷级数，所求结果是这个无穷级数的4倍

#6

这问题有点意思，楼主你怎么想到问这样的问题的？

#7

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    Label1: TLabel;
    Label2: TLabel;
    Edit1: TEdit;
    Edit2: TEdit;
    procedure Edit1KeyPress(Sender: TObject; var Key: Char);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Edit1KeyPress(Sender: TObject; var Key: Char);
var
  k:integer;
  p,t,s,e:real;
  m:string;
begin
  if key=#13 then
    begin
      k:=strtoint(edit1.Text);
      if k>50 then
        begin
          showmessage('数值太大可能会使程序失去响应');
          edit1.Text:='18';
          k:=18;
        end;
      p:=0; s:=2;
      e:=exp(k*ln(0.1));
      repeat
        t:=s;
        p:=sqrt(2+p);
        s:=s*2/p;
      until abs(t-s)<e;
      m:=format('%20.18f',[s]);
      edit2.Text:=copy(m,0,k+2);
      edit1.SelStart:=0;
      edit1.SelLength:=length(edit1.Text);
      edit1.SetFocus;
    end;
end;

end.
说明一下：我是菜鸟，刚开始学delphi。别骂我啊！这是我书上的一道例题。

#8


我听说有人按 PI 的数位编出乐谱，再用计算机演奏出来很好听的

大家有兴趣可以做一个啊
做出来别忘了给我一份  ljlblue@163.com

#9

long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(;b-c;)
  f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
  for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

#10

你自己改成Pascal的
Okay？

#11

up

#12

这个程序我上学时做过习题,现在源码没有了,不过可以提两点建议.
1.直接按arctg45的级数展开是不可取的,因为它收敛太慢,正确的级数是两个反正切的差,具体我记不清了,但是你可以在<初中代数的教学参考书中找到>,我当时是在那里找的.
2.必须定义一种可以表示任意精度(也可说是任意长度整数)的数据结构,计算机本身是不提供的,需要你去设计.还要能过进行加\减\乘\除.实现时你可以研究一下汇编,看计算机是怎么实现的,然后将这种计算扩展到任意长.我当时做的就在Dephi中嵌的汇编.

完成这两件事,你的题就完成了
估计要几个小时的时间才行.

#13

我说的几小时是你编程的时间,不是计算机计算的时间
我算到2000位也用不了几秒种

自然对数e简单些,直接展开收敛也很快

#14


to： mingyeh(插死陈文登) 

   太历害了！佩服！
   短短两行代码可以算出那么多位PI！

   可是不知道你是依据什么公式计算的，什么原理，能达到多少位的精度
   望赐教！

#15


to: wan_ming(wan ming) 
*********************************************************

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TForm1 = class(TForm)
    GroupBox1: TGroupBox;
    Label1: TLabel;
    Label2: TLabel;
    Edit1: TEdit;
    Edit2: TEdit;
    procedure Edit1KeyPress(Sender: TObject; var Key: Char);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Edit1KeyPress(Sender: TObject; var Key: Char);
var
  k:integer;
  p,t,s,e:real;
  m:string;
begin
  if key=#13 then
    begin
      k:=strtoint(edit1.Text);
      if k>50 then
        begin
          showmessage('数值太大可能会使程序失去响应');
          edit1.Text:='18';
          k:=18;
        end;
      p:=0; s:=2;
      e:=exp(k*ln(0.1));
      repeat
        t:=s;
        p:=sqrt(2+p);
        s:=s*2/p;
      until abs(t-s)<e;
      m:=format('%20.18f',[s]);
      edit2.Text:=copy(m,0,k+2);
      edit1.SelStart:=0;
      edit1.SelLength:=length(edit1.Text);
      edit1.SetFocus;
    end;
end;

end.

这个算法的依据是:
   置初值  p=0,s=2
      重复计算：
          p=sqrt(2+p);    s=s*2/p;
   每重复一次 s 越接近于 PI

但这个算法最多只能达到高级语言本身数据类型(如 c 的 double )的精度，达不到我提出的要求

#16

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char * argv[])
{
  long * pi, * t, m, n, r, s;
  int t0[][3] = {48, 32, 20, 24, 8, 4}, k0[][3] = {1, 1, 0, 1, 1, 1};
  int n0[][3] = {18, 57, 239, 8, 57, 239}, d, i, j, k, p, q;

  d = (argc > 1) ? (((i = atoi(argv[1])) < 0) ? 0 : i) : 1000;
  q = (argc > 2) ? 1 : 0;
  printf("pi= %s%d * arctg(1/%d) %s %d * arctg(1/%d) %s %d * arctg(1/%d) [%s]\n",
    k0[q][0] ? "" : "-", t0[q][0], n0[q][0], k0[q][1] ? "+" : "-", t0[q][1],
    n0[q][1], k0[q][2] ? "+" : "-", t0[q][2], n0[q][2], q ? "Stomer" : "Gauss");
  if ((t = (long *)calloc((d += 5) + 1, sizeof(long))) == NULL) return 1;
  if ((pi = (long *)calloc(d + 1, sizeof(long))) == NULL) return 2;
  for (i = d; i >= 0; i--) pi[i] = 0;
  for (p = 0; p < 3; p++) {
    for (k=k0[q][p], n=n0[q][p], t[i=j=d]=t0[q][p], i--; i >= 0; i--) t[i] = 0;
    for (r = 0, i = j; i >= 0; i--) {
      r = (m = 10 * r + t[i]) % n;
      t[i] = m / n;
      k ? (pi[i] += t[i]) : (pi[i] -= t[i]);
    }
    while (j > 0 && t[j] == 0) j--;
    for (k = !k, s = 3, n *= n; j > 0; k = !k, s += 2) {
      for (r = 0, i = j; i >= 0; i--) {
        r = (m = 10 * r + t[i]) % n;
        t[i] = m / n;
      }
      while (j > 0 && t[j] == 0) j--;
      for (r = 0, i = j; i >= 0; i--) {
        r = (m = 10 * r + t[i]) % s;
        m /= s;
        k ? (pi[i] += m) : (pi[i] -= m);
      }
    }
  }
  for (n = i = 0; i <= d; pi[i++] = r) {
    n = (m = pi[i] + n) / 10;
    if ((r = m % 10) < 0) r += 10, n--;
  }
  printf("pi= %ld.", pi[d]);
  for (i = d - 1; i >= 5; i--)
    printf("%ld%s", pi[i], ((m = d - i + 5) % 65) ? ((m % 5) ? "" : " ") : "\n");
  printf("%sDIGITS: %d\n", (m % 65) ? "\n" : "", d - 5);
  return 0;
}

#17

!!!!!!!!!!!!!!!!
long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(;b-c;)
  f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
  for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

#18

上面的代码我也觉得奇怪，b没附初值这段这么来的，for(;b-c;)
  f[b++]=a/5;
这也能编译通过？？

#19

#define N 2800
long a=10000,b,c=N,d,e,f[N+1],g;
main()
{
while(b<c)
  f[b++]=a/5;

while(c>0){
    d=0; g=c*2; b=c;
    do{
       d+=f[b]*a; f[b]=d%--g;
       d/=g--;    b--; 
       if(b) d*=b; else break;
    }while(1);
    printf("%.4d",e+d/a);
    e=d%a;
    c-=14;
}
}

从 mingyeh(插死陈文登) 改进的（可读性好一些），VC 6 编译通过，没有问题，int 型默认初值为 0

#20

我试了，dcyu(Dd）和 mingyeh(插死陈文登) 的算法在 10000 位以内都是对的，10000以后就不知道了

#21

拉马努金的改进公式：
设f(n)=(13591409+545140134n)/640320^(3*n)*(-1)^n*(6n)!/(n!)^3/(3n)!
则1/pi=12/640320^(3/2)*(f(0)+f(1)+.......)
这个级数每增加一项能提高14位小数的精度

迭代方法：
y[0]=Sqrt(2)-1
z[n]=(1-y[n-1]^4)^(1/4)
y[n]=(1-z[n])/(1+z[n])
a[0]=6-4*Sqrt(2)
a[n]=(1+y[n])^4*a[n-1]-2^(2n+1)*y[n]*(1+y[n]+y[n]^2)
8次迭代能精确到小数点178814位

#22

计算e非常简单，而且收敛很快，用公式：
e=1+1/1!+1/2!+1/3!+1/4!+...

#23

我在书店里看到一本小册子专门讲如何算pi的

#24

long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(;b-c;)
  f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
  for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

这种算法的数据理论依据, 请参见这篇文章:
Stanley Rabinowitz, Stan wagon
A Spigot Algorithm for the Digits of pi
American Mathematical Monthly, 
vol. 102, Iss. 3 (Mar. 1995), pp195-203.

该期刊电子版网址为:
http://uk.jstor.org/journals/maa.html

==================================================
有关pi值计算的一些算法, 也可参见:
http://numbers.computation.free.fr/Constants/TinyPrograms/tinycodes.html
http://hades.ph.tn.tudelft.nl/Internal/PHServices/Documentation/MathWorld/math/math/p/p303.htm
http://www1.physik.tu-muenchen.de/~gammel/matpack/html/Mathematics/Pi.html

#25

呵呵，我有原代码，可是只能算几万的，算不了几亿的，算几亿靠我的
塞杨二566恐怕不行，！！！
如果只是几百万的话，还没问题，到几亿，而你又不想等的话，就必须靠
fft（算是FFT的变种吧！）具体贴起来可就多了。

#26

有二次收敛算法，就是如果这次精度是1000位，下一次计算可以达到2000位精度，就是比较复杂，要是几亿位，恐怕编程的时间要远远超过最后算出PI的时间.
sqrt()表示平方根运算
x0=sqrt(2)
pi0=2+sqrt(2)
y0=sqrt(sqrt(2))

x(i+1)=1/2*(sqrt(x(i)) + 1/sqrt(x(i)))
pi(i+1)=pi(i)*(x(i+1)+1)/(y(i)+1)
y(i+1)=(y(i)*sqrt(x(i+1)) + 1/sqrt(x(i+1))) / (y(i) + 1)

希望没有打错。
算法用到了平方根运算，乘除法，必须运用FFT算法。

#27

e就好算多了，直接求1/n!的和就可以了

#28

还有一种用概率来求的。用园和矩形，然后随机生成点。根据点落的在矩形内的多少。然后依据几何概型来求！

#29

恐怕编程的时间要远远超过最后算出PI的时间.
呵呵，这倒不至于，不过确实非常的麻烦！
而且很多都没办法，不得不用汇编写（用C的话效率就要低一倍以上，没办法。还的稍微利用一点硬件的特性！
哎算几亿，太麻烦了，
算个2万吧，又太容易。哎。

#30

去年大一的时候，因为老师悬赏，我写了一个能计算pi的程序。
用turbo pascal 7.0写的，当时整整花了一周的时间，什么事都没干。
程序能计算到8000位，1000位的时间大概在15分钟。当然由于当时急着交作业，也就没对算法进行优化。总的思路如下：
1.要定义4个高精度的计算，加、减、乘、除。
2.然后找一个收敛比较快的公式，（现在我忘了）
3.按照公式计算，只不过所有的运算都必须是高精度运算。
由于在dos下，内存的空间有限，因此数组最多只能到8000个元素，如果在windows中用delphi或者是c++builder可以计算得更多。
我写了200行的源代码，不过我觉得这挺有意思的，至少能够说明自己编程的水平。有兴趣的话，或者需要源程序，我们可以联系.
我的电邮是supercctv@sina.com

#31

http://www.jason314.com/
这个网站把PI介绍的很详细
我用其中的Bailey-Borwein-Plouffe算法的改进算法经过优化后计算PI小数点后100,000位仅用1秒钟。
还有更优的计算公式，不过我不大会实现。

#32

以100000000为例，我说的方法内存占用大于pi的内存占用的8倍，而pi占用
100000000/8*0.3010，所以总的内存占用为330MB左右，要是没有大的内存，就免了。还有，总共需要26步迭代，每次要是用到1分钟的时间，你的程序就很快了。记住，计算乘法使用的FFT算法，最大的运算位数只是几千位，要是计算一亿位，就会需要特殊的算法.

#33

vgrtgfg

#34

参考了手边的几本书, 比较的结果是在众多的展开式中, 以下式收敛较快

pi = 24 * arctan(1/8) + 8 * arctan(1/57) + 4 * arctan(1/239)

其中, 在计算时将 arctan 函数用级数展开.

不知大家有何看法?

#35

大家不用忙着贴代码，先听我慢慢说。下面的公式一般我们叫做马庭公式：
pi = 24 * arctan(1/8) + 8 * arctan(1/57) + 4 * arctan(1/239)
如果不是要求太高，用上面的公式就行了，原因就是因为编程比较简单。去年，我曾经花了几个月去研究这个问题。我找了不少资料，参照了这个网站的信息：http://www.jason314.com/ ，后来改进了站长的汇编程序，如果没记错，速度可以在几小时内算到一百万位。但这已经是马庭公式这类型的算法的极限。我国内下载的程序也没有超越这个速度的。这类的公式还有更快的，如下，但是我想速度也不会有什么突破。
Pi/4 = 22*arctan(1/28)+2*arctan(1/443)-5*arctan(1/1393)-10*arctan(1/11018)
Pi/4 = 44*arctan(1/57)+7*arctan(1/239)-12*arctan(1/682)+24*arctan(1/12943)

#36

现在回答楼主的问题，真正能算PI的程序我只见过两个。
第一个是流行测试 CPU 的 SuperPI 程序。这个程序是由日本人用 F2C 把   Fortran 的源程序改成 C 的，改得适合在PC机上使用，使用的就是 Gauss-Legendre algorithm （高斯－勒让德）算法。而那个Fortran程序就是由近年来一直刷新世界纪录的那两个日本人  Yasumasa KANADA (Associate Professor)  和  D.Takahashi  写的，其实是刷新自己的记录而已。他们使用的是Gauss-Legendre algorithm算法和Borwein's 4-th order convergent algorithm算法，也许Fortran更适合科学计算。具体的算法我就不贴出来了，如果大家有兴趣可以找一下，简单的说就是用迭代法。关键问题已经不是寻找算法，而是怎样在更好的在计算机里实现算法，如何存储结果等。

#37

Introduction of Kanada Lab.

Field of Research:
 High Performance Computing, Virtual Workstation System, Research on Research, Computer Algebra

Name of Laboratory:
 Laboratory for Computational tools

Location:
 University of Tokyo. Information science.

Name of staff:
 Yasumasa KANADA (Associate Professor)
 Akira YOSHIOKA  (Research Associate)

网站在东京大学里，要找还是可以的。

#38

连续的回复不能超过三次！

#39

第二个是PIFast，现在最新的版本是4.1，是一个法国人写的。日本人的程序是在大型机上运行的，而PIFast完完全全是在PC机运行的，主要是用了FFT，所以速度很快，快得你不相信，比PC版的SuperPI快很多。但是具体算法就不是很清楚，大家可以到他的网站看看。
http://numbers.computation.free.fr/Constants/constants.html

#40

网址我快不记得了，不过我下载的那个程序我特意也superPI比了一下，
算的少的话，看不出差距，越大差距越大，
那个程序比super pi 快多了！
不过他用了那个Apf233，里面一堆汇编文件，我好不容易才在VC6。0下编译通
过，运行。

#41

忘了告诉大家，目前PC机上流行的最快的圆周率计算程序是PiFast，PC机上的最高计算记录12,884,901,372位。使用Chudnovsky公式计算，速度最快。还可以使用Ramanujan公式，速度慢大约40%。PiFast可以利用磁盘缓存，突破物理内存的限制进行超高精度的计算，最高计算位数可达240亿位，并提供基于Fabrice Bellard公式的验算功能。它除了计算圆周率，还可以计算e和sqrt(2)。

#42

还可以用
多特魔德法！！！！

#43

说 错了!和和 !嘿
摩特卡罗法!

#44

up
标签：

相关文章