冪乗法
(Power Algorithm)

   整数を法とした計算の中で,

181006655297358 (mod 181006655297359)

のように,非常に大きい冪を取る計算をする必要があることがしばしばあります。この例はフェルマーの素数判定法で出てきました。この計算結果は法とする数以下の数であることが分かっていますから,工夫すれば Tiny Basic 等の普通の言語で計算することが出来ます。しかも比較的高速に計算する方法として冪乗法が良く知られています。

 ここではその冪乗法(Power Algorithm)を説明します。
 
冪乗法
 

 一般に an (mod m) の計算をする最も素朴な方法は, an を計算して,その結果を m で割った時の余りを求めるものです。しかしこの方法は n が小さい時のみ有効な方法です。 実際,an はnが大きくなるとそれに従って急速に大きな数になります。ですからすぐに計算精度をオーバーして正確な値を計算することは出来ません。

 次に考えられるのは,a1, a2, と順次計算をし,その度にm で割った時の余りを求め,それに a を掛けて行くという計算です。例えば,

4 (mod 9)

を計算するのに,

1≡ 5 (mod 9)
≡ 51×5=5×5=25≡ 7(mod 9)
≡ 5×5=7×5=35≡ 8(mod 9)
≡ 5×5=8×5=40≡ 4(mod 9)


と計算する方法です。これですと,計算に表れる数は常に,法とする整数の2乗以下ですから,法が余り大きくなければ正確な計算が可能です。

 しかしこの方法は重大な欠点を持っています。それは計算ステップが指数に表れる数と同じ程度の回数になることです。この方法で,

181006655297358 (mod 181006655297359)

を計算するには,181006655297358 回もの計算ステップを実行しなくてはなりません。 この回数は膨大で,高速な計算機でも実行は現実的ではありません。


 冪乗法はこの方法の計算回数を減らす工夫をしたものです。例えば,

21 (mod 9)

を計算するとします。上の方法だと,21(20?)ステップの計算です。冪乗法はこの計算に補助項

≡7(mod 9),5≡4(mod 9),5≡7(mod 9), 516≡4(mod 9)

を計算します。この計算は各々その前の項の2乗ですから,計算ステップは4回と考えられます。これを使って,

20 (mod 9)≡516×5×5≡4×4×5≡8(mod 9)

を求めます。これらの計算ステップは2回ですから,合計6回ステップで計算が出来たことになります。

 この計算は指数20を2進法で展開し,その各桁が1か0かによって対応する補助的項を掛けています。

冪乗法のプログラム

 これをプログラムにすると以下のようになります。ここで説明のため行番号を付けています。


01 Function PowerMod(a,n,m)
02    pw = 1
03    While n >= 1
04       If (n mod 2) = 1 then
05          pw = a * pw mod m
06       End If
07       a = a^2 mod m
08       n = n \ 2
09    Wend
10    PowerMod = pw
11 End Function

 この関数は PowerMod(a,n,m) として,a^n mod m を返します。このプログラムでは上でのaの冪である補助項を 再び a に入れています。

2行でのpw は最後に値として返す power です。

3行から9行が計算の主要部です。4行では指数nが2進展開で桁が1かどうかの判定をしています。1の場合,そのときの補助項a を掛けます。

7行では,補助項 a を2乗し次のステップの補助項にします。

8行では,指数nを2で割った商をとり,2進展開で1つ上の桁の計算に移ります。


 このプログラムを使って Tiny Basic で実際にPowerMod(123,4567890123,456789)を計算させると,

1234567890123  ≡ 108531 (mod 456789)

が得られます。

乗法への応用

 上の冪乗法による高い指数の法計算は,簡単で強力なものですが,少し不満が残ります。それは上の方法では,PowerMod(a,n,m) の計算途中で m^2 程度の数が表れることです。このことは,上のプログラムは m が扱える整数の最大数の平方根以下でないと,正しい結果を返さないことを意味します。

 例えば, Tiny  Basic での整数計算は64ビット計算で行われますから,桁数は19桁,従って,m は9桁以下の数でないと正確さは保障されません。もっと大きなmについても計算できないでしょうか。

 実は冪乗法の考え方を加法的に応用し,乗法に付いての方法とすると乗法が大きいmについて可能になります。


 例えば,13×11(mod 15)の計算を考えます。普通に算すると,13×11=143が出てきますから,15^2程度の精度が必要です。しかし,これを

11=2^3+2+1と2進展開し,

補助項

×13≡13 (mod 15),
×13≡11(mod 15),
×13≡2×11≡7(mod 15),
×13≡2×7≡14(mod 15)

を計算して,

13×11≡13×2+13×2+13×1≡8(mod 15)

を求めると,途中の計算では全て30以下,即ち m の2倍以下の数しか表れません。

 ですから,この方法を適用すれば,最大整数の平方根でなく最大整数の1/2程度の m まで計算でkることになります。

 この方法はよく考えると,冪乗法での積を加法に置き換えたものに他ならないことが分かります。ですからそのプログラムも冪乗法のものと驚くほど良く似ています。

加法的 Power Algorithm を使った乗法

 整数 m を法とした乗法の計算


01 Function MulMod(a,n,m)
02    mu = 0
03    While n >= 1
04       If (n mod 2) = 1 then
05          mu = (a + mu) mod m
06       End If
07       a = a * 2 mod m
08       n = n \ 2
09    Wend
10    MulMod = mu
11 End Function

Function PowerMod との類似性に注目してください。

このプログラムを上のPowerMod に組み込むと,m が最大整数の1/2程度の数までの計算が可能になります。

 上のプログラムでは a が負だとうまく動きません。そこで以下では,負の数が入力された場合,mod  m での正の数に直してから計算をしています。


Function PowerMod(a,n,m)
   a = a mod m
   if a < 0 then a = a + m
   pw = 1
   While n >= 1
      If (n mod 2) = 1 then
         pw = MulMod(a,pw,m)
      End If
      a = MulMod(a,a,m)
      n = n \ 2
   Wend
   PowerMod = pw
End Function

Function MulMod(a,n,m)
   a = a mod m
   if a < 0 then a = a + m
   mu = 0
   While n >= 1
      If (n mod 2) = 1 then
         mu = (a + mu) mod m
      End If
      a = a * 2 mod m
      n = n \ 2
   Wend
   MulMod = mu
End Function
まとめ

 冪乗法は上のように簡単なものですが,m を法とする計算では強力です。この方法の計算回数はおよそa,n,mの最大数の2を底とする対数の2乗に比例する程度のステップで す。従って大きな数に対しても実用的です。上のプログラムPoweMod.tbtはここです。 これを使って冒頭の例を計算すると

181006655297358 ≡167696422262194 (mod 181006655297359)

が直ちに得られます。