情報基礎A 「Cプログラミング」(ステップ8・「大きな数」の計算)

再帰的な関数の応用として、「大きな数」を計算する方法を考える。

1.大きな数が登場する計算

コンピュータが扱うことのできる数値の範囲は当然ながら限られている。 その範囲は、コンピュータやOSの種類によっても異なるが、C言語でint型(符号付きの整数) で表現可能な範囲は以下のプログラムを実行して確かめることができる:

ヘッダファイル limits.h には最大・最小値を含むいくつかの値が定義されている。 なお、実数型の最大・最小値は float.h に定義されている。

#include <stdio.h>
#include <limits.h>
main() {
  printf("int型の最小値 %d\n",INT_MIN) ;
  printf("int型の最大値 %d\n",INT_MAX) ;
}

よく普及している32ビットマシンの場合、表現できる整数の範囲は-2147483648〜2147483647 である。 そして、これはかなり窮屈と言わざるを得ない。 というのは、例えば、「組み合わせ」を階乗を使って $$ {n \choose k} = \frac{n!}{(n-k)! k!} $$ で計算しようとすると、たかだか$n=13$の場合ですら$13! = 6227020800 > 2147483647$となって、扱える数値の範囲を超えてしまう。

大きな数が登場する別の例として、自然数のべき乗の剰余の計算がある。 $a,p,n$を正の整数とすると $$ a^b \mod n $$ ($a$の$b$乗を$n$で割った余り)の計算は、素数の判定(フェルマーテスト)や、 RSA暗号など、情報に関係した分野では、しばしば登場する。 この場合も$b$の値が大きくなると$a^b$はすぐに表現可能な整数の範囲を超えてしまうだろう。

以下では、この大きい数の剰余の計算方法を考察してみよう。

2.小さな数を使った計算に分解する

$x$を$n$で割った商$q$と剰余$r$に分解したとき($x = q n + r$)、それを$s$乗した数の剰余について $$ (q n + r)^s = r^s \mod n $$ が成り立つことはすぐにわかる (上式の左辺を展開すると、$r^s$の項以外は全て$n$のべき乗がかかる、 すなわち、$(q n + r)^s = q^s n^s + s q^{s-1} n^{s-1} r + \cdots + s q n r^{s-1} + r^s$ であることに注意)。

また、2つの数$x = q_1 n + r_1, y= q_2 n + r_2$について、その積の剰余について $$ x y = r_1 r_2 \mod n $$ が成り立つことも自明であろう。

これらを使うと、以下の手順で、剰余の値を保ったまま、$a^b$を変形することができる:

  1. $a^k > n$となるような$k$を探す
  2. 上で求めた$k$を使って $$ a^b = (a^k)^{s} \times a^t $$ と「分解」する。ここで$b = s k + t$、すなわち$s$は$b$を$k$で割った商、$t$はその剰余である。
  3. ここで、$a^k$を$n$で割った余りを$r$と置く、すなわち、$a^k = q n + r$とすると、上の議論から、$(a^k)^s = r^s \mod n$
  4. こうして、$a^b$の剰余を、より小さな2つの数($r^s, a^t$)の剰余の計算に「分解」することができた: $$ a^b = \left( r^s \mod n \right) \left( a^t \mod n \right) \mod n $$

上記の「分解」のプロセスを忠実になぞって、$a^b \mod n$を与えるCの関数int mpower(int a, int b, int n) を書くと、以下のようなコードになるだろう:

int mpower(int a, int b, int n)
{
  int k,a1,a2,b1,b2 ;
  if (b==0) return 1 ;
  if (b==1) return a%n ;
  a1=a ;
  k=1 ;
  while (a1<n && k<b) {
    a1=a1*a ;
    k=k+1 ;
  }
  a1=a1%n ;
  b1=b/k ;
  a2=a ;
  b2=b%k ;
  return (mpower(a1,b1,n)*mpower(a2,b2,n))%n ;
}

計算が進んで、$a^1$のパターンに至ったら、さらなる分解は行わず、素直に剰余を計算する(return a%n)。 また、(計算したいのは$n>1$の場合なので)明らかに$a^0 \mod n = 1$である(return 1)。

icon-pc 練習: 動作確認

上で設計した関数が正しく動作するかどうか、プログラムを完成させて、動作確認しなさい。

上の例では、大きな数を$a^b = a^{ks +r}$と分解したが、(効率は悪くても構わないので)他の分解方法を用いたアルゴリズムを考案し、それを実装してみなさい。

icon-hint ヒント

通常の電卓などでは検算ができないので、例えば、Online Big Number Calculatorの計算欄に

(2^100)%1234

のように入力して、結果を較べるとよい(上の例は$2^{100} \mod 1234$)。

icon-pc 練習: 組み合わせの計算

int型の整数計算を使って、できるだけ大きな$n, k$に対しても$n \choose k$を計算することができるよう、 Cの関数int combination(int n, int k)を設計しなさい。

icon-hint ヒント

階乗を使って計算する場合、途中で大きな数が出現しないよう、適宜、約分しながら、計算を進める必要があるだろう。 最大公約数を見つけるには、ユークリッドの互除法が使える。

別のアプローチとして、二項定理パスカルの三角形も手がかりを与えてくれるかもしれない。