[数学] n乗の計算

$a^n$ の計算について整理します。この記事の大部分は、実数のべき乗に限らず、結合法則を満たす積が定義されているもの全てに適用できます。実数以外では、行列のべき乗や多項式のべき乗などが多いでしょうか。

なお、pythonでは組み込み関数のpow()があり、基本的な$a^n$ あるいは$a^n \pmod{m}$といった計算は自分で実装する必要はありません。

バイナリ法

再帰による実装

$a^{101}$を計算するには次のようにします。
・$a^{101}$ を計算するには、$x = a^{50}$ として、$ax^2$ を計算すればよい。
・$a^{50}$ を計算するには、$x = a^{25}$ として、$x^2$を計算すればよい。
・$a^{25}$ を計算するには、$x = a^{12}$ として、$ax^2$ を計算すればよい。
・$a^{12}$ を計算するには、$x = a^6$ として、$x^2$ を計算すればよい。
・$a^6$ を計算するには、$x = a^3$ として、$x^2$ を計算すればよい。
・$a^3$ を計算するには、$x = a$ として、$ax^2$ を計算すればよい。

したがって、$a$ の値が与えられると、2種の計算$x\mapsto x^2$, $x\mapsto ax^2$ を合計 $6$ 回繰り返すことで、$a^{101}$ が計算できることになります。

次々とより小さい値の場合の計算結果を利用するため、再帰関数による実装が簡潔です。

def power(a, n):
    if not n:
        return 1
    x = power(a, n >> 1)
    x *= x
    return a * x if n & 1 else x

ループによる実装

$a^{101}$を計算するには次のようにします。
・2進法により、$101 = 64 + 32 + 4 + 1$と分解できる。したがって、$a^{101} = a^{64}\cdot a^{32}\cdot a^{4}\cdot a$ と計算できます。
・$a^{2n} = (a^{n})^2$ を繰り返し用いると、$a^2$, $a^4$, $a^8$, $a^{16}$, $a^{32}$, $a^{64}$が順に計算できます。

これら2種の計算をループしながら同時処理してあげると、次のようになります。

def power(a,n):
    x = 1
    while n:
        if n & 1:
            x *= a
        n >>= 1
        a *= a
    return x

計算量について

2種の実装を紹介しましたが、どちらも本質的な計算量は変わりません。$n$ を$2$ 進法で見て、bitが立っていると1回、ひとつ桁数が増えるときに1回の掛け算をすることになります。掛け算の回数、if 分岐の回数などはどちらの実装も一致します。$2$進法表記した場合の桁数の $1$ ～ $2$ 倍の掛け算を行うため、計算量は $O(\log n)$ になります。（この際 2数の積をとる操作のコストは定数と考えています。実際には巨大な桁数の自然数などを扱うのであれば、事情は変わってきます。）

while ループより再帰呼出しの方がコストが高いので、後者の方がわずかにパフォーマンスで優れます。私は前者の実装をすることの方が多いです。（後者の実装は、2種の値を同時に更新していくぶん、少し頭を使う量が多いようなイメージでｗ）

また掛け算の回数は変わりませんが、例えば行列のべき乗 $A^n$ の計算において $A$ が疎行列であるならば、前者の実装の方が1回あたりの積のコストは安く済みますね。

事前計算による計算量削減

何度も繰り返し $n$ 乗計算を行う場合には、条件次第では、計算量を削減できます。

平方分割・k乗根分割による方法

例えば、同じ $a$ に対してさまざまな $1\leqq n\leqq 10^7$ に対する $a^n$ を計算することを考えます。このとき、次のような方法があります。

・$1\leqq k\leqq 10^3$ に対して $a^k$ を計算しておく。
・$1\leqq k\leqq 10^4$ に対して $a^{1000k}$ を計算しておく。

合計 $10^4 + 10^3$ の回数の事前計算を行っています。この事前計算のもとで、$a^n$ の計算は
\[ a^{1234567} = a^{1234000} \cdot a^{567}\]
のように行うことができます。バイナリ法を使うだけでは $30$ 回以上の掛け算が必要なところ、たった 1 度の掛け算で済んでいます。

一般には、$AB \geq N$となる$A, B$をとれば、$O(A+B)$の計算回数を行うことで、$1\leq n\leq N$に対する$a^n$ が1件あたり $O(n)$ で計算できます。上のように$A = 1000$ などと決め打っても十分な場合が多いです。また、$AB\geq N$ の条件下で $A+B$ の最小値は$A = \sqrt{N}$ で実現されるため、$A = \lfloor \sqrt{N} \rfloor$ などととることに理論的な妥当性があります。特にそのような場合に平方分割といいます。最終的には

$O(2\sqrt{N})$の事前計算をもとに、1つあたりの$n$乗を$1$回の乗算で計算する

ことに成功しています。なお、同様の考え方を$k$ 乗根分割で行うと、

$O(kN^{1/k})$の事前計算をもとに、1つあたりの$n$乗を$k-1$回の乗算で計算する

ことも可能になります。事前計算量・メモリ容量と、1件あたりの$a^n$ の計算量のトレードオフがあり、目的に応じて $k$ を選ぶのも良いでしょう。なお、$k$ 乗根分割の最も極端な場合として、$k = 1$の場合があります。全部事前計算しておくと、全部乗算なしに値を取り出せる…当たり前（笑）

原始根の利用

平方分割による事前計算は、$a$の値が固定されてしまうデメリットがあります。
小さい素数 $p$ を法とする $a^n \pmod{p}$ を、さまざまな $a,n$ に対して計算したいというのであれば、また異なる方法をとることができます。なお、ここでの考え方は $p$ が素数べきのとき（これも原始根がある）や、合成数のとき（中国剰余定理で素数べきに帰着）にも利用可能ですが、本記事では省略します。

例えば、$a^n\pmod{13}$ の計算を考えてみます。$\pmod{13}$では、$0$ を除くあらゆる整数が $2^k$ と合同になります。

$2^k$ には$1$ から $12$ までの整数が全て現れており、$2^{12} \equiv 1\pmod{13}$ が成り立ちます。一般に$p$ を素数とするとき、ある $r$ が存在して、全ての $1\leqq a\leqq p-1$ は $r^k$ の形で表すことができます。この $r$ を原始根といいます。

さえ、上の表さえあれば、$a^n\pmod{13}$ の計算は、人間ですら一瞬です。例えば$123^{4567}\pmod{13}$ を計算するには、次のようにすればよいです：
\[
 123^{4567} \equiv (123 \pmod{13})^{4567\pmod{12}} \equiv 6^7 \equiv (2^{5})^7\equiv 2^{35}\equiv 2^{11}\equiv 7\pmod{13}. 
\] より一般の mod についても、このような表を計算しておく（$O(p)$の事前計算）ことで、あらゆる$a, n$ に対する $a^n\pmod{p}$ が $O(1)$ の計算時間で計算できます。

周期の利用