圆周率

kenmark 发表于 2008-3-9 15:42:48

一般的计算方法
1. Machin公式（马青公式）
pi=16arctan (1/5)-4arctan (1/239)
这个公式由英国天文学教授John Machin约翰·马青于1706年发现。他利用这个公式计算到了100位的圆周率。马青公式每计算一项可以得到1.4位的十进制精度。因为它的计算过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。
　　还有很多类似于马青公式的反正切公式。在所有这些公式中，马青公式似乎是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，马青公式就力不从心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用FFT（Fast Fourier Transform）算法。FFT可以将两个大数的乘除运算时间由O（n2）缩短为O（nlog（n））。
关于FFT算法的具体实现和源程序，请参考Xavier Gourdon的主页

2. Ramanujan公式（拉马努金公式）
\pi=\frac{9801}{2\sqrt2\sum_{n=0}^\infty\frac{(4n)!}{4^{4n}(n!)^4}\frac{(1103+26390n)}{99^{4n}}}
　　1914年，印度数学家Ramanujan拉马努金在他的论文里发表了一系列共14条圆周率的计算公式。这个公式每计算一项可以得到8位的十进制精度。1985年Gosper用这个公式计算到了圆周率的17,500,000位。
　　1989年，大卫·丘德诺夫斯基和格雷高里·丘德诺夫斯基兄弟将拉马努金公式改良，这个公式被称为丘德诺夫斯基公式，每计算一项可以得到15位的十进制精度。1994年丘德诺夫斯基兄弟利用这个公式计算到了4,044,000,000位。丘德诺夫斯基公式的另一个更方便于计算机编程的形式是：
\pi=\frac{426880\sqrt{10005}}{\sum_{n=0}^\infty\frac{(6n)!(545140134n+13591409)}{(n!)^3(3n)!(-640320)^{3n}}}

3. AGM（Arithmetic-Geometric Mean）算法
　　Gauss-Legendre高斯-勒让德公式：
初值：
a=x=1,b=\frac{1}{\sqrt2},c=\frac{1}{4}
重复计算：
y=a,a=\frac{a+b}{2},b=\sqrt{by},c=c-x(a-y)^2,x=2x
最后计算：
\pi=\frac{(a+b){}^2}{4c}

　　这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就够了。1999年9月，日本的高桥大介和金田康正用这个算法计算到了圆周率的206,158,430,000位，创出新的世界纪录。

4. Borwein波尔文四次迭代式：
　　这个公式由Jonathan Borwein乔纳森·波尔文和Peter Borwein彼得·波尔文于1985年发表，它四次收敛于圆周率。
a_0=6-4\sqrt2,y_0=\sqrt2-1
重复计算：
y_{n+1}=\frac{1-(1-y_n^4)^{\frac{1}{4}}}{1+(1-y_n^4)^{\frac{1}{4}}}
a_{n+1}=a_n(1+y_{n+1})^4-2^{2n+3}y_{n+1}(1+y_{n+1}+y_{n+1}^2)
最后计算：
\pi=\frac{1}{a_n}

5. Bailey-Borwein-Plouffe算法
\pi=\sum_{n=0}^\infty\frac{1}{16^n}(\frac{4}{8n+1}-\frac{2}{8n+4}-\frac{1}{8n+5}-\frac{1}{8n+6})
这个公式简称BBP公式，由David Bailey, Peter Borwein和Simon Plouffe于1995年共同发表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算前面的n-1位。这为圆周率的分布式计算提供了可行性。1997年，Fabrice Bellard找到了一个比BBP快40％的公式：
\pi=\frac{1}{64}\sum_{n=0}^\infty\frac{-1}{1024^n}(-\frac{32}{4n+1}-\frac{1}{4n+3}+\frac{256}{10n+1}-\frac{64}{10n+3}-\frac{4}{10n+5}-\frac{4}{10n+7}-\frac{4}{10n+9})

kenmark 发表于 2008-3-9 15:52:44

圆周率的计算历史

时间　　　纪录创造者　小数点后位数       所用方法
前2000 古埃及人 0
前1200　中国    0
前500 《圣经》　　 0 （周三径一）
前250　　阿基米德　　　 3
263 　　　刘徽　　　　　5                   古典割圆术
480          祖冲之             7
1429          Al-Kashi          14
1593          Romanus       15
1596          鲁道夫          20                   古典割圆术
1609          鲁道夫          35
1699          夏普                71                   夏普无穷级数
1706          马青             100                   马青公式
1719（法）德·拉尼       127（112位正确）夏普无穷级数
1794（奥地利）乔治·威加 140                   欧拉公式
1824（英）威廉·卢瑟福 208（152位正确）勒让德公式
1844          Strassnitzky & Dase          200
1847          Clausen          248
1853          Lehmann       261
1853       Rutherford       440
1874       威廉·山克斯       707(527位正确)
20世纪后
年月          纪录创造者                   所用机器       小数点后位数
1946             （英）弗格森　　　　　　    　　　　　　　　　620
1947 1       （英）弗格森　　　       　　　　　　　　　　　710
1947 9       Ferguson & Wrench 　　　　　       　　　　　　808
1949             Smith & Wrench 　　　　　　　　　       　　　1,120

1949             Reitwiesner et al　　　　　　ENIAC          2,037
1954          Nicholson & Jeenel　　　　　NORC　　　　3,092
1957                Felton          Pegasus          7,480
1958 1             Genuys          IBM704　       10,000
1958 5             Felton          Pegasus          10,021
1959                Guilloud          IBM 704          16,167
1961             Shanks & Wrench       IBM 7090       100,265
1966          Guilloud & Filliatre    IBM 7030       250,000
1967          Guilloud & Dichampt    CDC 6600       500,000
1973          Guilloud & Bouyer       CDC 7600    1,001,250
1981             Miyoshi & Kanada    FACOM M-200    2,000,036
1982                Guilloud                         2,000,050
1982                Tamura       MELCOM 900II    2,097,144
1982          Tamura & Kanada    HITACHI M-280H    4,194,288
1982          Tamura & Kanada    HITACHI M-280H    8,388,576
1983    Kanada, Yoshino & Tamura HITACHI M-280H    16,777,206
198510             Gosper       Symbolics 3670    17,526,200
1986 1             Bailey          CRAY-2       29,360,111
1986 9       Kanada & Tamura HITACHI S-810/20 33,554,414
198610       Kanada & Tamura HITACHI S-810/20 67,108,839
1987 1Kanada, Tamura & Kubo et al    NEC SX-2    134,217,700
1988 1       Kanada & Tamura HITACHI S-820/80 201,326,551
1989 5          Chudnovskys CRAY-2 & IBM-3090/VF480,000,000
1989 6          Chudnovskys       IBM 3090    525,229,270
1989 7       Kanada & Tamura HITACHI S-820/80 536,870,898
1989 8          Chudnovskys       IBM 3090 1,011,196,691
198911       Kanada & Tamura HITACHI S-820/801,073,741,799
1991 8          Chudnovskys                      2,260,000,000
1994 5          Chudnovskys                      4,044,000,000
1995 8    Takahashi & Kanada HITACHI S-3800/4804,294,967,286
199510    Takahashi & Kanada                   6,442,450,938
1997 7    Takahashi & Kanada                   51,539,600,000
1999 4    Takahashi & Kanada                   68,719,470,000
1999 9    Takahashi & Kanada HITACHI SR8000206,158,430,000
2002       Takahashi Team                   1,241,100,000,000

kenmark 发表于 2008-3-9 15:57:31

Mathtype->Tex转换效果太差了！！
2楼罗列的表格也无法正常～5～～～～
一般我都是用PiFast算着玩的呵呵，大家也可以下下来玩玩
此外，写一下Mathtype->Tex转换的经验
1.Preference->Translators->Tex -- LaTeX2.09 and later后面两个都不要
2.转换后复制转换后的文本到剪贴板，去除开头的\[和结束的\]
3.使用查找替换，去除所有的空格
4.查找"\sum\limits_"替换成"\sum_"
5.可以看了，但是还有很多冗余的括号，这些只能自己去了
终于全部改完了，手动修改公式真是累死了，总共修改了不少于40次，5～～～终于完成了！！！！

kenmark 发表于 2008-3-9 15:58:07

圆周率的最新计算纪录

1、新世界纪录
　　圆周率的最新计算纪录由日本人金田康正的队伍所创造。他们于2002年算出π值1,241,100,000,000 位小数，这一结果打破了他们于1999年9月18日创造的206,000,000,000位小数的世界纪录。

2、个人计算圆周率的世界纪录
　　在一个现场解说验证活动中，一名59岁日本老人Akira Haraguchi将圆周率π算到了小数点后的83431位，这名孜孜不倦的59岁老人向观众讲解了长达13个小时，最终获得认同。这一纪录已经被收入了Guinness世界大全中。据报道，此前的纪录是由一名日本学生于1995年计算出的，当时的精度是小数点后的42000位。

kenmark 发表于 2008-3-9 15:59:05

G++编译器中的运算程序（转自百度）

微机WindowsXP中Dev-cpp中的运算程序（30000位）（C++）
#include < cstdlib.h >
#include < iostream.h >
#include < fstream.h >
#define N 30010
using namespace std;
void mult (int *a,int b,int *s)
{
for (int i=N,c=0;i>=0;i--)
{
      int y=(*(a+i))*b+c;
      c=y/10;
      *(s+i)=y%10;
}
}
void divi (int *a,int b,int *s)
{
for (int i=0,c=0;i<=N;i++)
{
      int y=(*(a+i))+c*10;
      c=y%b;
      *(s+i)=y/b;
}
}
void incr(int *a,int *b,int *s)
{
for (int i=N,c=0;i>=0;i--)
{
      int y=(*(a+i))+(*(b+i))+c;
      c=y/10;
      *(s+i)=y%10;
}
}
bool eqs(int *a,int *b)
{
int i=0;
while (((*(a+i))==(*(b+i)))&&(i<=N)) i++;
return i>N;
}
int main(int argc, char *argv[])
{
int lpi,lls,lsl,lp;

int *pi=lpi,*ls=lls,*sl=lsl,*p=lp;
for (int i=0;i<=N;i++)*(pi+i)=*(ls+i)=*(sl+i)=*(p+i)=0;
memset(pi,0,sizeof(pi));
memset(ls,0,sizeof(ls));
memset(sl,0,sizeof(sl));
memset(p,0,sizeof(p));
*pi=*ls=*sl=1;
for (int i=1;true;i++)
{
   mult(ls,i,sl);
   divi(sl,2*i+1,ls);
   incr(pi,ls,p);
   if (eqs(pi,p)) break;
   int *t;
   t=p;
   p=pi;
   pi=t;
   if (i%50==0) cout << i << "";
}
cout << endl;
mult(p,2,pi);
ofstream fout("pi.txt");
fout << *pi << ".";
for (int i=1;i<=N;i++)
{
   fout << *(pi+i);
   if (i%10==0) fout << " ";
   if (i%80==0) fout << endl;
}
return EXIT_SUCCESS;
}
注：①运行时会有数据弹出，那是无关紧要的，只为了加快了感觉速度；
   ②最后的txt文本里有30010位，其中最后10位可能是错的。
   ③程序中的<，>是大写的请改成小写。

gxqcn 发表于 2008-3-9 19:15:50

原帖由 kenmark 于 2008-3-9 15:57 发表 https://bbs.emath.ac.cn/static/image/common/back.gif
Mathtype->Tex转换效果太差了！！
2楼罗列的表格也无法正常～5～～～～
一般我都是用PiFast算着玩的呵呵，大家也可以下下来玩玩
此外，写一下Mathtype->Tex转换的经验
1.Preference->Translators->Tex -- LaTeX2 ...

MathType->TeX 转换效果确实不理想，尤其是多了许多不必要的括号。

表格输入是支持的，发帖时切换输入模式到“所见即所得模式”，点击"添加表格"按钮。但排版好也不算容易。

nyy 发表于 2023-9-8 15:08:34

现在的世界纪录，都是用Chudnovsky公式

lihpb00 发表于 2023-9-8 22:15:31

我只会最简单的反三角函数泰勒级数展开式

页: [1]

数学研发论坛's Archiver

圆周率

圆周率的计算历史

圆周率的最新计算纪录

G++编译器中的运算程序（转自百度）