バレットの還元アルゴリズム

合同算術において、バレットの還元アルゴリズムは P.D. バレットによって 1986 年に開発されたアルゴリズムである。以下の式を計算することを考える。

$c=amodn$

愚直な方法は、高速な除算アルゴリズムを使用する方法である。バレットの還元アルゴリズムは、除算を乗算に置き換えることにより、$n$ が定数で $a<n^2$ であるときにこの計算を高速化する。

着想

$n$ の浮動小数点数としての逆数を $s=1/n$ とおくと、$\lfloor x\rfloor$ を床関数として

$amodn=a-\lfloor as\rfloor n$

となる。$s$ が十分な精度で計算される限り、これは誤差なく成り立つ。

1ワード版のバレットの還元アルゴリズム

バレットが最初に考案したのは、上記のアルゴリズムの整数版であり、特に値がマシンのワードに収まる場合に適用できる。

上記の方法で $amodn$ を整数の範囲内で計算する場合、自明に対応するのは $n$ による除算を用いる方法である。

func reduce(a uint) uint {
    q := a / n  // 小数部分は暗黙に切り捨てられる
    return a - q * n
}

しかし、除算の計算には時間がかかり、CPUによっては定数時間で完了しない命令が使用されうる。そこで、バレットの方法では $2^k$ による除算なら低コストな右シフトとして実現できることに注目し、$1/n$ を $m/2^k$ で近似する。

$2^k$ が与えられたとき、最良の $m$ を見つけるために以下の式を考える。

$\frac{m}{2^k}=\frac{1}{n}⟺m=\frac{2^k}{n}$

$m$ を整数にするために、$2^k/n$ をどうにか整数に丸める必要がある。最も近い整数に丸めるのが最適だが、$m/2^k$ が $1/n$ より大きくなるとアンダーフローを生じうるため、一般には $m=\lfloor 2^k/n\rfloor$ を用いる。

すると、前述の関数は以下のように書き換えられる。

func reduce(a uint) uint {
    q := (a * m) >> k // ">> k" は k ビットの右シフト
    return a - q * n
}

ここで、$m/2^k\le 1/n$ なので、この関数における q は小さめに評価される。これによりaは $[0,n)$ に必ずしも収まらず、$[0,2n)$ の値を取りうる。そのため、値が $n$ を超えた場合は $n$ を引くことによってこれを補正する。

func reduce(a uint) uint {
    q := (a * m) >> k
    a -= q * n
    if n <= a {
        a -= n
    }
    return a
}

$m/2^k$ は近似にすぎないので、近似が正当となる $a$ の範囲を考える必要がある。$1/n$ の近似誤差は

$e=\frac{1}{n}-\frac{m}{2^k}$

なので、qの誤差は$ae$である。$ae<1$、すなわち $a<1/e$ である限り、この還元は正当である。$a\ge 1/e$ である場合も必ずしもこの関数が誤った値を返すわけではないが、$a$ を先述の範囲内に収めることで一般の場合においての正当性を保証できる。

$k$ を大きく取ると、還元が正当となる $a$ の範囲は広がる一方、他の計算においてオーバーフローを生じる可能性がある。

例

16 ビット整数のシステム上で $n=101$ としたときを考える。

アルゴリズムが意味をなすためには、少なくとも $k$ は $7$ 以上である必要がある。さもなくば、剰余を必要としない小さい値に対してしか還元が正しく動作しなくなってしまう。$k=7$ のときは、$m=\lfloor 2^k/n\rfloor =\lfloor 128/101\rfloor =1$ であり、$k=8$ のときは $m=\lfloor 256/101\rfloor =2$ である。後者においては $1/101$ を $2/256$ で近似することになるので、前者における $1/128$ による近似と全く等価になる。$k=9$ とすれば $m=\lfloor 512/101\rfloor =5$ となり、よりよい近似を得る。

$k=7$ と $k=9$ において、アルゴリズムが正しい値を与える入力の範囲を考える。前者においては、$e=1/n-m/2^k=1/101-1/128=27/12928$ なので、$a<1/e$ すなわち $a<478.81$ である。$a$ は整数なので、実質的な上限は 478 である。（実際は 504 まで正しく動く。）

$k=9$ においては、$e=1/101-5/512=7/51712$ なので $a$ の実質的な上限は 7387 である。（実際は 7473 まで正しく動く。）

近似がより良くなる次の $k$ は 13 であり、$81/8192$ を与える。もっとも、計算途中に現れる $ax$ が $810\le a$ でオーバーフローすることに注意すれば、この問題設定では $k=7$ とした方が良い。

特定の k に対する証明

$2^{k_0}>n$ なる最小の整数 $k_0$ をとり、$k$ を $k_0+1$ とすると、これは上記の式で説明した「意味をなす」$k$ である。前述のコードに対応して次の値を定義する。

• $q=\lfloor \frac{ma}{2^k}\rfloor$
• $r=a-qn$.

床関数の定義から、$q$ は整数で、$r\equiv a(\mathrm{mod}n)$ である。また、$a\le 2^k$ なら $r<2n$ となる。このとき、前述のコードを数式で表すと

$amodn=\{\begin{array}{cc}r& \text{if }r<n\\ r-n& \text{otherwise}\end{array}$

このとき、$r<2n$ であることが以下のように示せる。

$a\le 2^k$ なら

$\frac{a}{2^k}\cdot (2^{k}modn)<n$

である。$n\cdot \frac{mamod{2}^k}{2^k}<n$ が $a$ にかかわらず成立するので、以下の式変形を得る。

$\frac{a\cdot (2^{k}modn)}{2^k}+n\cdot \frac{mamod{2}^k}{2^k}<2n$

$a-(a-\frac{a\cdot (2^{k}modn)}{2^k})+\frac{n\cdot (mamod{2}^k)}{2^k}<2n$

$a-\frac{a}{2^k}\cdot (2^k-(2^{k}modn))+\frac{n\cdot (mamod{2}^k)}{2^k}<2n$

$a-\frac{na}{2^k}\cdot \left(\frac{2^k-(2^{k}modn)}{n}\right)+\frac{n\cdot (mamod{2}^k)}{2^k}<2n$

$a-\frac{na}{2^k}\cdot \lfloor \frac{2^k}{n}\rfloor +\frac{n\cdot (mamod{2}^k)}{2^k}<2n$

$a-\frac{nma}{2^k}+\frac{n\cdot (mamod{2}^k)}{2^k}<2n$

$a-\left(\frac{ma-(mamod{2}^k)}{2^k}\right)\cdot n<2n$

$a-\lfloor \frac{ma}{2^k}\rfloor \cdot n<2n$

$a-qn<2n$

$r<2n$

マルチワード版のバレットの還元アルゴリズム

バレットはこのアルゴリズムをRSA暗号の実装のために考案した。しかし、RSA暗号においては入力値は基本的にワードサイズを超える。そのためバレットは上記のアルゴリズムのマルチワード版も与えた。詳細については Handbook of Applied Cryptography を参照のこと。

多項式に対するバレットのアルゴリズム

バレットのアルゴリズムは多項式除算に対しても適用できる。これは、多項式を反転して X 進数を用いることで計算される。

関連項目

• 類似のアルゴリズムとして、モンゴメリ乗算がある。

参考文献

• テンプレート:Cite book
• テンプレート:Cite book