美しきBarrett乗算よ（詠嘆）

はじめに

Barrett乗算を実装してみて、その美しさに気付いたので共有します。とくに、以下の2点が顕著です。

• ゼロ除算を除く任意のu64 % u32を高速に計算できる。当然、u32 * u32 % u32も計算できる。
• Lemireの高速剰余アルゴリズムとマジックナンバーを共有している。

他にも優れた点はたくさんあります。

注意事項

• 多倍長の演算を避けるため、64bitマシンで32bitの乗算を考えます。
• 1での除算には追加の処理が必要なので除外します。対応の一例も紹介します。

Barret乗算のここがすごい

オーバーヘッドが小さい

struct Barrett {
    // 64bitで持っておくと便利
    modulus: u64,
    magic: u64,
}

impl Barrett {
    fn new(modulus: u32) -> Self {
        // 後述のように、1除算は特別な対応が必要。必要ないので切り捨てる。
        assert!(modulus > 1);

        let modulus = modulus as u64;
        // 詳細は後述。1除算の場合のみ、オーバーフローして0になる。
        let magic = (u64::MAX / modulus).wrapping_add(1);

        Self { modulus, magic }
    }
}

Barrett乗算はマジックナンバーを1つしか必要としません。Montgomery乗算とPlantard乗算は2つ必要です^[1]。さらに、マジックナンバーの計算は除算と加算をそれぞれ1回だけです。2回使えば損失はほとんどゼロです。

時間的にも空間的にもオーバーヘッドは僅かです。

あらゆる32bit剰余乗算ができる

実装例の前にアルゴリズムの概要を説明します。Barrett乗算のマジックナンバーは、法を$M \in [1, 2^{32})$として、次のように定義されます。

$$R = \left\lceil \frac{2^{64}}{M} \right\rceil = \left\lfloor \frac{2^{64} - 1}{M} \right\rfloor + 1$$

任意の整数$A \in [0, 2^{64})$について、次が成り立ちます。

$$\left\lfloor \frac{A}{M} \right\rfloor \le \frac{A R}{2^{64}} < \left\lfloor \frac{A}{M} \right\rfloor + 2$$

つまり、$\lfloor A / M \rfloor$を高々1の誤差で得られます。

証明

定義より、次式が成り立ちます。

$$\begin{align*} A &= B M + C \, (0 \le B, 0 \le C \lt M)\\ R M &= 2^{64} + D \, (0 \le D \lt M) \end{align*}$$

したがって、

$$\begin{align*} A R &= (B M + C) R = B (2^{64} + D) + C R\\ B D + C R &\le (B + R)(M - 1)\\ &= (A - C) + (2^{64} + D) - (B + R)\\ &= (A + 2^{64}) + D - (B + C + R) \end{align*}$$

明らかに、$A R \ge B 2^{64}$が成り立ちます。また、$A < 2^{64}$かつ$D < M < R$より、$A R < (B + 2) 2^{64}$です。

以上より、Barrett乗算は次のように書けます。

impl Barrett {
    fn residue64(&self, x: u64) -> u64 {
        // `magic + x + 0`の上位64bit。widening_mul()はまだnightly。
        let quot = self.magic.carrying_mul(x, 0);
        // `quot`は32bitなので、オーバーフローしない。
        // アンダーフローしたら、`quot`を1だけ大きく見積もっていたことになる。
        let (residue, b) = x.overflowing_sub(quot * self.modulus);

        if b { residue.wrapping_add(self.modulus) } else { residue }
    }

    fn mul(&self, a: u32, b: u32) -> u32 {
        self.residue64(a as u64 * b as u64) as u32
    }
}

1除算の対応

1除算の場合、$R = 2^{64}$がオーバーフローしてゼロになってしまいます。Compiler Explorerに実装例があります。

整数をそのまま使える

Montgomery乗算やPlantard乗算では、剰余と一対一対応した別の値を内部的に利用します。一方、Barrett乗算ではそのような工夫は必要ありません。上記の事実と合わせて次のような性質を導けます。

• 特別な初期化処理が要らない。
• 自由に法を変えることができる。

Lemireの高速剰余アルゴリズムをゼロコストで使える

Barrett乗算は商を近似することで剰余を計算します。これには正規化のコストがかかります。Lemireの高速剰余アルゴリズムでは、商を介することなく剰余を直接求めます。おもしろいことに、マジックナンバーは同じです。

64bitマシンではu32 % u32を乗算2回で計算できます。さらに、乗算1回でu32 % u32 == 0を計算できます。

impl Barrett {
    fn residue32(&self, x: u32) -> u32 {
        // `residue64()`では商を近似するために、`hi`を利用していた。
        // 剰余を計算するには`lo`を使えば良い。
        let lo = self.magic.wrapping_mul(x as u64);
        // `lo * modulus + 0`の上位64bit
        lo.carrying_mul(self.modulus, 0).1 as u32
    }

    fn can_divide(&self, x: u32) -> bool {
        let lo = self.magic.wrapping_mul(x as u64);
        // 1除算では`magic`がオーバーフローしているので、
        // `lo <= self.magic.wrapping_sub(1)`とする。
        lo < self.magic
    }
}

Barret乗算のここがちょっと……

作業メモリが大きい

Barret乗算乗算はu32 * u32 % u32の計算をするために、96bitのスペースを利用します。そのため、32bit環境では多倍長の計算をせざるを得ません。一方、Montgomery乗算は64bitのスペースで良いため、mul_high()を利用して効率的な実装が可能です。64bit環境でu64 * u64 % u64をするときも同じです。

次のような使い分けが可能です。なお、数値は64bit環境での値です。

手法	法	特徴
Barrett乗算	$[2, 2^{32})$	偶数を法にとれる
Plantard乗算	$\lesssim 2^{31.3}$の奇数	速い
Montgomery乗算	u64の奇数	大きな法をとれる

乗法逆元の計算に除算が必要

Montgomery乗算やPlantard乗算では法が奇数なので、2の乗法逆元が存在します。このため、拡張Binary GCDアルゴリズムを高速化できます^[2]。Barrett乗算では偶数の法が許されるので、このような高速化はできません。拡張Binary GCDは比較的遅く、拡張Euclid互除法の方が速い傾向にあります。

まとめ

作業メモリの大きさが致命的ですが、Barrett乗算は美しいアルゴリズムです。Lemireの高速除算アルゴリズムとのマジックナンバーの一致に気付いたときは、感動してしまいました。

以下はコメントを受けて改稿しました。詳細はそちらを確認してください。(2026-04-27)

最後に少しだけ一般論を扱います。u{P} % u{Q}の除算で、マジックナンバーを$\lceil 2^{R} / M\rceil$とします。このとき、商を高々1の誤差で求めるには$R = \max(P, 2Q)$、剰余を厳密に求めるには$R = P + Q$が必要最小限です。$P \ge Q$という（自然な）仮定をおくと、剰余のマジックナンバーが商のそれよりも大きくなります。とくに、$Q = P$のときは$1.5$倍になります。ところで、商を誤差なしで求めるにも$R = P + Q$が必要です^[3]。商の誤差は修正が効くので精度を妥協してもよいという訳です。

脚注

1. Plantard乗算で入力を規格化するためにLemireのアルゴリズムを利用するなら追加で1つ必要。

2. Euclid互除法に対するBinary GCDと同程度の高速化が期待されます。

3. コメント欄を参照してください。