BBP 式を使用して円周率を計算するコード

BBP 式は、円周率の N 番目の結果を直接取得できると主張していますが、取得速度はどの位置でも同じであると考えていましたが、bellard のホームページ (http://bellard.org/pi/) からダウンロードしてください。改良した C プログラムをテストしたところ、N が大きくなるほど計算時間が長くなりました。PHP コードに変換した後、元の C プログラムの計算も 1K 桁目を計算するのに 1 分以上かかりました。遅い。

このプログラムは一度に 9 つの数値のみを出力します (囧rz、より多くの数値を出力するように変更する方法がわかりません)

/**
 * 円周率計算 (BBP)
 * @author Moyo
 * @url http://moyo.uuland.org/code/php-pi-calc/
 * @version 1.0
 * @日付 2013.01.12
*/

class pi
{
 public static function calc($__N__)
 {
 $n = (int)$__N__;
 $av = $ a = $vmax = $N = $num = $den = $k = $kq = $kq2 = $t = $v = $s = $i = 0;
 $sum = 0.0;
 $N = (int) (($n + 20) * log(10) / log(2));
 $sum = 0;
 for ($a = 3; $a  {
 $vmax = (int)(log(2 * $N) / log($a));
 $av = 1;
 for ($i = 0; $i  {
 $av = ($av * $a);
 }
 $s = 0;
 $num = 1;
 $den = 1;
 $v = 0;
 $ kq = 1;
 $kq2 = 1;
 for ($k = 1; $k  {
 $t = $k;
 if ($kq >= $a )
 {
 do
 {
 $t = (int)($t / $a);
 $v --;
 }
 while (($t % $a) == 0);
 $kq = 0 ;
 }
 $kq ++;
 $num = self::mul_mod($num, $t, $av);
 $t = (2 * $k -1);
 if ($kq2 >= $ a)
 {
 if ($kq2 == $a)
 {
 do
 {
 $t = (int)($t / $a);
 $v ++;
 }
 while (($t % $a) == 0);
 }
 $kq2 -= $a;
 }
 $den = self::mul_mod($den, $t, $av);
 $kq2 += 2;
 if ($ v > 0)
 {
 $t = self::inv_mod($den, $av);
 $t = self::mul_mod($t, $num, $av);
 $t = self::mul_mod ($t, $k, $av);
 for ($i = $v; $i  {
 $t = self::mul_mod($t, $a, $ av);
 }
 $s += $t;
 if ($s >= $av)
 {
 $s -= $av;
 }
 }
 }
 $t = self::pow_mod(10 , ($n - 1), $av);
 $s = self::mul_mod($s, $t, $av);
 $sum = (double)fmod((double)$sum + (double)$ s / (double)$av, 1.0);
 }
 return array(
 'n' => $n,
 'v' => sprintf('%09d', (int)($sum * 1e9) )
 );
 }
 プライベート静的関数 next_prime($n)
 {
 do
 {
 $n ++;
 }
 while (!self::is_prime($n));
 return $n;
 }
プライベート静的関数 is_prime($n)
 {
 $r = $i = 0;
 if (($n % 2) == 0)
 {
 return 0;
 }
 $r = (int)(sqrt ($n));
 for ($i = 3; $i  {
 if (($n % $i) == 0)
 {
 return 0;
 }
 }
 return 1;
 }
 プライベート静的関数 mul_mod($a, $b, $m)
 {
 return fmod((double)$a * (double)$b, $m);
 }
プライベート静的関数 inv_mod($x, $y)
 {
 $q = $u = $v = $a = $c = $t = 0;
 $u = $x;
 $v = $y;
 $ c = 1;
 $a = 0;
 do
 {
 $q = (int)($v / $u);
 $t = $c;
 $c = $a - $q * $c;
 $a = $t;
 $t = $u;
 $u = $v - $q * $u;
 $v = $t;
 }
 while ($u != 0);
 $a = $ a % $y;
 if ($a < 0)
 {
 $a = $y + $a;
 }
 return $a;
 }
 プライベート静的関数 pow_mod($a, $b, $m)
 {
 $r = $aa = 0;
 $r = 1;
 $aa = $a;
 while (1)
 {
 if ($b & 1)
 {
 $r = self::mul_mod( $r, $aa, $m);
 }
 $b = $b >> 1;
 if ($b == 0)
 {
 Break;
 }
 $aa = self::mul_mod($ aa, $aa, $m);
 }
 return $r;
 }
}

?>