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: モンゴメリ冪剰余サンプルコード     プログラムメモ

前回同様、除算 a = x / y があったとして
式を変形して、

a = (x * (2^(2k)-1)/y + 2^k) >> 2k  (k:=被除数のbit数)

とすれば、逆数を求める部分((2^(2k)-1)/y)の回復法もまた特殊な被除数のため、簡単なものに変形できる。

Bigint iDiv(Bigint x, Bigint y){
  int xb = x.length(); // 被除数のbit数
  int yb = y.length(); // 除数のbit数
  int n = xb * 2 - 1; // ループ回数

  Bigint a(0), m(0);
  while (n--){
    m |= 1;
    a <<= 1;
    if (m >= y){
      m -= y;
      a |= 1;
    }
    m <<= 1;
  }
  a *= x;
  m = 1;
  m <<= xb;
  a += m;
  return a >> (xb*2-1);
}

この方法は、被除数bit数の倍のループとなり、これが最大の減算回数になります。
しかし誤差が出ないので乗算は１回だけで済みます。
前回の方法はbit差に依存したループ回数なので多桁にも応用できますが、この方法だと被除数が大きいとループ回数が多くなってしまいます。

桁差が小さい時に有利な方法

ニュートン法を使った２進多桁の整数除算

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

---

Tags: プログラムメモ

除算 a = x / y があったとして
回復法を使うと被除数のbit数分の加算・減算を繰り返し行うことになる
式を変形して、

a = (x * ((2^k)/y) >> k  (k:=被除数のbit数)

とすれば、逆数を求める部分((2^k)/y)の回復法は特殊な被除数のため、簡単なものに変形できる。

Bigint iDiv(Bigint x, Bigint y){
  Bigint a(0), m(1);
  int xb = x.length(); // 被除数のbit数
  int yb = y.length(); // 除数のbit数
  int n = xb - yb + 1; // ループ回数(bit差+1)

  m <<= yb-1;
  while (n--){
    a <<= 1;
    if (m >= y){
      m -= y;
      a |= 1;
    }
    m <<= 1;
  }
  a *= x;
  a >>= (xb - 1);

  // 誤差修正(ここでのaは,正解より1or2小さい場合がある)
  m = y * a;
  while(x > m){
    m += y;
    if (x > m)
      a += 1;
  }
  return a;
}

と書ける。
計算量は最大でもbit差程度の加減算と２回の乗算なので
bit差が数100のオーダーならNewton法を使うより有利かも。

ニュートン法を使った２進多桁の整数除算

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

---

Tags: プログラムメモ