大きな数を扱う

大きい数=工夫or対数!

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

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

integer(4バイト整数型±2 147 483 647)を用いた計算
定義のまま計算 工夫して計算 対数を利用して計算
最大のn,r \(_{12}C_6\) \(_{29}C_{14}\) \(_{33}C_{16}\)
理由 階乗の計算によるオーバーフロー “割り切れる”ための制約 (限界値)

integer(8バイト整数型±9 223 372 036 854 775 807)を用いた計算
計算方法 定義のまま計算 工夫して計算 対数を利用して計算
最大のn,r \(_{20}C_{10}\) \(_{61}C_{30}\) \(_{48}C_{24}\)
理由 階乗の計算によるオーバーフロー “割り切れる”ための制約 倍精度型変数の有効桁数が足りない
定義のまま計算

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

program main
  implicit none
  integer::n

  write(6,*)"huge",huge(n)

  stop
end program main

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

function nCr_fact(n,r)
  !sikinote
  implicit none
  integer::nCr_fact,fact,n,r
  external::fact
 
  nCr_fact=fact(n)/(fact(n-r)*fact(r))
 
  return
end function nCr_fact
function fact(n)
  implicit none
  integer::fact,i,n
 
  if(n.le.-1)write(6,*)"####warning#### parameter of factorial has negative value"
  fact=1
  do i=2,n
     fact=fact*i
  enddo

  return
end function fact

( ※余談ですが、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}
\)
として計算します。こうすれば階乗みたいに大きな数を計算する必要がなくなります。
関数を書けばこうなります。

function nCr(n,r)
  !sikinote
  implicit none
  integer::n,r,i,r0,nCr

  r0=n-r
  if(r.le.r0)r0=r

  nCr=1
  do i=1,r0
     nCr=nCr*(n-r0+i)
     nCr=nCr/i
  enddo
 
  return
end function nCr

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

nCr=nCr*(n-r0+i)

の文のところで最大値が決まってしまうためです。
定式化すれば、
\(_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]\)
です。プログラムを書けば

function nCr_log(n,r)
  !sikinote
  implicit none
  integer::n,r,i,r0,nCr2
  double precision::tmp
 
  r0=n-r
  if(r.le.r0)r0=r

  tmp=0.d0
  do i=n-r0+1,n
     tmp=tmp+log(dble(i))
  enddo
  do i=2,r0
     tmp=tmp-log(dble(i))
  enddo

  nCr2=nint(exp(tmp))  

  return
end function nCr_log

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

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


「大きな数を扱う」への2件のフィードバック

  1. 今日は、宇井と申します。
    nCrについてのコメントです。
    配列を用意して
    1,0,0,0,0,・・・
    と初期化します。
    全体を1つずらして加算するをくりかえせば
    1,1,0,0,0,・・・
    1,2,1,0,0,・・・
    1,3,3,1,0,・・・
    と計算されて行きます。
    逆演算も可能でnを自在に増減出来ますよ。
    nが大きい時に多少時間がかかりますが、nとrをずらしながら使用する場合に便利です。

    1. こんにちは。
      パスカルの三角形によるnCrの求め方ですね。
      そういった方法もありましたか…気づきませんでした。
      コメントありがとうございます。
      複数の、さほど大きくないnCrを求める際には確かに有利ですね。

      デメリットがあるとすると
      ・プログラムが煩雑になるという点
      でしょうか。高速に計算できると思いますが
      nCrのプログラムが数式上で直観的ではなくなる点がつらいですね。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です