繰り返し二乗法によるべき乗(pow(x,n))の計算のアルゴリズム

2020年2月24日数学bit演算,剰余,競プロ,べき乗,ダブリング,繰り返し二乗法,2進数

多くのプログラミング言語でサポートされてる \(x^n\) を計算する関数のアルゴリズムです。

Input : pow(2,4) (x=2, n=4)
Output : 16

べき乗は非常に大きな数値になるため、オーバーフローの原因になることもあります。競技プログラミングなどでは素数で割った余りを計算させることが多く、自身でべき乗計算の実装を行う必要がある場合もしばしばです。

ここでは、べき乗を高速に計算する「繰り返し二乗法」を中心に解説します。

アルゴリズム

整数 \(x\) と \(n\) が与えられた時に、\(x^n\) を計算するアルゴリズムは以下のようになります。

繰り返し二乗法を用いた \(x^n\) のアルゴリズム:

  1. ans = 1 として、以下を \(n\) を2進数で見た時の全ての桁について行う
    • \(n\) の \(i\) 桁目 が 1 ならば \(x^{2^i}\) を ans にかける
  2. ans が答え

※ 繰り返し二乗法では、\( x^1, x^2, x^4, x^8, x^{16}, \cdots \) を上手く組み合わせて計算することを考えます。\(n\) を2進数でみて、下から \(i\) 桁目が 1 なら \(x^{2^i}\) をかけます(\(i\) は 0-indexed)。

※ 繰り返し二乗法はダブリングの一種と捉えることもできます。ダブリングは色々な問題に応用でき、最近共通祖先(LCA)の計算などにも利用されます。

計算量

  • 時間計算量: \(O(\log n)\)
  • 空間計算量: \(O(1)\)

\(n\) のビット数だけループを繰り返すので、時間計算量は \(O(\log n)\) となり、単純に \(n\) 回の掛け算を計算するよりも高速です。

C++での実装例

通常の pow(x,n)

以下のように \(n\) を一桁ずつ削っていくことで、最下位ビットを見るだけですべての桁を確認することができます。

#include <bits/stdc++.h>
using namespace std;

long long pow(long long x, long long n) {
    long long ret = 1;
    while (n > 0) {
        if (n & 1) ret *= x;  // n の最下位bitが 1 ならば x^(2^i) をかける
        x *= x;
        n >>= 1;  // n を1bit 左にずらす
    }
    return ret;
}

int main() {
    long long x, n;
    cin >> x >> n;
    cout << pow(x, n) << endl;
    return 0;
}

剰余の形での pow(x,n)

競技プログラミングなどで大きな素数で割ったあまりなどを求める場合には、以下のように掛け算や足し算をする度に剰余の計算をします。

#include <bits/stdc++.h>
using namespace std;
const int MOD = 1000000007;

long long pow(long long x, long long n) {
    long long ret = 1;
    while (n > 0) {
        if (n & 1) ret = ret * x % MOD;  // n の最下位bitが 1 ならば x^(2^i) をかける
        x = x * x % MOD;
        n >>= 1;  // n を1bit 左にずらす
    }
    return ret;
}

int main() {
    long long x, n;
    cin >> x >> n;
    cout << pow(x, n) << endl;
    return 0;
}

繰り返し二乗法の考え方

「\(1\) に \(x\) を \(n\) 回掛け算する」ような単純なアルゴリズムでは、時間計算量が \(O(n)\) となってしまうため、あまり高速とは言えません。

繰り返し二乗法では、

$$x^1, x^2, x^4, x^8, x^{16}, \cdots $$

を上手く組み合わせて計算することを考えます。

例えば、\(3^8\) の計算では、普通に計算すると 3 を 8 回掛け算することになりますが、

  • \(3^1 = 3\)
  • \(3^2 = 3^1 \cdot 3^1 = 9\)
  • \(3^4 = 3^2 \cdot 3^2 = 81\)
  • \(3^8 = 3^4 \cdot 3^4 = 6561\)

などのように考えると、たったの 4 回で計算をすることができるのです。

\(n\) を 2のべき乗で考える

先程の例から考えると、\(x^n\) について、 \(n\) を以下のように 2 のべき乗の形で表すと上手く考えられそうです。

$$ n = 2^{k_1} + 2^{k_2} + 2^{k_3} + \cdots$$

このようにすると、

\begin{align}
x^n &= x^{2^{k_1} + 2^{k_2} + 2^{k_3} + \cdots} \\
&= x^{2^{k_1}} \cdot x^{2^{k_2}} \cdot x^{2^{k_3}} + \cdots
\end{align}

となるので上手く計算ができます。

例えば、\(3^6\) の計算では、

$$ 6 = 2^{2} + 2^{1}$$

となるので、

\begin{align}
3^6 &= 3^{(2^{2} + 2^{1})} \\
&= 3^{2^{2}} \cdot 3^{2^{1}} \\
&= 3^{4} \cdot 3^{2}
\end{align}

より \(3^{4}\) と \(3^{2}\) が分かれば簡単に計算が可能です。

この \(6\) を2進数で表すと、 \(110_2\) となることに注目すると、

  • n を2進数でみて、下から i 桁目が 1 なら \(x^{2^i}\) をかける(i は 0-indexed)

ことを繰り返せば良いことが分かります。

練習問題

  • [AOJ] NTL_1 Power : 繰り返し二乗法でないと時間制限に間に合いません。剰余の形で求める必要もあります。