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

２^{181006655297358} (mod 181006655297359）

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

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

冪乗法

　

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

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

５⁴ (mod 9）

を計算するのに，

５¹≡ ５ (mod 9）
５^２≡ ５¹×５＝５×５＝２５≡ ７(mod 9）
５^３≡ ５^２×５＝７×５＝３５≡ ８(mod 9）
５^４≡ ５^３×５＝８×５＝４０≡ ４(mod 9）

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

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

２^{181006655297358} (mod 181006655297359）

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

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

５^２１ (mod 9）

を計算するとします。上の方法だと，２１（２０？）ステップの計算です。冪乗法はこの計算に補助項

５^２≡７(mod 9），５^４≡４(mod 9），５^８≡７(mod 9）， ５^１６≡４(mod 9）

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

５^２０ (mod 9）≡５^１６×５^４×５^１≡４×４×５≡８（mod 9）

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

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

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

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 に入れています。

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

３行から９行が計算の主要部です。４行では指数ｎが２進展開で桁が１かどうかの判定をしています。１の場合，そのときの補助項a を掛けます。

７行では，補助項 a を２乗し次のステップの補助項にします。

８行では，指数ｎを２で割った商をとり，２進展開で１つ上の桁の計算に移ります。

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

123^4567890123  ≡ 108531 (mod 456789)

が得られます。

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

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

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

　例えば，１３×１１（mod 15）の計算を考えます。普通に算すると，１３×１１＝１４３が出てきますから，１５＾２程度の精度が必要です。しかし，これを

１１＝２＾３＋２＋１と２進展開し，

補助項

２^０×１３≡１３ （mod 15），
２^１×１３≡１１（mod 15），
２^２×１３≡２×１１≡７（mod 15），
２^３×１３≡２×７≡１４（mod 15）

を計算して，

１３×１１≡１３×２^３＋１３×２＋１３×１≡８（mod 15）

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

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

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

　整数 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の最大数の２を底とする対数の２乗に比例する程度のステップで す。従って大きな数に対しても実用的です。上のプログラムPoweMod.tbtはここです。 これを使って冒頭の例を計算すると

２^{181006655297358} ≡167696422262194　(mod 181006655297359）

が直ちに得られます。