２．べき乗

ｐ^ｎの計算

べき乗計算のためのアルゴリズムを素直に書けば、以下のようになる。

w=1
for i=1 to n
  w=w*p
next i
print w

もちろん、これでまちがっていない。
しかし、このアルゴリズムの処理時間はｎの値に比例する。
従って、ｎが極端に大きい場合、例えば、ｎ＝10¹⁰⁰のような場合には、完全にお手上げである。
このループの回数を減らす方法を考える。

例えば、２¹⁰⁰の計算を考えてみる。 上のやり方なら、掛算を100回やらなければならない。 しかし、

２¹⁰⁰＝(２⁵⁰)²

であるから、２⁵⁰を上の方法で計算してやると、掛算回数は50回。
よって、２¹⁰⁰は51回の掛算回数で求められる。
これは、最初の100回に比べて、ほぼ半分の回数である。

ならば、２⁵⁰自体についても同じようにやってみる。

２⁵⁰＝(２²⁵)²

だから、掛算回数は、27回となる。

２²⁵＝(２¹²)²×２

とすると、掛算回数は、12＋１＋１＋１＋１＝16回。これをずっと押し進めると、

２¹⁰⁰＝(((((2²×2)²)²)²×2)²)²

となり、掛算回数は、８回となる。

これをアルゴリズムとして記述すると、

1. ｗ＝１とする（解を格納するための変数）
2. べきｎを２進展開する
3. 上位ビットから順に見ていく
   ・１ならｗ＝ｗ2×ｐ
   ・０ならｗ＝ｗ2
4. 最下位ビットまで計算したときのｗの値が求める解、ｗ＝ｐnとなる

２進展開を格納するために配列を取っておく必要があるが、その大きさは、

    n = 2x
log n = log 2x = x log 2
∴  x = int((log n)/(log 2))

となる。
また掛算回数は、最大２ｘ＝２int((log ｎ)/(log ２))である。
この値は、ｎの桁数の約２倍と思ってよい。
このアルゴリズムにより、べき数がひじょうに大きい場合のべき乗を計算することができる。

プログラムは、以下のとおり。

 10   ' power
 20   print "p^n"
 30   input "p=";P:W=P:I%=0
 40   input "n=";N:dim A%(int(log(N)/log(2)))
 50   while N>1
 60     N=N\2:inc I%:if res=0 then A%(I%)=0 else A%(I%)=1
 70   wend
 80   for J%=I% to 1 step -1
 90     if A%(J%)=0 then W=W^2 else W=P*W^2
100   next J%
110   print W

UBASICでは、桁数の許す限り、p^n で計算が可能である。

modpow(p, n, m)

ｐ^ｎのかわりに、ｐ^ｎをｍで割った余りが知りたいという場合も多い。
基本的には、上記プログラムで最後の掛算のところで、ｍで割った余りを取ればよい。
プログラムは、以下のとおり。

 10   ' modpow
 20   print "p^n@m"
 30   input "p=";P:W=P:I%=0
 40   input "n=";N:dim A%(int(log(N)/log(2))):input M
 50   while N>1
 60     N=N\2:inc I%:if res=0 then A%(I%)=0 else A%(I%)=1
 70   wend
 80   for J%=I% to 1 step -1
 90     if A%(J%)=0 then W=W^2@M else W=P*W^2@M
100   next J%
110   print W

UBASICでは、modpow(p,n,m) という組み込み関数になっている。
この、２進展開して各ビットの値により処理を決める、 というテクニックは、後の Lucas 列の計算でも使用する。