Barrett modular reduction アルゴリズム

モジュラ計算を効率よく行うアルゴリズムの一つに Barrett modular reduction がある。
これを解説してみる。
xとmが与えられたとき、
$r = x \mod m$ を計算したい。

$b \gt 3$ であるbを選び、mをb進数で以下のように表現するとkが定まる。
$m = (m_{k-1} ... m_1 m_0)_b$
ここで $m_{k-1} \ne 0$
$x \lt b^{2k}$ と仮定する。このアルゴリズムではこの場合にしか適用できない。
xは以下のように表現できる。
$x = (x_{2k-1} ... x_1 x_0)_b$

$\mu = floor(b^{2k} / m)$ を事前に計算しておく。
floor(a)はaの小数部を切り捨てて整数にする関数である。

アルゴリズム手順

STEP1:
$q_1 \leftarrow floor(x / b^{k-1})$
$q_2 \leftarrow q_1 \cdot \mu$
$q_3 \leftarrow floor(q_2 / b^{k+1})$

STEP2:
$r_1 \leftarrow x \mod b^{k+1}$
$r_2 \leftarrow q_3 \cdot m \mod b^{k+1}$
$r \leftarrow r_1 - r_2$

STEP3:
もし $r \lt 0$ の場合、 $r \leftarrow r + b^{k+1}$

STEP4:
$r \ge m$ である限り、 $r \leftarrow r - m$ を実行する。

STEP5:
rを返す。

STEP1

$x = Qm + R, 0 \le R \lt m$
となるQ,mが存在する。
このとき、 $Q - 2 \le q_3 \le Q$ が満たされる。

まず $q_3 \le Q$ を確かめる。
以下が成り立つ。
$q_1 \le x / b^{k-1}$
$q_3 \le q_2 / b^{k+1}$
$\mu \le b^{2k} / m$
つなげると、
$q_3 \le q_2 / b^{k+1} = q_1 \mu / b^{k+1} \le q_1 b^{2k} / m \cdot b^{k+1} \le x / b^{k-1} b^{2k} / m \cdot b^{k+1} = x / m = Q + R/m$
$0 \le R/m \lt 1$ で、 $q_3$ とQは整数なので、 $q_3 \le Q$ が成り立つ。

次に $Q - 2 \le q_3$ を確かめる。
以下が成り立つ。
$x/b^{k-1} \le q_1 + 1$
$q_2/b^{k+1} \le q_3 + 1$
$b^{2k}/m \le \mu + 1$
つなげると、
$q_3 + 1 \ge q_2/b^{k+1} = q_1 \mu / b^{k+1} \ge q_1 (b^{2k}/m - 1) / b^{k+1} \ge (x/b^{k-1} - 1) (b^{2k}/m - 1) / b^{k+1} = x/ m - x / b^{2k} - b^{k-1}/m + 1/b^{k+1}$
右辺の４つの式を評価する。
$x / m \ge Q$

$b^{2k} \gt x$ である。変形して、
$- x/b^{2k} \gt - 1$

$b^{k-1} \le m$ が成り立つので変形して、
$- 1 \le - b^{k-1}/m$

右辺の4つの式を評価すると
$x/ m - x / b^{2k} - b^{k-1}/m + 1/b^{k+1} \gt Q - 1 - 1 + 1/b^{k+1} \gt Q - 1 - 1 + 1/b^{k+1}$
最後の2項は1に満たず $q_3, Q$ はともに整数なので、 $q_3 + 1 \ge Q -1$ が成り立つ。

STEP2, 3

$r_1, r_2$ は0以上 $b^{k+1}$ 未満なので、
$0 \le r_1 \lt b^{k+1}$
$0 \le r_2 \lt b^{k+1}$

$- b^{k+1} \lt - r_2 \le 0$ 加えて、
$- b^{k+1} \lt r_1 - r_2 \lt b^{k+1}$ (式A)

$r_1 - r_2 = (x-q_3 m) \mod b^{k+1} = (Qm + R - q_3 m) \mod b^{k+1} = (Q-q_3)m + R \mod b^{k+1}$ (式B)

$q_3 \le Q$ より
$0 \le (Q-q_3)m + R$
$R \lt m$ なので
$(Q-q_3) m + R \lt (Q- q_3) m + m$

$Q-2 \le q_3$ より
$Q - q_3 \le 2$

$(Q-q_3)m + m \le 2m + m = 3m$
したがって、
$0 \le (Q-q_3)m + R \lt 3m$ (式C)

(式A),(式B)より
$r_1 - r_2 \ge 0$ であれば、
$r_1 - r_2 = (Q-q_3)m + R$ となる。
$r_1 - r_2 \lt 0$ であれば、
$r_1 - r_2 + b^{k+1} = (Q-q_3)m + R$ となる。

STEP4, 5

求めたいのはRであるから $(Q-q_3)m + R$ は正なのでmを引いていけばRが求められる。
(式C)より、mを引く回数は最大で2回で済むことが分かる。

どこが効率的になっているのか

bは3より大きければよいので、通常は2の冪乗に取る。
STEP1のbの冪乗の整数の割り算は右へのビットシフト演算と同一なので、計算効率がよい。
STEP2のbの冪乗でのモジュラ計算はビットマスク演算と同一なので、計算効率がよい。
STEP1の $q_2$ の計算は $q_1 \cdot \mu$ の計算を全て行う必要はなく、 $q_3$ の右ビットシフトで消える項は計算を省略できる。(計算量1)
STEP2の $q_3 \cdot \mu$ の計算も同様に全て行う必要はなく、モジュラ計算で消える項は計算を省略できる。(計算量2)

マルチ精度乗算

マルチ精度乗算の繰り上がりについての評価をまずしておく。
$0 \le y \lt b^{k+1}$ の範囲にあるyについて $y = y_k b^{k} + y_{k-1} b^{k-1} + ... + + y_1 b + y_0$ と表現できる。
zも同じ範囲にあるとする。 $y \cdot z$ (式A)を計算する。
(式A) =
$y_0 z_0$
$+ b (y_1 z_0 + y_0 z_1)$
$+ b^{2} (y_2 z_0 + y_1 z_1 + y_0 z_2)$
$+ ...$
$+ b^{k-2} (y_{k-2} z_0 + y_{k-1} z_1 + ... + y_0 z_{k-2})$ = (式B)
$+ b^{k-1} (y_{k-1} z_0 + y_{k-2} z_1 + ... + y_0 z_{k-1})$
$+ b^{k} (...)$
$+ b^{k+1} (...)$
$+ b^{k+2} (...)$
$+ ...$

k<bであるとき、(式B)は $b^{k+1}$ より小さくなるつまり、 $b^{k+1}$ で割り算したとき、0になることを示す。
これを示せば最初の項から(式B)の項までは計算を省略できるということになる。

$y_i \le b-1, z_i \le b-1$ であることと $b \gt 3$ を使って(式B)を評価する。
$(式B) \le b^{k-2} ( (b-1)(b-1) + (b-1)(b-1) + ... + (b-1)(b-1) )$ = $b^{k-2} (b-1)^{2} (k-1) \lt b^{k+1} (1 - \frac {1} {b})^{3} \lt b^{k+1} (1 - \frac {1} {3})^{3} \lt b^{k+1}$
これで示せた。

計算量1

xは $b^{2k}$ 未満なので $q_1 = floor( {x} / b^{k-1})$ は $(y_k ... y_1 y_0)_b$ と表現できる。
mは $b^{k}$ 未満なので $\mu = floor(b^{2k} / m)$ は $(z_k ... z_1 z_0)_b$ と表現できる。
したがって $q_2 = q_1 \mu$ は(式A)のようになる。
$b^{k+1}$ で割り算したときに0になる(式A)とそれ以前のすべての項は計算する必要がない。
単精度乗算の項数は $(k+1)(k+1) - (1 + 2 + ... + (k-1)) = (k+1)^{2} - k(k-1)/2 = (k+1)^{2} - \tbinom {k}{2}$ となる。
計算すべき和は
$b^{k-1} (y_{k-1} z_0 + y_{k-2} z_1 + ... + y_0 z_{k-1})$
$+ b^{k} (y_{k} z_0 + y_{k-1} z_1 + ... + y_0 z_{k})$
$+ b^{k+1} (y_{k} z_1 + y_{k-1} z_2 + ... + y_1 z_{k})$
$+ ...$
$+ b^{2k-1} (y_{k} z_{k-1} + y_{k-1} z_{k})$
$+ b^{2k} y_{k} z_{k}$
である。

計算量2

$q_1, \mu$ はともに $b^{k+1}$ 未満なので、
$q_2 = q_1 \mu \lt b^{2k+2}$ である。
$q_3 = floor(q_2 / b^{k+1}) \lt floor(b^{2k+2} / b^{k+1}) = floor(b^{k+1}) = b^{k+1}$
したがって、
$q_3 = (y_{k} ... y_0)_b$ と表現できる。
mは $m = (m_{k-1} ... m_0)_b$ と表現される。
$q_3 m$
$= y_0 m_0$
$+ b (y_0 m_1 + y_1 m_0)$
$+ b^{2} (y_0 m_2 + y_1 m_1 + y_2 m_0)$
$+ ...$
$+ b^{k-1} (y_0 m_{k-1} + ... + y_{k-1} m_0)$
$+ b^{k} (y_1 m_{k} + ... + y_{k} m_1)$
$+ b^{k+1} (...)$ (式C)
$+ ...$