大きい数＝工夫or対数!

大きい数を扱う…例えば階乗を計算しろ、という問題や組み合わせ\(_nC_r\)を計算しろ、という問題に時々あたります。

まず、典型的な例として組み合わせ
\(\displaystyle _nC_r=\frac{n!}{(n-r)!r!}\)
を計算しましょう。
この計算はそのまま定義のまま計算しようとすると階乗が出てくるので、途中の値は物凄くでかい数字になるのですが、最後の答えは小さい値になるはずです。それなのに計算できないのはおかしくない!?工夫次第で何とかなるものです。
下の表は３つの計算方法についてまとめたものです。

定義のまま計算

例えばfortranでintegerの宣言により定義した場合、階乗が正しく計算できる範囲(オーバーフローしない範囲)は、\(12!=479\ 001\ 600\)
まで計算できます。
変数の型がいくつまで計算できるかは、組み込み関数hugeを使うことで確認できます。

これで確認すると、整数型は
\(\pm 2\ 147\ 483\ 647\)
までの値なら代入することができるため、
\(13!=6\ 227\ 020\ 800\)は計算できないことがわかります。
一応関数はこんな感じ。

( ※余談ですが、8バイト型整数で宣言、すなわちinteger(8)で宣言すると
\(9\ 223\ 372\ 036\ 854\ 775\ 807\)
まで計算できます(19桁、\(20!=2\ 432\ 902\ 008\ 176\ 640\ 000\)まで)。)

なので定義通りの計算はn=12までが限界です。実用的じゃありません。

工夫して計算

\(\displaystyle _nC_r=\frac{n!}{(n-r)!r!}\)を工夫します。

\(
\begin{align}
\displaystyle _nC_r&=\frac{n!}{(n-r)!r!} \\
&=\frac{n(n-1)(n-2)\cdots(n-r+2)(n-r+1)}{r(r-1)(r-2)\cdots2\cdot 1} \\
&=\frac{n}{r}\cdot \frac{n-1}{r-1}\cdot \frac{n-2}{r-2}\cdots \frac{n-r+2}{2}\cdot \frac{n-r+1}{1} \\
&=\prod_{i=1}^r\frac{n-r+i}{i}
\end{align}
\)
として計算します。こうすれば階乗みたいに大きな数を計算する必要がなくなります。
関数を書けばこうなります。

で\(_nC_r\)を求めることができます。
上のプログラムでは、最大\(_{29}C_{14}=77\ 558\ 760\)まで求めることができます。
これは絶対に割り切れることを前提にしているためにおこるもので、

の文のところで最大値が決まってしまうためです。
定式化すれば、
式\(_nC_r\)に対し、\(_{n-1}C_{r-1}\cdot n\)が扱える整数値を超えないこと
と表現できます。

\(_{29}C_{14}=77\ 558\ 760\)が計算できるのは、
\(_{29-1}C_{14-1}*29=1\ 085\ 822\ 640\)
であるので計算可能です。次の\(_{30}C_{15}\)を計算するためには
\(_{30-1}C_{15-1}*30=2\ 326\ 762\ 800\)が
integerの最大の整数値\(\pm 2\ 147\ 483\ 647\)
を超えてしまうのでこれ以上計算はできません。
※もしもinteger(8)を用いると\(_{61}C_{30}\)まで計算できます。

対数を利用して計算

もう一つ対数を利用した方法を紹介します。
対数を利用すると、
\(\displaystyle _nC_r=\frac{n!}{(n-r)!r!}\)は、\(S= _nC_r\)とおくと
\(
\begin{align}
\displaystyle \ln(S)&=\ln\left[\frac{n!}{(n-r)!r!}\right] \\
&=\ln(n!)-\ln\{(n-r)!\}-\ln(r!) \\
&=\sum_{i=n-r+1}^n\ln(i)-\sum_{i=2}^{r}\ln(i)
\end{align}
\)
と書けます。よって、
\(\displaystyle S=\exp\left[\sum_{i=n-r+1}^n\ln(i)-\sum_{i=2}^{r}\ln(i)\right]\)
です。プログラムを書けば

となります。この場合は\(_{n}C_{r}\)が整数型で扱える範囲を越えなければいいという条件だけになるので、
integerでは\(_{33}C_{16}\)までokです。

integer(8)では\(_{66}C_{33}\)までできます…が、上のプログラムそのままではできません。
おそらく、組み込み関数nint(実数型を整数型に四捨五入して代入する)がinteger(8)に対応していない(かも)ということと、倍精度の有効桁数がだんだんと足らなくなっていくので厳しいです。後者による問題によって、確認した範囲では\(_{48}C_{24}\)位が限界です。

Tweet