使用BBP公式计算Pi的代码

BBP公式号称可以直接获取Pi的第N位结果，本来我以为这种获取的速度在任意一个位置都是相同的呢，但是从bellard的主页（http://bellard.org/pi/）下载一个改进的C程序测试后发现N越大，计算时间越长，我转换成php代码后，计算第1K位时会超过1分钟的时间，原C程序计算的时候也是很慢

这个程序一次只会输出9个数字（囧rz，不知道怎样改能多输出几位）

2. /**
3. * 圆周率计算(BBP)
4. * @author Moyo
5. * @url http://moyo.uuland.org/code/php-pi-calc/
6. * @version 1.0
7. * @date 2013.01.12
8. */
10. class pi
11. {
12. public static function calc($__N__)
13. {
14. $n = (int)$__N__;
15. $av = $a = $vmax = $N = $num = $den = $k = $kq = $kq2 = $t = $v = $s = $i = 0;
16. $sum = 0.0;
17. $N = (int)(($n + 20) * log(10) / log(2));
18. $sum = 0;
19. for ($a = 3; $a {
20. $vmax = (int)(log(2 * $N) / log($a));
21. $av = 1;
22. for ($i = 0; $i {
23. $av = ($av * $a);
24. }
25. $s = 0;
26. $num = 1;
27. $den = 1;
28. $v = 0;
29. $kq = 1;
30. $kq2 = 1;
31. for ($k = 1; $k {
32. $t = $k;
33. if ($kq >= $a)
34. {
35. do
36. {
37. $t = (int)($t / $a);
38. $v --;
39. }
40. while (($t % $a) == 0);
41. $kq = 0;
42. }
43. $kq ++;
44. $num = self::mul_mod($num, $t, $av);
45. $t = (2 * $k -1);
46. if ($kq2 >= $a)
47. {
48. if ($kq2 == $a)
49. {
50. do
51. {
52. $t = (int)($t / $a);
53. $v ++;
54. }
55. while (($t % $a) == 0);
56. }
57. $kq2 -= $a;
58. }
59. $den = self::mul_mod($den, $t, $av);
60. $kq2 += 2;
61. if ($v > 0)
62. {
63. $t = self::inv_mod($den, $av);
64. $t = self::mul_mod($t, $num, $av);
65. $t = self::mul_mod($t, $k, $av);
66. for ($i = $v; $i {
67. $t = self::mul_mod($t, $a, $av);
68. }
69. $s += $t;
70. if ($s >= $av)
71. {
72. $s -= $av;
73. }
74. }
75. }
76. $t = self::pow_mod(10, ($n - 1), $av);
77. $s = self::mul_mod($s, $t, $av);
78. $sum = (double)fmod((double)$sum + (double)$s / (double)$av, 1.0);
79. }
80. return array(
81. 'n' => $n,
82. 'v' => sprintf('%09d', (int)($sum * 1e9))
83. );
84. }
85. private static function next_prime($n)
86. {
87. do
88. {
89. $n ++;
90. }
91. while (!self::is_prime($n));
92. return $n;
93. }
94. private static function is_prime($n)
95. {
96. $r = $i = 0;
97. if (($n % 2) == 0)
98. {
99. return 0;
100. }
101. $r = (int)(sqrt($n));
102. for ($i = 3; $i {
103. if (($n % $i) == 0)
104. {
105. return 0;
106. }
107. }
108. return 1;
109. }
110. private static function mul_mod($a, $b, $m)
111. {
112. return fmod((double)$a * (double)$b, $m);
113. }
114. private static function inv_mod($x, $y)
115. {
116. $q = $u = $v = $a = $c = $t = 0;
117. $u = $x;
118. $v = $y;
119. $c = 1;
120. $a = 0;
121. do
122. {
123. $q = (int)($v / $u);
124. $t = $c;
125. $c = $a - $q * $c;
126. $a = $t;
127. $t = $u;
128. $u = $v - $q * $u;
129. $v = $t;
130. }
131. while ($u != 0);
132. $a = $a % $y;
133. if ($a {
134. $a = $y + $a;
135. }
136. return $a;
137. }
138. private static function pow_mod($a, $b, $m)
139. {
140. $r = $aa = 0;
141. $r = 1;
142. $aa = $a;
143. while (1)
144. {
145. if ($b & 1)
146. {
147. $r = self::mul_mod($r, $aa, $m);
148. }
149. $b = $b >> 1;
150. if ($b == 0)
151. {
152. break;
153. }
154. $aa = self::mul_mod($aa, $aa, $m);
155. }
156. return $r;
157. }
158. }
160. ?>

复制代码