バイナリ法について色々考える機会があったのでご紹介します。

バイナリ法とは

バイナリ法は、べき剰余演算を効率良く行うアルゴリズムです。RSA暗号のように2048bit（10進数で617桁）といった巨大な数値でも、素早く計算してくれます。

バイナリ法には左向きバイナリ法と右向きバイナリ法の2種類があり、2進展開した指数のビットの評価順が異なります。左向きバイナリ法は文字通り左向き（最小ビットから最大ビット）に評価します。

アルゴリズム

左向きバイナリ法のアルゴリズムは次の通りです。

例

アルゴリズム中の $S$, $T$ の変化を表で確認してみましょう。

表を見ると、ビットが立ったタイミングで $S$ が変化していることがわかると思います。

$T$ が $c$ の○乗部分の値を更新し、ビットが立っている場合にのみ $S$ に「$ST \mod n$」の結果を代入していくことで、最終的に $S$ に「$c^{d} \mod n$」の値が格納される仕組みになっています。

アルゴリズムだけ見ると最初はよく分からないと思うので、ぜひ自分でも簡単な例を考えて表を作ってみてください。実際に書くと、中の動きがよく理解できてくると思います。

実装

C++で実装。実装にあたり、多倍長演算ライブラリであるGMP（GNU MP）を使用しました。

GMPの導入についてはこちらの記事をご覧ください↓

【C/C++対応】UbuntuにGMPをインストールする方法 | IB-Note

C/C++で多倍長演算を行いたいときにとても便利ですが、インストールが若干めんどくさいGMP。久しぶりにUbuntuでインストールしようとしたら思いのほか手間取ってしまったので、備忘録も兼ねてメモしておきます。

itblogger-note.blogspot.com

以下、ソースコードです。

#include <iostream>
#include <random>
#include <gmpxx.h>
#include <chrono>

#define SIZE 2048

using namespace std;
using namespace std::chrono;

mpz_class gen_rand(void);
mpz_class mod_exp(mpz_class c, mpz_class d, mpz_class n);

int main() {
  mpz_class c = gen_rand();
  mpz_class d = gen_rand();
  mpz_class n = gen_rand();
  mpz_class ans;

  //処理時間計測
  system_clock::time_point start, end;

  start = system_clock::now();
  ans = mod_exp(c, d, n);
  end = system_clock::now();

  double time = static_cast<double>(duration_cast<microseconds>(end - start).count() / 1000.0);
  
  cout << "c=" << c << endl;
  cout << "d=" << d << endl;
  cout << "n=" << n << endl;
  cout << "ans=" << ans << endl;
  cout << "time: " << time << "[ms]\n";
  
  return 0;
}

mpz_class gen_rand(void) {
  mpz_class x;
  random_device rnd;
  gmp_randclass rand(gmp_randinit_default);

  rand.seed(rnd());
  x = rand.get_z_bits(SIZE);

  return x;
}

mpz_class mod_exp(mpz_class c, mpz_class d, mpz_class n) {
  //dを2進数に変換
  mpz_class d_dec;
  int i, tmp;
  mpz_class d_bin[2048] = {0};

  d_dec = d;
  for (i = 0; d_dec > 0; i++) {
    d_bin[i] = d_dec % 2;
    d_dec /= 2;
  }
  tmp = i;

  //左向きバイナリ法
  mpz_class S = 1, T = c;
  
  for (i = 0; i <= tmp; i++) {
    if (d_bin[i] == 1) {
      S = (S * T) % n;
    }
    T = (T * T) % n;
  }
  return S;
}

プログラムでは、最初にgen_rand()によって2048bitの整数をランダムに生成し、それらを $c$, $d$, $n$ としてmod_exp(c, d, n)で計算しています。

私の実行環境では、処理時間は平均5.737ms程度でした。素朴なべき剰余演算プログラムだともはや（遅すぎて）計測不可能なことを思うと、バイナリ法がいかにすごいか分かります。

考察

なぜバイナリ法はこれほど早いのでしょうか？ざっくりと計算量について考えてみます。

まず、バイナリ法は $d$ を2進展開する時点でループ回数が半分ずつに減るため $O(\log N)$。続けてべき剰余演算を行いますが、このときのループ回数は $d$ のビット数（$=d^{\prime}$）となりますが $d^{\prime} \ll d$ であり、全体では $O(\log N)$ とみなせます。

一方、素朴なべき剰余演算プログラムだと $d$ が大きくなるのに比例してループ回数も増えるため、$O(N)$ です。この計算量の差は圧倒的ですね。

実際、2048bitだと本来であれば10進数で617桁分のループを回さなければならないところ、2進展開すると2048回（4桁分）のループで済んでしまいます。バイナリ法、恐るべしですね…！