Pythonプログラミング（ステップ８・「大きな数」の計算）

１．大きな数が登場する計算

コンピュータが扱うことのできる数値の範囲は当然ながら限られている。 その範囲は、コンピュータやOSの種類によっても異なるが、 よく普及している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$はすぐに表現可能な整数の範囲を超えてしまうだろう。

それではあまりに窮屈ということで、Python 2のlong型やPython 3のint型では、扱える桁数に事実上の上限は無い。 しかしながら、コンピュータのハードウェアには「大きな数」を直接扱う機能は無いので、 大きな数の計算には時間がかかったりメモリーを消費したりする。

以下では、この大きい数の剰余を、なるだけ効率的に計算する方法を考察してみよう。

２．小さな数を使った計算に分解する

$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$ であることに注意）。

また、２つの数$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$の剰余を、より小さな２つの数（$r^s, a^t$）の剰余の計算に「分解」することができた： $$ a^b = \left( r^s \mod n \right) \left( a^t \mod n \right) \mod n $$

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

def mpower(a, b, n):
    if b==0:
        return 1
    if b==1:
        return a%n
    a1=a
    k=1
    while a1<n and 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)。

　練習： 動作確認

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

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

　ヒント

通常の電卓などでは検算ができないので、Pythonをインタラクティブモードで起動して

$ python
Python x.x.x ....
>>> (2**10000000)%1234 [Enter]
1024

のように結果を得ることができる(この場合、$2^{10000000} \mod 1234$）。 こういう計算が簡単にできるのなら、プログラミングは必要無さそうに思えるかもしれないが、べき指数 10000000 をさらに大きな数にすると、 直接的な計算では、結果が得られなくなってしまうだろう（計算時間とメモリ不足）。

また、ネット上にも、 例えば、Online Big Number Calculator といったサービスがある。こちらの場合は、計算欄に

(2^10000000)%1234

のように入力して、結果を較べることができる。

　練習： 組み合わせの計算

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

　ヒント

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

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