找回密码
 欢迎注册
查看: 9056|回复: 7

[欣赏] 圆周率

[复制链接]
发表于 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})$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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
1985  10                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
1986  10          Kanada & Tamura     HITACHI S-810/20     67,108,839
1987   1  Kanada, 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/VF  480,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
1989  11          Kanada & Tamura     HITACHI S-820/80  1,073,741,799
1991   8             Chudnovskys                        2,260,000,000
1994   5             Chudnovskys                        4,044,000,000
1995   8       Takahashi & Kanada   HITACHI S-3800/480  4,294,967,286
1995  10       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 SR8000  206,158,430,000
2002           Takahashi Team                       1,241,100,000,000
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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~~~终于完成了!!!!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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位。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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[N+1],lls[N+1],lsl[N+1],lp[N+1];
   
    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位可能是错的。
       ③程序中的<,>是大写的请改成小写。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2008-3-9 19:15:50 | 显示全部楼层
原帖由 kenmark 于 2008-3-9 15:57 发表
Mathtype->Tex转换效果太差了!!
2楼罗列的表格也无法正常~5~~~~
一般我都是用PiFast算着玩的呵呵,大家也可以下下来玩玩
此外,写一下Mathtype->Tex转换的经验
1.Preference->Translators->Tex -- LaTeX2 ...


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

表格输入是支持的,发帖时切换输入模式到“所见即所得模式”,点击"添加表格"按钮。但排版好也不算容易。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-8 15:08:34 | 显示全部楼层
现在的世界纪录,都是用Chudnovsky公式
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2023-9-8 22:15:31 | 显示全部楼层
我只会最简单的反三角函数泰勒级数展开式
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-4-20 05:56 , Processed in 0.087851 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表