Hundredsoft

ウィキぺディアでモンゴメリ乗算を検索すると、わかり易く詳しい説明が載っている。
残念ながらこの説明全体をプログラミングしたものが見つからなかったのでC++ソースを書いてみた。

// モンゴメリ乗算剰余演算 c = (a*b) mod nを求める
// gcd(n,r) == 1 となること。(nは奇数であれば良い)
//
Bigint Montgomery(const Bigint& a, const Bigint& b, const Bigint& n) const{
  // n を法とする
  // Bit長を調べる(nb)
  int nb = n.length();

  // r: nより大きな２のべき
  Bigint r(1);
  r <<= nb;

  // r2 = r^2 mod n
  Bigint r2(r*r);
  r2 %= n;

  // n*n2 ≡ -1 mod r
  Bigint n2 = 0; /* 求めるN' */
  Bigint t = 0;
  Bigint vi = 1;
  int ni = nb;
  while (ni--){ /* Rのトップビットを除いたビット数分繰り返す */
    if ((t & 1) == 0){
      /* ゼロになっているビットがあったら、N'のその部分を1にする（NはRと互いに素なので必ず奇数）*/
      t += n; /* 掛け算だが、二進数一桁の掛け算なので実質は足し算 */
      n2 += vi; /* N'のその部分を1にする */
    }
    t >>= 1; /* 必ず端数が出るが切り捨てる */
    vi <<= 1; /* Rは2の冪なので、絶対端数は出ない */
  }
  // ここまでで、今後計算に必要になる、r2,n,nb,n2が得られる。
  // つまりnを法とするモンゴメリクラスを作成するなら
  // 引数nをコンストラクタとするクラスを作成し、
  // r2,n,nb,n2をメンバ変数とする。


  // aのモンゴメリ表現をam, bのモンゴメリ表現をbmとする
  // t は作業領域
  t = a * r2;
  Bigint am = t * n2;
  am.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  am *= n;
  am += t;
  am >>= nb; // 1/Rの意味
  if (am >= n) am -= n;

  t = b * r2;
  Bigint bm = t * n2;
  bm.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  bm *= n;
  bm += t;
  bm >>= nb; // 1/Rの意味
  if (bm >= n) bm -= n;

  // cm: am * bm のモンゴメリリダクション
  t = (am * bm);
  Bigint cm = t * n2;
  cm.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  cm *= n;
  cm += t;
  cm >>= nb; // 1/Rの意味
  if (cm >= n) cm -= n;

  // cmのモンゴメリリダクション
  t = cm;
  Bigint c = t * n2;
  c.hibitmask(nb); // mod Rの意味,(nb)bit以上を0クリア
  c *= n;
  c += t;
  c >>= nb; // 1/Rの意味
  if (c >= n) c -= n;

  return c;
}

ここでは、できるだけwikipediaの説明と合致するように書きましたが、
実利用の際には、r2 = r^2 mod n の計算が重いので、
（これが早ければ、そもそもモンゴメリ演算の必要がない）
nを法とするモンゴメリクラスを作成するのが良いです。
つまりコード中にも書いてあるように、引数nをコンストラクタとしたモンゴメリクラスを作成し、r2,n,nb,n2はメンバ変数にして、モンゴメリリダクション, 乗算剰余演算, 加減算剰余演算, 冪剰余演算等は、モンゴメリクラスのMethodにします。

※ここで書いた乗算は「乗算剰余演算」の項の上段にある方法です。

(注)ここでのBigintクラスは、多倍長整数を格納する仮想のクラスであり、実在するものではありません。

---

Tags: モンゴメリ冪剰余サンプルコード     プログラムメモ