【競プロ】繰り返し二乗法と再帰関数

🕒 2019/11/01 🔄 2023/05/01

ここでは、 $n$ 乗の計算を速く行う、繰り返し二乗法について見ていきます。いろいろ書き方はありますが、再帰関数を使う方法を見ます。競プロ記事の一覧はこちら。

📘 目次

整数乗を再帰関数を使って計算する

ここでは、 $a$ の $n$ 乗を $p$ で割った余りを求めることを考えていきます。 $a,n$ は正の整数とします。 $p$ は2以上の整数ですが、一旦、 $10^9+7$ としておきます。これは、競プロでよく使われる設定です。

整数の整数乗でも見ましたが、 $a^n$ とは $a$ を $n$ 回掛ければいいのでした。また、余りの計算でも見たように、 $a^n$ を $p$ で割った余りを求めるには、掛けるたびに余りを考えればいいのでした。

そのため、for文を使うと、次のように求められます。

#include <iostream>
using namespace std;

int main() {
  long long a, n, p = 1e9 + 7, ans = 1;
  cin >> a >> n;
  for (int i = 0; i < n; i++) ans = (ans * a) % p;
  cout << ans;
  return 0;
}

例えば 2 3 と入力すれば、 $2^3=8$ が返ってくることになります。また、2 1000000006 と入力すれば、少し時間はかかりますが、1 と返ってきます。実は、2 の部分を $10^9+7$ 未満の別の正の整数に変えても、1 が返ってきます。これは壊れているのではなくて、実際に 1 なんです。この不思議な現象については、別の機会に見ることになります。

このコードを再帰関数を使って書き直しましょう。和や階乗と再帰関数で見た内容を踏まえれば、 $a^n$ は、 $a$ と $a^{n-1}$ との積なので、次のようなコードを思いつくかもしれません。

#include <iostream>
using namespace std;

int modPow(long long a, long long n, long long p) {
  return (n == 1)? a: a * modPow(a, n - 1, p) % p;
}

int main() {
  long long a, n, p = 1e9 + 7;
  cin >> a >> n;
  cout << modPow(a, n, p);
  return 0;
}

n が 1 なら a をそのまま返し、1 より大きいなら、 a と n - 1 乗したものとの積を返す、という流れです。

これで、 $a^n$ を $p$ で割った余りが計算できます。しかし、この計算では、 $n$ 乗を計算するのに、 $n$ 回程度の計算が必要となります。また、再帰の深さは、1万程度なら問題なくても、もっと深くなると制限を超えてしまいます。先ほどのように n に 1000000006 を代入すると答えは返ってこないでしょう。

実はこの計算回数・再帰の深さはもっと減らすことができます。以下でその方法を見ていきます。

繰り返し二乗法

例えば、紙を8枚に分けるとき、1枚1枚切っていく人はいないと思うんですよね。まずは半分に切って、重ねてさらに半分に切って、重ねてもう一度半分に切れば、8枚に分けることができます。1枚1枚切っていけば7回ですが、重ねて切れば3回で済みます。

この発想で、 $2^8$ を速く計算する方法を考えてみましょう。これを式で書けば次のようになります。\[ 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2 \]これを、初めから順番に掛けずに、前半4つと後半4つに分ければどうでしょうか。\[ (2\times 2\times 2\times 2) \times (2\times 2\times 2\times 2) \]こうすると、前半の4つを計算すれば、後半の4つは使いまわせます。つまり、\[ (2\times 2\times 2\times 2)^2 \]を計算すればいいことになります。

さらに、このカッコの中も、前半と後半に分ければ、\[ \{ (2\times 2)^2 \}^2 \]と変形できます。 $2\times 2$ を計算し、それを2つ掛け合わせ、その答えをさらに２つ掛け合わせることで、 $2^8$ が計算できる、ということです。たったの3回で答えが求められるようになります。

他の指数ではどうでしょうか。例えば、 $2^{100}$ であれば、同じようにして $2^{50}$ を2回掛ければいいですね。50回掛けたもの同士を掛ければ、100回掛けたことになるからです。また、 $2^{50}$ は、 $2^{25}$ を2回掛ければいいでしょう。しかし、 $2^{25}$ のように、指数が奇数になってしまうと困りますね。こうなると、同じものを2回掛ける、というテクニックが使えません。

しかし、分解の仕方を変えれば大丈夫です。今考えている方法では、偶数が出てくればうれしいので、 $2^{25}$ を $2\times 2^{24}$ というように分解すればいいですね。こうすれば、後半の部分は、また小さく分解することができます。

ここまでの内容をもとに、先ほどの modPow(a, n, p) を書き直してみましょう。

n が偶数であれば、 $a^n$ を $a^{\frac{n}{2} } \times a^{\frac{n}{2} }$ だと考えればいいのでしたね。このことから、 modPow(a, n / 2, p) を2回掛けるようにすればいいことがわかるでしょう。 n が奇数なら、 $a^n$ を $a\times a^{n-1}$ と分けて、a * modPow(a, n - 1, p) を計算すればいいです。

こうすると、指数に対応する部分は、偶数なら半分になり、奇数なら $1$ 減ることになります。考えるべき指数はどんどん値が小さくなり、最終的にはすべて $1$ にたどりつきます。modPow(a, 1, p) は、ただの a ですね。

こうして、次のように書き直すことで、 $a^n$ を $p$ で割った余りが速く求められるようになります。

#include <iostream>
using namespace std;

int modPow(long long a, long long n, long long p) {
  if (n == 0) return 1; // 0乗にも対応する場合
  if (n == 1) return a % p;
  if (n % 2 == 1) return (a * modPow(a, n - 1, p)) % p;
  long long t = modPow(a, n / 2, p);
  return (t * t) % p;
}

int main() {
  long long a, n, p = 1e9 + 7;
  cin >> a >> n;
  cout << modPow(a, n, p);
  return 0;
}

具体的にどのように処理されるかを見てみましょう。

2 5 と入力したとします。つまり、 $2^5$ を計算する、ということです。このとき、modPow(2, 5, p) では、 n % 2 == 1 が TRUE なので、2 * modPow(2, 4, p) を実行しようとします。これは、 $2 \times 2^4$ を求めようとする処理です。

modPow(2, 4, p) では、if 文をすり抜けて、t = modPow(2, 2, p) を実行しようとします。ここでもう一度関数が呼び出されます。この呼び出しのときには、modPow(2, 1, p) を実行しようとして、また呼び出されます。このタイミングでは、n == 1 だから 2 が返るようになります。

$2^1=2$ だとわかったことから、 $2^2=(2^1)^2$ が計算され、さらに $2^4=(2^2)^2$ が計算され、最後に $2^5=2\times 2^4$ が計算され、最終的な答えが得られます。

なお、if (n == 0) return 1; は0乗にも対応したい場合に書きます。 $a^n$ の $n$ が正の整数だけなら不要です。

今までの再帰関数よりも動きが複雑ですが、計算の回数がどんどん半分になっていくため、メリットはすごく大きいです。 $2^{10}=1024$ なので、だいたい1000乗を計算する回数が10回程度まで減ることになります。 $10^9$ 乗程度なら、少なければ30回ほどで計算できます。これはすごいスピードアップです。 $10^9$ 乗を計算する機会は競プロでは出てきますが、 $10^9$ 回掛けていると絶対に間に合いません。ここで見たテクニックを使って、計算回数を減らす必要があります。

2 1000000006 と入力してみましょう。先ほど再帰関数を使ったときは結果が返ってきませんでしたが、今回は返ってきます。しかも、for 文を使っていた時よりも、だいぶ速いはずです。

この内容を踏まえれば、AOJ NTL_1_B Powerができるでしょう。

ここで見たように、偶数乗の計算を計算するときに、指数を半分にして2回掛けることで計算回数を減らしていく方法を、繰り返し二乗法といいます。バイナリ法ともいいます。ここで見た書き方はあまり標準的ではないですが、根底にある考え方は同じです。