Pebble Coding

ソフトウェアエンジニアによるIT関連技術や数学の備忘録

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}
右辺の4つの式を評価する。
 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)
 + ...

(式C)以降の部分は、モジュラ計算で消えてしまうので、計算を省略できる。
計算する項の数は (1 + 2 + ... + k) + k = k(k+1) / 2 + k = \tbinom {k+1} {2} + kとなる。

おまけ

mのMSBの値 b^{k-1}の桁の値が小さく2mが b^{k}よりも小さい場合は、 b^{k-1}の桁のr1-r2の引き算は省略できる。

プライバシーポリシー

お問い合わせ

スポンサーリンク