有哪位兄弟知道super pi的算法？

yangjio4849

如题。我想编一个fortran程序，测试CPU浮点运算速度，这还是我们的作业呢。

nuclearweapon

看看级数

yangjio4849

是反三角函数的级数吗？我就奇怪如何将一个fortran里的数存到那么多位？有效数字为何会那么多？

yangjio4849

我找到反正弦函数的展开了，但是展开得极慢，50000000000项得到的pai=  3.14155814，还差得远呢。有收敛快的吗？
程序如下，非常简单的一个：
      program arcsin
      double precision pai,xishu
      print *,'tell me the number of terms'
      read (*,*) n
      pai=2.0+1.0/3.0
      xishu=1.0
      do 10 i=2,n
       xishu=xishu*(2*i-1)/2.0/i
       pai=pai+xishu*1.0/(2*i+1)
10    continue
      print *,'pai=',pai
      end

还有，我就不明白windows下的superpai是如何保存pai的数值的，好像fortran里面没有什么可以达到几百万位的有效数字啊。

soloforce

用傅立叶级数展开的办法可行；但本人采用1/4倍率收敛的级数仍然比较缓慢。
可以肯定superpi不仅仅是简单的采用级数展开计算的

我这里有个级数展开式：

1. 
   
2. pi=3 + (1*1)/(2*3*4) * [ 3 + (3*3)/(4*5*4) *[ 3 + (5*5)/(6*7*4) * [ 3 + (7*7)/(8*9*4) * [ ... ] ] ] ]
   

复制代码

soloforce

Post by yangjio4849
我找到反正弦函数的展开了，但是展开得极慢，50000000000项得到的pai=  3.14155814，还差得远呢。有收敛快的吗？
程序如下，非常简单的一个：
      program arcsin
      double precision pai,xishu
      print *,'tell me the number of terms'
      read (*,*) n
      pai=2.0+1.0/3.0
      xishu=1.0
      do 10 i=2,n
       xishu=xishu*(2*i-1)/2.0/i
       pai=pai+xishu*1.0/(2*i+1)
10    continue
      print *,'pai=',pai
      end

还有，我就不明白windows下的superpai是如何保存pai的数值的，好像fortran里面没有什么可以达到几百万位的有效数字啊。

你这简单累加绝对不行的，没有任何一种语言可以达到你要求的精度。
可以考虑分段计算，即把小数点后面每n位放在一指定数组段中。

soloforce

贴两个算pi的小程序，采用普遍的傅立叶级数展开方法，不同之处是采用的级数。速度是绝对不如superpi的，仅供参考。使用时请保留著作权信息。

程序一：

1. 
   
2. /*************************************************************************************************************
   
3. 
   
4.         pi值计算程序 pi_v1
   
5. 
   
6.         作者：welans @ linuxsir.cn ,2005.5. Copyright reserved
   
7. 
   
8.         测试平台：Linux2.6.11-nitro2 + gcc3.4.3 @ AthlonXP1600+, 1G DDR266
   
9. 
   
10. 
    
11.         简介：
    
12. 
    
13.                 计算圆周率pi到小数点后约5.5万位精度
    
14. 
    
15.                 根据公式 pi=3 + (1*1)/(2*3*4) * [ 3 + (3*3)/(4*5*4) *[ 3 + (5*5)/(6*7*4) * [ 3 + (7*7)/(8*9*4) * [ ... ] ] ] ]       
    
16. 
    
17.                 每STEP计算8位，每计算8位需占用14个unsigned long long int型数组单元
    
18. 
    
19.        
    
20.         本地速度评估：5000位小于1秒,50000位小于50秒( 比起super pi 差远了 :-< )
    
21. 
    
22. 
    
23. 
    
24.         输入：需计算的圆周率位数，要求小于55000，不保证55000以后的精度。
    
25. 
    
26.         输出：计算得到的圆周率结果保存到文件'pi.dat'中。
    
27. 
    
28.        
    
29. 
    
30.         编译：gcc pi.c -o pi -O3
    
31. ************************************************************************************************************/
    
32. 
    
33. #include <stdio.h>
    
34. 
    
35. 
    
36. 
    
37. #define MIN_LEN 8
    
38. 
    
39. #define WIDTH        8
    
40. 
    
41. #define STEP 14
    
42. 
    
43. 
    
44. 
    
45. int ACC_CORRECT(int acc)
    
46. 
    
47. {
    
48. 
    
49.         if(acc<MIN_LEN || acc% WIDTH)
    
50. 
    
51.         {
    
52. 
    
53.                 printf("Accuracy length should be a multiple of %d and no less than %d!\n",WIDTH,MIN_LEN);
    
54. 
    
55.                 return 0;
    
56. 
    
57.         }
    
58. 
    
59.         else return 1;
    
60. 
    
61. }
    
62. 
    
63. 
    
64. 
    
65. int main()
    
66. 
    
67. {
    
68. 
    
69.          int ACCURACY_LEN;           // default accuracy length.
    
70. 
    
71.          int MAX;                           // size of the array.
    
72. 
    
73.          unsigned long long int *a;                           // pointer of the array.
    
74. 
    
75.          unsigned long long int i,m,t,k,d,g=30000000;
    
76. 
    
77.        
    
78. 
    
79.         FILE *fp=NULL;
    
80.         if((fp=fopen("pi.dat","wb"))==NULL)
    
81. 
    
82.         {
    
83. 
    
84.                 printf("Can not create 'pi.dat'!\n");
    
85. 
    
86.                 return -1;
    
87. 
    
88.         }
    
89. 
    
90. 
    
91.         do
    
92. 
    
93.         {
    
94. 
    
95.                 printf("Input accuracy length:");
    
96. 
    
97.                 scanf("%d",&ACCURACY_LEN);
    
98. 
    
99.         }
    
100. 
     
101.         while( !ACC_CORRECT(ACCURACY_LEN));
     
102. 
     
103. 
     
104. 
     
105.         MAX=(ACCURACY_LEN*STEP)/WIDTH;
     
106. 
     
107.         a= malloc(MAX*sizeof(unsigned long long int));
     
108. 
     
109. 
     
110. 
     
111.         // main procedure for computing pi.
     
112. 
     
113.         for(i=0;i<MAX;i++) a[i]=30000000;
     
114. 
     
115.         for(m=MAX;m>0;m-=STEP)
     
116. 
     
117.         {
     
118. 
     
119.                 for (t=0,i=m-1;i>0;i--)
     
120. 
     
121.                 {
     
122. 
     
123.                         t+=a[i]*100000000;
     
124. 
     
125.                         k=i*(i<<4)+(i<<3);
     
126. 
     
127.                         a[i]=t%k;
     
128. 
     
129.                         d=(i<<1)-1;
     
130. 
     
131.                         t=t/k*d*d;
     
132. 
     
133.                 }
     
134. 
     
135.                 fprintf(fp,"%.8d ",g+t/100000000);
     
136. 
     
137.                 g=t%100000000;
     
138. 
     
139.         }       
     
140. 
     
141.         //end of computing pi
     
142. 
     
143.         fprintf(fp,"\n");
     
144.         printf("\nResults saved to 'pi.dat'!\n");
     
145. 
     
146.         fclose(fp);
     
147. 
     
148.         free(a);
     
149. 
     
150.         return 0;
     
151. 
     
152. }
     
153. 
     
154. 
     

复制代码

程序二：

1. 
   
2. /*************************************************************************************************************
   
3. 
   
4.         pi值计算程序 pi_v2
   
5. 
   
6.         作者：welans @ linuxsir.cn ,2005.5. Copyright reserved
   
7. 
   
8.         测试平台：Linux2.6.11-nitro2 + gcc3.4.3 @ AthlonXP1600+, 1G DDR266
   
9. 
   
10. 
    
11.         简介：
    
12. 
    
13.                 计算圆周率pi到小数点后约30万位精度
    
14. 
    
15.                 根据公式 pi=2+1/3*(2+2/5*(2+3/7*(2+4/9*(2+5/11  .....))))))
    
16. 
    
17.         本地速度评估：5000位近1秒,50000位小于70秒,10万位小于5分钟。效率略低于piv1。
    
18. 
    
19. 
    
20. 
    
21.         输入：需计算的圆周率位数，理论精度约30万位。
    
22. 
    
23.         输出：计算得到的圆周率结果保存到文件'pi.dat'中。
    
24. 
    
25.        
    
26. 
    
27.         编译：gcc pi.c -o pi -O3
    
28. ************************************************************************************************************/
    
29. 
    
30. #include <stdio.h>
    
31. 
    
32. 
    
33. 
    
34. #define MIN_LEN        8
    
35. 
    
36. #define WIDTH                8
    
37. 
    
38. #define STEP 27
    
39. 
    
40. 
    
41. 
    
42. int ACC_CORRECT(int acc)
    
43. 
    
44. {
    
45. 
    
46.         if(acc<MIN_LEN || acc% WIDTH)
    
47. 
    
48.         {
    
49. 
    
50.                 printf("Accuracy length should be a multiple of %d and no less than %d!\n",WIDTH,MIN_LEN);
    
51. 
    
52.                 return 0;
    
53. 
    
54.         }
    
55. 
    
56.         else return 1;
    
57. 
    
58. }
    
59. 
    
60. 
    
61. 
    
62. int main()
    
63. 
    
64. {
    
65. 
    
66.         unsigned int ACCURACY_LEN;
    
67. 
    
68.         unsigned int MAX ;
    
69. 
    
70.         unsigned long long int i,m,t,d,g=20000000;
    
71. 
    
72.         unsigned long long int *a;
    
73. 
    
74. 
    
75.         FILE *fp=NULL;
    
76.         if((fp=fopen("pi.dat","wb"))==NULL)
    
77. 
    
78.         {
    
79. 
    
80.                 printf("Can not create 'pi.dat'!\n");
    
81. 
    
82.                 return -1;
    
83. 
    
84.         }
    
85. 
    
86. 
    
87.         do
    
88. 
    
89.         {
    
90. 
    
91.                 printf("Input accuracy length:");
    
92. 
    
93.                 scanf("%d",&ACCURACY_LEN);
    
94. 
    
95.         }
    
96.         while(!ACC_CORRECT(ACCURACY_LEN));
    
97. 
    
98. 
    
99.         MAX=ACCURACY_LEN*STEP/WIDTH;
    
100. 
     
101.         a= malloc(MAX*sizeof(unsigned long long int));
     
102. 
     
103. 
     
104.         for(i=0;i<MAX;i++)a[i]=20000000;
     
105. 
     
106.         for(m=MAX;m>0;m-=STEP)
     
107. 
     
108.         {
     
109. 
     
110.                 for (t=0,i=m-1;i>0;i--)
     
111. 
     
112.                 {
     
113. 
     
114.                         t+=a[i]*100000000;
     
115. 
     
116.                         d=(i<<1)+1;
     
117. 
     
118.                         a[i]=t%d;
     
119. 
     
120.                         t=t/d*i;
     
121. 
     
122.                 }
     
123. 
     
124.                 fprintf(fp,"%.8d ",g+t/100000000);
     
125. 
     
126.                 g=t%100000000;
     
127. 
     
128.         }
     
129. 
     
130.        
     
131. 
     
132.         printf("\nResult saved to file 'pi.dat'!\n");
     
133. 
     
134.         fclose(fp);
     
135. 
     
136.         free(a);
     
137. 
     
138.         return 0;
     
139. 
     
140. }
     

复制代码

yangjio4849

Post by welans
贴两个算pi的小程序，采用普遍的傅立叶级数展开方法，不同之处是采用的级数。速度是绝对不如superpi的，仅供参考。使用时请保留著作权信息。

程序一：

1. 
   
2. /*************************************************************************************************************
   
3. 
   
4.         pi值计算程序 pi_v1
   
5. 
   
6.         作者：welans @ linuxsir.cn ,2005.5. Copyright reserved
   
7. 
   
8.         测试平台：Linux2.6.11-nitro2 + gcc3.4.3 @ AthlonXP1600+, 1G DDR266
   
9. 
   
10. 
    
11.         简介：
    
12. 
    
13.                 计算圆周率pi到小数点后约5.5万位精度
    
14. 
    
15.                 根据公式 pi=3 + (1*1)/(2*3*4) * [ 3 + (3*3)/(4*5*4) *[ 3 + (5*5)/(6*7*4) * [ 3 + (7*7)/(8*9*4) * [ ... ] ] ] ]       
    
16. 
    
17.                 每STEP计算8位，每计算8位需占用14个unsigned long long int型数组单元
    
18. 
    
19.        
    
20.         本地速度评估：5000位小于1秒,50000位小于50秒( 比起super pi 差远了 :-< )
    
21. 
    
22. 
    
23. 
    
24.         输入：需计算的圆周率位数，要求小于55000，不保证55000以后的精度。
    
25. 
    
26.         输出：计算得到的圆周率结果保存到文件'pi.dat'中。
    
27. 
    
28.        
    
29. 
    
30.         编译：gcc pi.c -o pi -O3
    
31. ************************************************************************************************************/
    
32. 
    
33. #include <stdio.h>
    
34. 
    
35. 
    
36. 
    
37. #define MIN_LEN 8
    
38. 
    
39. #define WIDTH        8
    
40. 
    
41. #define STEP 14
    
42. 
    
43. 
    
44. 
    
45. int ACC_CORRECT(int acc)
    
46. 
    
47. {
    
48. 
    
49.         if(acc<MIN_LEN || acc% WIDTH)
    
50. 
    
51.         {
    
52. 
    
53.                 printf("Accuracy length should be a multiple of %d and no less than %d!\n",WIDTH,MIN_LEN);
    
54. 
    
55.                 return 0;
    
56. 
    
57.         }
    
58. 
    
59.         else return 1;
    
60. 
    
61. }
    
62. 
    
63. 
    
64. 
    
65. int main()
    
66. 
    
67. {
    
68. 
    
69.          int ACCURACY_LEN;           // default accuracy length.
    
70. 
    
71.          int MAX;                           // size of the array.
    
72. 
    
73.          unsigned long long int *a;                           // pointer of the array.
    
74. 
    
75.          unsigned long long int i,m,t,k,d,g=30000000;
    
76. 
    
77.        
    
78. 
    
79.         FILE *fp=NULL;
    
80.         if((fp=fopen("pi.dat","wb"))==NULL)
    
81. 
    
82.         {
    
83. 
    
84.                 printf("Can not create 'pi.dat'!\n");
    
85. 
    
86.                 return -1;
    
87. 
    
88.         }
    
89. 
    
90. 
    
91.         do
    
92. 
    
93.         {
    
94. 
    
95.                 printf("Input accuracy length:");
    
96. 
    
97.                 scanf("%d",&ACCURACY_LEN);
    
98. 
    
99.         }
    
100. 
     
101.         while( !ACC_CORRECT(ACCURACY_LEN));
     
102. 
     
103. 
     
104. 
     
105.         MAX=(ACCURACY_LEN*STEP)/WIDTH;
     
106. 
     
107.         a= malloc(MAX*sizeof(unsigned long long int));
     
108. 
     
109. 
     
110. 
     
111.         // main procedure for computing pi.
     
112. 
     
113.         for(i=0;i<MAX;i++) a[i]=30000000;
     
114. 
     
115.         for(m=MAX;m>0;m-=STEP)
     
116. 
     
117.         {
     
118. 
     
119.                 for (t=0,i=m-1;i>0;i--)
     
120. 
     
121.                 {
     
122. 
     
123.                         t+=a[i]*100000000;
     
124. 
     
125.                         k=i*(i<<4)+(i<<3);
     
126. 
     
127.                         a[i]=t%k;
     
128. 
     
129.                         d=(i<<1)-1;
     
130. 
     
131.                         t=t/k*d*d;
     
132. 
     
133.                 }
     
134. 
     
135.                 fprintf(fp,"%.8d ",g+t/100000000);
     
136. 
     
137.                 g=t%100000000;
     
138. 
     
139.         }       
     
140. 
     
141.         //end of computing pi
     
142. 
     
143.         fprintf(fp,"\n");
     
144.         printf("\nResults saved to 'pi.dat'!\n");
     
145. 
     
146.         fclose(fp);
     
147. 
     
148.         free(a);
     
149. 
     
150.         return 0;
     
151. 
     
152. }
     
153. 
     
154. 
     

复制代码

程序二：

1. 
   
2. /*************************************************************************************************************
   
3. 
   
4.         pi值计算程序 pi_v2
   
5. 
   
6.         作者：welans @ linuxsir.cn ,2005.5. Copyright reserved
   
7. 
   
8.         测试平台：Linux2.6.11-nitro2 + gcc3.4.3 @ AthlonXP1600+, 1G DDR266
   
9. 
   
10. 
    
11.         简介：
    
12. 
    
13.                 计算圆周率pi到小数点后约30万位精度
    
14. 
    
15.                 根据公式 pi=2+1/3*(2+2/5*(2+3/7*(2+4/9*(2+5/11  .....))))))
    
16. 
    
17.         本地速度评估：5000位近1秒,50000位小于70秒,10万位小于5分钟。效率略低于piv1。
    
18. 
    
19. 
    
20. 
    
21.         输入：需计算的圆周率位数，理论精度约30万位。
    
22. 
    
23.         输出：计算得到的圆周率结果保存到文件'pi.dat'中。
    
24. 
    
25.        
    
26. 
    
27.         编译：gcc pi.c -o pi -O3
    
28. ************************************************************************************************************/
    
29. 
    
30. #include <stdio.h>
    
31. 
    
32. 
    
33. 
    
34. #define MIN_LEN        8
    
35. 
    
36. #define WIDTH                8
    
37. 
    
38. #define STEP 27
    
39. 
    
40. 
    
41. 
    
42. int ACC_CORRECT(int acc)
    
43. 
    
44. {
    
45. 
    
46.         if(acc<MIN_LEN || acc% WIDTH)
    
47. 
    
48.         {
    
49. 
    
50.                 printf("Accuracy length should be a multiple of %d and no less than %d!\n",WIDTH,MIN_LEN);
    
51. 
    
52.                 return 0;
    
53. 
    
54.         }
    
55. 
    
56.         else return 1;
    
57. 
    
58. }
    
59. 
    
60. 
    
61. 
    
62. int main()
    
63. 
    
64. {
    
65. 
    
66.         unsigned int ACCURACY_LEN;
    
67. 
    
68.         unsigned int MAX ;
    
69. 
    
70.         unsigned long long int i,m,t,d,g=20000000;
    
71. 
    
72.         unsigned long long int *a;
    
73. 
    
74. 
    
75.         FILE *fp=NULL;
    
76.         if((fp=fopen("pi.dat","wb"))==NULL)
    
77. 
    
78.         {
    
79. 
    
80.                 printf("Can not create 'pi.dat'!\n");
    
81. 
    
82.                 return -1;
    
83. 
    
84.         }
    
85. 
    
86. 
    
87.         do
    
88. 
    
89.         {
    
90. 
    
91.                 printf("Input accuracy length:");
    
92. 
    
93.                 scanf("%d",&ACCURACY_LEN);
    
94. 
    
95.         }
    
96.         while(!ACC_CORRECT(ACCURACY_LEN));
    
97. 
    
98. 
    
99.         MAX=ACCURACY_LEN*STEP/WIDTH;
    
100. 
     
101.         a= malloc(MAX*sizeof(unsigned long long int));
     
102. 
     
103. 
     
104.         for(i=0;i<MAX;i++)a[i]=20000000;
     
105. 
     
106.         for(m=MAX;m>0;m-=STEP)
     
107. 
     
108.         {
     
109. 
     
110.                 for (t=0,i=m-1;i>0;i--)
     
111. 
     
112.                 {
     
113. 
     
114.                         t+=a[i]*100000000;
     
115. 
     
116.                         d=(i<<1)+1;
     
117. 
     
118.                         a[i]=t%d;
     
119. 
     
120.                         t=t/d*i;
     
121. 
     
122.                 }
     
123. 
     
124.                 fprintf(fp,"%.8d ",g+t/100000000);
     
125. 
     
126.                 g=t%100000000;
     
127. 
     
128.         }
     
129. 
     
130.        
     
131. 
     
132.         printf("\nResult saved to file 'pi.dat'!\n");
     
133. 
     
134.         fclose(fp);
     
135. 
     
136.         free(a);
     
137. 
     
138.         return 0;
     
139. 
     
140. }
     

复制代码

不错啊，大哥是自己编写的吗？好厉害。
我测试了一下，v1，50000位大概40.973s，v2同样位数大概54.671s，Athlon3000+和512M内存。什么时候把这程序改成fortran的看看。

yangjio4849

这么看来流行的super pi软件是用一种更加高效的级数来做的了？真想知道啊。

soloforce

我不认同。
收敛越快的级数，往往伴随着复杂度的增加，使单步计算的时间复杂度也上升了。
比如：我的v1级数收敛速度是v2的两倍，不过计算消耗时间仅仅略少于v2。

《C数值算法》（[A]Numerical Recipes in C，数值计算权威著作）一书中最后一章非典型数值算法中也给出了一个算Pi的程序，程序本身和采用的递推式很复杂，但其程序在我机器上运行效率远低于v1和v2。(限于许可证协议，我不能贴出该程序)

我觉得可以考虑并行计算和分布计算。

附件是一个网络上的图片，不知道其正确性，供参考