まずは普通に再帰関数で。
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";
}
結果
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進数演算するよう変更。
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";
}
結果
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!程度の計算でも数十秒かかります。
$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
ある程度の桁数をまとめて計算できるようにしてみます。
echo 999999999 * 999999999; // 999999998000000001
echo 9999999999 * 9999999999; // 9.999999998E+19
9桁同士までなら大丈夫そうなのでとりあえず9桁で。
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に階乗関数gmp_fact()があるのでこれを使ってみます。
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!の場合
// 1000000! = 82639316883312400623 …中略… 25617650584125440000 …中略… 0000000000
// 0.3422455
10000000!の場合
// 10000000! = 12024234005159034561 …中略… 03073025741946880000 …中略… 0000000000
// 6.7120242
階乗を使いたいけれど自分で書きたいわけではないのなら、gmp_fact()を使うのが良さそうですね。
追記
5000!に約4.2秒かかっていた自前計算版(test3.php)も、乗算処理を階乗計算用としてループ内に入れて配列/文字列変換処理コストを節約することで、同じアルゴリズムのままでも3倍くらい早くなりました。
それでもgmp_fact()には到底及びませんが。
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