php 使用BBP公式计算Pi的代码

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

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

<?php

/**
 * 圆周率计算(BBP)
 */

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 <= (2 * $N); $a = self::next_prime($a))
        {
            $vmax = (int)(log(2 * $N) / log($a));
            $av = 1;
            for ($i = 0; $i < $vmax; $i ++)
            {
                $av = ($av * $a);
            }
            $s = 0;
            $num = 1;
            $den = 1;
            $v = 0;
            $kq = 1;
            $kq2 = 1;
            for ($k = 1; $k <= $N; $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 < $vmax; $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))
        );
    }
    private static function next_prime($n)
    {
        do
        {
            $n ++;
        }
        while (!self::is_prime($n));
        return $n;
    }
    private static function is_prime($n)
    {
        $r = $i = 0;
        if (($n % 2) == 0)
        {
            return 0;
        }
        $r = (int)(sqrt($n));
        for ($i = 3; $i <= $r; $i += 2)
        {
            if (($n % $i) == 0)
            {
                return 0;
            }
        }
        return 1;
    }
    private static function mul_mod($a, $b, $m)
    {
        return fmod((double)$a * (double)$b, $m);
    }
    private static function 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;
    }
    private static function 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;
    }
}

?>
//该片段来自于http://outofmemory.cn