PHPで階乗計算

まずは普通に再帰関数で。

~~~php:test1.php
<?php
function factorial($n) {
    if($n > 0) return $n * factorial($n - 1);
    return 1;
}

for($n = 1; $n <= 180; ++$n) {
    echo "$n! = ", factorial($n), "\n";
}
~~~

結果

~~~txt
1! = 1
2! = 2
3! = 6
4! = 24
5! = 120
6! = 720
7! = 5040
8! = 40320
9! = 362880
10! = 3628800
11! = 39916800
12! = 479001600
13! = 6227020800
14! = 87178291200
15! = 1307674368000
16! = 20922789888000
17! = 355687428096000
18! = 6402373705728000
19! = 121645100408832000
20! = 2432902008176640000
21! = 5.1090942171709E+19
22! = 1.1240007277776E+21
23! = 2.5852016738885E+22
24! = 6.2044840173324E+23
        :
       中略
        :
168! = 2.5260757449732E+302
169! = 4.2690680090047E+304
170! = 7.257415615308E+306
171! = INF
172! = INF
173! = INF
        :
~~~

当然ですが、21!から指数表記に、171!以降はINFに。
ちょっと物足りないので、文字列と配列で1桁ずつ10進数演算するよう変更。

~~~php:test2.php
<?php
function factorial($n) {
    if($n > 0) return intMul($n, factorial($n - 1));
    return 1;
}

function intMul($a, $b) {
    $a = "$a";
    $b = "$b";
    $r = [];
    $al = strlen($a);
    $bl = strlen($b);

    for($i = 0; $i < $al; ++$i) {
        $_a = $al - 1 - $i;
        for($j = 0; $j < $bl; ++$j) {
            $_b = $bl - 1 - $j;
            $k = $i + $j;
            if(!isset($r[$k])) $r[$k] = 0;
            $r[$k] += $a[$_a] * $b[$_b];
            if($r[$k] > 9) {
                $k1 = $k + 1;
                if(!isset($r[$k1])) $r[$k1] = 0;
                $r[$k1] += (int)($r[$k] / 10);
                $r[$k] %= 10;
            }
        }
    }
    return strrev(implode('', $r));
}

for($n = 1; $n <= 180; ++$n) {
    echo "$n! = ", factorial($n), "\n";
}
~~~

結果

~~~txt
1! = 1
2! = 2
3! = 6
4! = 24
5! = 120
6! = 720
7! = 5040
8! = 40320
9! = 362880
10! = 3628800
11! = 39916800
12! = 479001600
13! = 6227020800
14! = 87178291200
15! = 1307674368000
16! = 20922789888000
17! = 355687428096000
18! = 6402373705728000
19! = 121645100408832000
20! = 2432902008176640000
21! = 51090942171709440000
22! = 1124000727777607680000
23! = 25852016738884976640000
24! = 620448401733239439360000
        :
       中略
        :
180! = 20089606249913429965 …中略… 986039603200000000000000000000000000000000000000000000
~~~

しっかり表示されるものの、1桁ずつ計算しているので結構遅いです。
私の普段使いの数年前のノートPCだと、5000!程度の計算でも数十秒かかります。

~~~php
$n = 5000;
$start = hrtime(true);
$result = factorial($n);
$end = hrtime(true);
echo "$n! = $result\n", ($end - $start) / 1e9, "\n";

// 5000! = 422857792660554 …中略… 5509378334720000  …中略… 0000000000
// 47.1204944
~~~

ある程度の桁数をまとめて計算できるようにしてみます。

~~~php
echo 999999999 * 999999999; // 999999998000000001
echo 9999999999 * 9999999999; // 9.999999998E+19
~~~

9桁同士までなら大丈夫そうなのでとりあえず9桁で。

~~~php:test3.php
<?php
function factorial($n) {
    if($n > 0) return intMul($n, factorial($n - 1));
    return 1;
}

define('DIGIT', 9);
define('NUMBER', 10 ** DIGIT);
define('F', '%0'. DIGIT. 'd');

function intMul($a, $b) {
    $a = "$a";
    $b = "$b";

    $a_ = $a;
    $ra = [];
    $asn = strlen($a) / DIGIT;
    for($i = 0; $i < $asn; ++$i) {
        array_push($ra, (int)substr($a_, -DIGIT));
        $a_ = substr($a_, 0, -DIGIT);
    }

    $b_ = $b;
    $rb = [];
    $bsn = strlen($b) / DIGIT;
    for($i = 0; $i < $bsn; ++$i) {
        array_push($rb, (int)substr($b_, -DIGIT));
        $b_ = substr($b_, 0, -DIGIT);
    }

    $c = [];
    for($i = 0; $i < count($ra); ++$i) {
        for($j = 0; $j < count($rb); ++$j) {
            $p = $i + $j;
            if(!isset($c[$p])) $c[$p] = 0;
            $c[$p] += $ra[$i] * $rb[$j];
            if($c[$p] >= NUMBER) {
                if(!isset($c[$p + 1])) $c[$p + 1] = 0;
                $c[$p + 1] += (int)($c[$p] / NUMBER);
                $c[$p] %= NUMBER;
            }
        }
    }

    for($i = 0; $i < count($c) - 1; ++$i) $c[$i] = sprintf(F, $c[$i]);
    return implode('', array_reverse($c));
}

$n = 5000;
$start = hrtime(true);
$result = factorial($n);
$end = hrtime(true);
echo "$n! = $result\n", ($end - $start) / 1e9, "\n";

// 5000! = 422857792660554 …中略… 5509378334720000  …中略… 0000000000
// 4.2777483
~~~

10倍くらい早くなりました。

自分で組むのはひとまず置いて、[GMP](https://www.php.net/manual/ja/book.gmp.php)に階乗関数[gmp_fact()](https://www.php.net/manual/ja/function.gmp-fact.php)があるのでこれを使ってみます。

~~~php:test4.php
<?php
$n = 5000;
$start = hrtime(true);
$result = gmp_fact($n);
$end = hrtime(true);
echo "$n! = ", gmp_strval($result), "\n", ($end - $start) / 1e9, "\n";

// 5000! = 422857792660554 …中略… 5509378334720000  …中略… 0000000000
// 0.0003974
~~~

更に10000倍くらい早くなりました。

**1000000!の場合**

~~~txt
// 1000000! = 82639316883312400623 …中略… 25617650584125440000  …中略… 0000000000
// 0.3422455
~~~

**10000000!の場合**

~~~txt
// 10000000! = 12024234005159034561 …中略… 03073025741946880000  …中略… 0000000000
// 6.7120242
~~~

階乗を使いたいけれど自分で書きたいわけではないのなら、`gmp_fact()`を使うのが良さそうですね。

**追記**
5000!に約4.2秒かかっていた自前計算版(test3.php)も、乗算処理を階乗計算用としてループ内に入れて配列/文字列変換処理コストを節約することで、同じアルゴリズムのままでも3倍くらい早くなりました。
それでも`gmp_fact()`には到底及びませんが。

~~~php:test5.php
<?php
function factorial($f) {
    if($f < 0 || $f !== (int)$f) return false;

    $d = 9;
    $dn = 10 ** $d;
    $df = "%0{$d}d";

    $c = [1];
    for($n = 1; $n <= $f; ++$n) {
        $ra = $c;
        $c = [0];
        $b_ = "$n";
        $rb = [];
        $bsn = strlen($n) / $d;
        for($i = 0; $i < $bsn; ++$i) {
            array_push($rb, (int)substr($b_, -$d));
            $b_ = substr($b_, 0, -$d);
        }
        for($i = 0; $i < count($ra); ++$i) {
            for($j = 0; $j < count($rb); ++$j) {
                $p = $i + $j;
                if(!isset($c[$p])) $c[$p] = 0;
                $c[$p] += $ra[$i] * $rb[$j];
                if($c[$p] >= $dn) {
                    $p1 = $p + 1;
                    if(!isset($c[$p1])) $c[$p1] = 0;
                    $c[$p1] += (int)($c[$p] / $dn);
                    $c[$p] %= $dn;
                }
            }
        }
    }
    for($i = 0; $i < count($c) - 1; ++$i) $c[$i] = sprintf($df, $c[$i]);
    return implode('', array_reverse($c));
}

$n = 5000;
$start = hrtime(true);
$result = factorial($n);
$end = hrtime(true);
echo "$n! = $result\n", ($end - $start) / 1e9, "\n";

// 5000! = 422857792660554 …中略… 5509378334720000 …中略… 00000000
// 1.3479136
~~~