🧮

有限体の実装3(Montgomery乗算の紹介)

2023/03/01に公開

初めに

今回は有限体の山場、Montgomery乗算を紹介します。
Montgomery乗算は普通の乗算の代わりとなる重要な演算です。最初にMontgomery乗算のPythonコードを紹介し、それが持つ数学的な性質を考察します。
記事全体の一覧は有限体の実装一覧参照。

Montgomery乗算

素数 p の有限体の元 x, y において乗算は xy \bmod{p} です( p で割った余り)。
しかし p で割る処理は重たいので、それを避ける方法の一つがMontgomery乗算です。

しばらく数学的な話が続きます。まずいくつか記号の準備をします。
L をCPUのレジスタサイズ(64ビットなら64, 32ビットなら32)、M=2^L とします。以下L=64と仮定します。
N は素数 pL ビット整数で表現したときの配列の大きさです。
Mp は互いに素なので M'M - p' p = 1 となる整数 0 < M' < p0 < p' < M が存在します。
p' \equiv -p^{-1} \pmod{M}, M' \equiv M^{-1} \pmod{p} という意味です。

次のコードで示す関数 mont(x, y) をMontgomery乗算と呼ぶことにします。
もっと簡潔に記述できる部分もありますが、C/C++で実装しやすいようにやや冗長に書いています。コード中のipはp'と読み替えてください。

def mont(x, y):
  MASK = 2**L - 1
  t = 0
  for i in range(N):
    t += x * ((y >> (L * i)) & MASK) # (A)
    q = ((t & MASK) * self.ip) & MASK # (B)
    t += q * self.p # (C)
    t >>= L # (D)
  if t >= self.p: # (E)
    t -= self.p
  return t

コードの説明をしましょう。

(A)

x, yがN桁のuint64_tの配列とすると (y >> (L * i)) & MASKy[i] を表します。
つまりC++なら x * y[i]mulUnitを使って計算し、tに足します。

(B)

(A)で計算したtの下位64ビットにp'(=ip)を掛けてその下位64ビットをqとします。
p'もuint64_tなのでここは

uint64_t q = t[0] * ip;

と同じです。

(C)

mulUnitを使って p と q を掛けて t に足します。

(D)

t を 64ビット右シフトします。

(E)

ループを抜けたら t が p より大きければ p を引いて値を返します。

Montgomery乗算の意味

さて、mont(x, y)が数学的に何をしているか考察しましょう。

定理 : mont(x, y) = x y M'^{N} \bmod{p} である。

{}\bmod{p}で一致することの証明

まず準備として整数 t に対する補助関数 f(t) = t + ((t * p') \bmod{M})p を用意します。この関数はPythonコードの(B), (C)に相当します。
{}\bmod{M} で考えると f(t) \equiv t + (t * p')p = t(1 + p'p) = t M' M なので f(t)M で割れます。
{}\bmod{p} で考えると f(t) \equiv t なので f(t)/M \equiv t/M \equiv t M' です( pM が互いに素なので)。

さてPythonのコードに戻ります。y = \sum_{i=0}^{N-1} y_i M^i とします。
最初 t = 0 で (A)から(D)まで実行すると t=x y_0 に対する f(x y_0)M' = x y_0 M' を計算したことになります。

次のループでは t = x y_0 M' について同様のことをすると t = (x y_0 M' + x y_1)M' = (x y_0 + x y_1 M)M'^2 となります({}\bmod{p} を省略)。
これを N 回繰り返すと t = (x y_0 + x y_1 M + \cdots + x y_{N-1} M^{N-1})M'^N = xy M'^{N} です。

(E)を実行するだけで両辺が一致することの証明

最後に mont(x, y) が x y M'^{N} \bmod{p} に一致することを確認しましょう。
i回目のループのときにステップ(A)に入る直前は t \le 2p-1 とします(i=0のときはt=0なので成立)。

  • (仮定) t \le 2p-1
  • (A)の右辺の最大値は (p-1)(M-1) なので t \le 2p-1 + (p-1)(M-1)
  • (B)の右辺の最大値は M-1
  • (C)の計算後は t \le 2p-1 + (p-1)(M-1) + (M-1)p=(2p-1)M
  • (D)の計算後は t \le (2p-1)M / M = 2p-1

よってループを何度繰り返しても t \le 2p-1 が成り立ちます。
したがって t \bmod{p} を計算するには p 以上だったら p を引く(E)だけで x y M'^{N} \bmod{p} を計算できます。

以上で定理の証明が終わりました。

Pythonによるサンプルコード

mont.pyにmont(x, y)のサンプルコードを置きました。
いくつかの素数 p に対して、上で紹介したmont(x, y)と x y M'^N \bmod{p} を計算する mont_explicit が同じ結果を得ることを確認しています。

p'の求め方

最後に素数 pM=2^{64} が与えられたときに M'M - p'p = 1 となる p' の計算方法を紹介します。
M=2^{64} なので p'p = M'M -1 \equiv M-1 \pmod{M}p'p の下位64ビットがすべて1となる p' を求めるということです。
下位64ビットしか見ないので素数 p の下位64ビット X だけを考えます。
X p'X を左シフトしながら p' のビットが立っているところだけ自身に加算(xor)しながら計算します。

X p'の計算

今求めたい p' は、積の下位64ビットがすべて1となるような値ですから、途中の変数をtとして、tのiビット目が0なら対応するp'のビットを立ててtにX<<iを加算する、ということをくりかえせばOKです。
p'が求まれば、M' = (1 + p' p)/MによりM'も求まります。

def getMontgomeryCoeff(p):
  X = p & (M-1)
  ret = 0
  t = 0
  for i in range(L):
    if (t & (2**i)) == 0:
      t += X << i
      ret += 1 << i
  return ret

まとめ

PythonによるMontgomery乗算のコードと数学的な性質の紹介をしました。次回はこの関数を用いて実際に計算する方法とC++による実装を紹介します。

GitHubで編集を提案

Discussion