mod 偶数の乗算にモンゴメリ乗算を使いたい？

導入

$\overline{M_x}M_x\bmod M_y=1$ とします。

$$x+M_x\left\lparen\overline{M_x}(y-x)\bmod M_y\right\rparen$$

これは Garner のアルゴリズムの計算式です。非負整数を与えてこの式を計算すると、次の条件をすべて満たす非負整数 $z$ が得られます。

• $z\leq x+M_x(M_y-1)$
• $z\bmod M_x=x$
• $z\bmod M_y=y$

$y=0$ とすることで、モンゴメリリダクション (Wikipedia) と呼ばれる操作に変わります。

• $z\leq x+M_x(M_y-1)$
• $z\bmod M_x=x$
• $z\bmod M_y=0$
• $z/M_y\bmod M_x=x\overline{M_y}\bmod M_x$ ( ただし $\overline{M_y}M_y\bmod M_x=1$ )

すなわち、 $x$ を $M_x$ で割った余りを変えずに、 $M_y$ で割り切れるようにする計算式になります。実用では、 $M_x$ を $31$ ビット以下の奇数、 $M_y=2^{32}$ とすることで、剰余算に相当する機能を、定数空間の前計算のもとで $32$ ビットどうしの乗算程度の基本的な演算で実装できます。

$M=2^bM_x$ とします。除数が偶数かもしれない場合（ $M$ で割った余りが知りたい場合）にもモンゴメリリダクションをなんとか利用できないでしょうか？

解決法

(1) CRT 復元

一つの方法は、 ${}\bmod 2^b$ と ${}\bmod M_x$ を分けて計算して最後に Garner のアルゴリズムで復元することです。乗算はシンプルで良いですが、加算が $2$ 回に分かれてしまいます。そもそもモンゴメリリダクションが Garner のアルゴリズムの特殊例だとすれば、そこをもっと活用できそうです。

(2) モンゴメリリダクションの変形

ここで、例えば $M_y=2^{32+b}$ ととって $y=2^{32}x\bmod M_y$ としてみます。

• $z\leq x+M_x(M_y-1)$
• $z\bmod M_x=x$
• $z\bmod M_y=2^{32}x\bmod M_y$

ちゃんと $2^{32}$ で割り切れる形になっていますね。よって、

• $z/2^{32}\bmod M_x=\overline{M_y^\prime}x\bmod M_x$ ( ただし $2^{32}\overline{M_y^\prime}\bmod M_x=1$ )
• $z/2^{32}\bmod 2^b=x\bmod 2^b$

を満たす整数 $z/2^{32}$ が得られるので、目的の除数が奇数の場合とほぼ同様のインターフェースが成り立つようになりました。

この方法の弱点は、 $64$ ビットどうしの乗算（の下位 $64$ ビットを求める操作）を必要とするところです。 掛け算の上位 $64$ ビットを計算する Barrett reduction よりはマシですが、通常のモンゴメリリダクションのように $32$ ビットどうしの乗算（ Intel AVX2 で SIMD に乗る）を使って実装しようとすると非効率的になります。

(3) 32 ビット整数の乗算に収める

ここからさらに $1$ 回乗算を使って、 $x\bmod M_x$ の値を変化させずに $x$ の下位 $b$ ビットを $0$ にします。 $d2^b\bmod M_x=1$ を満たす最小の非負整数 $d$ をとって $x$ を $x+d2^b(x\bmod 2^b)-(x\bmod 2^b)$ に変えます。この状態でさっきの計算（ $y=2^{32}x\bmod M_y$ としてリダクション）をすると、計算途中では下位 $b$ ビットが $0$ であることから ${}\bmod M_y$ の計算が ${}\bmod 2^{32}$ に変わるので、 $32$ ビットの乗算で実装できるようになります。

• $z\leq x+M_x(M_y-1)+2^{2b}(M_x-1)$
• $z\bmod M_x=x$
• $z\bmod M_y=2^{32}x\bmod M_y$

結局アルゴリズムは次に示すようなものになります。

$\text{precalc}(M_x, b):$
$\hspace{20px}r\leftarrow {M_x}^{-1}\bmod 2^{32}$
$\hspace{20px}d\leftarrow 2^{-b}\bmod M_x$
$\hspace{20px}M=2^bM_x$
$\text{reduction}(x):$
$\hspace{20px}p\leftarrow x\bmod 2^{b}$
$\hspace{20px}x^{\prime}\leftarrow \lfloor x/2^b\rfloor +pd$
$\hspace{20px}y\leftarrow 2^{32-b}p$
$\hspace{20px}t\leftarrow \left\lfloor\left\lparen x^{\prime}+M_x\left\lparen r(y-x^{\prime})\bmod 2^{32}\right\rparen\right\rparen /2^{32-b}\right\rfloor$
$\hspace{20px}\text{repeat 2 times : if }t\geq M\text{ : }t\leftarrow t-M$

雑に実測

乗算 mod の計算速度を計測します。色々面倒なので、適当に作った C++ プログラムを AtCoder のコードテストで走らせてコードテストの結果で返ってくる時間を記録しました。計測は $2$ 周やりました。基本型の剰余演算子を使わない場合は、 modint 構造体を定義して実装しました。計算内容は下の通りです。

#include <cstdio>

int main(){
    int n; scanf("%d", &n); // 300
    int m; scanf("%d", &m); // 2147483192 = 268435399(prime) * 2^3
    long long ans = 0;
    for(int i=1; i<=n; i++){
        int p = i, k = 1, interval = 2;
        for(int j=1; j<1'000'000; j++){
            k += interval; if(k >= m) k -= m;
            p = (long long)p * k % m;
        }
        ans += p;
    }
    printf("%lld\n", ans);
    if(ans != 304223001390){ printf("fatal\n"); fflush(stdout); }
    return 0;
}

atcoder::dynamic_modint の制約を超過していますが、正しく計算できているようなのでその件は無視します。

記録は次の通りです。

計測の条件がモンゴメリ乗算に有利である割には、あまり速くないですね。下手すると剰余演算子を使ったほうが速いケースも考えられそうです。

結論

$128$ ビット整数型を使える環境であれば atcoder::dynamic_modint を使うほうがよさそうです。

GitHub - atcoder/ac-library: AtCoder Library

AtCoder Library. Contribute to atcoder/ac-library development by creating an account on GitHub.

github.com

また、並列化したいわけでなければ、 $64$ ビットの乗算の下 $64$ ビットを求める操作を無理やり減らす意味はなさそうです。

メリットとして、実装に必要な演算が基本的なので、並列化で高速化させやすいです。また、 $128$ ビット整数型を使用すれば $64$ ビットの mod の上で高速に乗算することができます。

いい文献ないかなぁ