64bit数の素数判定

関連スライド:
https://speakerdeck.com/mizar/64bitshu-nosu-shu-pan-ding

素数判定問題

No.3030 ミラー・ラビン素数判定法のテスト https://yukicoder.me/problems/no/3030
実行時間制限: 1ケース 9.973秒 / メモリ制限: 509MB / ソースコードのサイズ制限: 64kiB

問題文
与えられた $n$ 個の正整数 $\{x_0,x_1,\cdots,x_{n-1}\}$ それぞれについて、それが素数かどうかを判定してください。

• $n$ は正整数。
• $1 \le n \le 10,000$
• $1 \le x_i \lt 2^{63}$ (64ビット符号つき整数に収まります)

出力
正整数のそれぞれについて、その値と素数か否か（素数なら 1 そうでないなら 0）を半角空白区切りで、入力順を維持したまま出力し、改行してください。

試し割りで解けるか？

• $2^{63}=9,223,372,036,854,775,808$ (922京3372兆0368億5477万5808)
• $2^{64}=18,446,744,073,709,551,616$ (1844京6744兆0737億9551万1616)
• $1$ は素数ではない、$2$ は素数、$4$ 以上の偶数は素数ではない
• それ以外のケース($3 \le x$ の奇数)なら、 $3 \le k \le \lfloor\sqrt{x}\rfloor$ の奇数 $k$ で全て割り切れなければ素数と判定できるが… $\lfloor\sqrt{2^{63}}/2\rfloor = 1518500249$ … 1万個の整数について、それぞれ約15億回試し割りする必要がある計算
• 仮に1秒で10億回試し割りできるとしても(除算は加減乗算より計算が重いので、実際にはより悪く見積もった方が良い)、10秒以内には解けない
• $3 \le k \le \sqrt{2^{63}}$ の全ての素数(約1.5億個)を表にしてソースコードに羅列も出来ないし(64kiB制限)、その素数表を生成する時間もない、$3 \le k \le \sqrt{2^{63}}$ の素数表が元からあったとしても試し割り法ではやっぱり時間が足りない

この記事における数学記号の補足説明

• $\operatorname{mod}$ : モジュロ(modulo)・剰余演算、割り算における余りの値を取得する演算記号。
  • ja.wikipedia.org - 剰余演算
  • ja.wikipedia.org - 除法
  • ja.wikipedia.org - 整数の合同
  • $\underset{\operatorname{mod}n}{\equiv}$ : 両辺の値が $\operatorname{mod}n$ において合同、つまり両辺の値をそれぞれ $n$ で割った余りが等しい事を示す
  • $\underset{\operatorname{mod}n}{\not\equiv}$ : 両辺の値が $\operatorname{mod}n$ において合同ではない、つまり両辺の値をそれぞれ $n$ で割った余りが等しくない事を示す
  • $\underset{\operatorname{mod}2^{64}}{+},\underset{\operatorname{mod}2^{64}}{-},\underset{\operatorname{mod}2^{64}}{\times}$ : それぞれの演算子（加算記号・減算記号・乗算記号）の左右にある値の入力およびその結果出力をそれぞれ、 $2^{64}$ で割った余りの値を使う事を示す
• $\gcd(a,b)$ : 最大公約数(greatest common divisor)。
  • ja.wikipedia.org - 最大公約数
• $\operatorname{lcm}(a,b)$ : 最小公倍数(least common multiple)。
  • ja.wikipedia.org - 最小公倍数
• $\min(a,b)$ : 最小値(minimum)。
• $\max(a,b)$ : 最大値(maximum)。
  • ja.wikipedia.org - 最大と最小
• $\operatorname{abs}(x)$ : 絶対値(absolute value)。
  • ja.wikipedia.org - 絶対値
• $\sqrt{x}$ : 平方根(square root)。
  • ja.wikipedia.org - 平方根
• $b^p$ : べき乗、冪乗(exponentiation)。
  • ja.wikipedia.org - 冪乗
• $\log(x)$ : 対数(logarithm)。
  • ja.wikipedia.org - 対数
• $\ln(x)$ : 自然対数(natural logarithm)。
  • ja.wikipedia.org - 自然対数
• $\lfloor x \rfloor$ : 床関数(floor function)。実数に対してそれ以下の最大の整数を返す関数。正の数の場合は小数点以下切り捨ての操作と等しい。
  • ja.wikipedia.org - 床関数と天井関数
• $a\simeq b$ : 近似(approximation)。ほぼ等しい。
  • ja.wikipedia.org - 近似
• $\mathcal{O}(\quad)$ : ランダウの記号(Landau symbol)。この記事では計算量の問題を大まかに述べるのに使っている。
  • ja.wikipedia.org - ランダウの記号
• $({a}|{n})$ : ルジャンドル記号(Legendre symbol)・ヤコビ記号(Jacobi symbol)・クロネッカー記号(Kronecker symbol)。 $\displaystyle{\left(\frac{a}{n}\right)}$ のように縦に並べた形式が取られる事も多いが、この記事では横に並べる形式で表記する。
  • ja.wikipedia.org - ルジャンドル記号
    • $n$ が 奇素数の場合に定義された記号。
  • ja.wikipedia.org - ヤコビ記号
    • ルジャンドル記号を拡張して $n$ が 正の奇数の場合に定義された記号。
  • en.wikipedia.org - Kronecker_symbol
    • ヤコビ記号を拡張して $n$ が任意の整数の場合に定義された記号。線形代数などで使われる「クロネッカーのデルタ」とは別物。

mod (剰余) の基本的な性質

同値律

• 反射律: 任意の整数 $a$ に関して
  $a\equiv a\ (\operatorname{mod} n)$;
• 対称律: 任意の整数の対 $a,b$ に関して
  $a\equiv b\ (\operatorname{mod} n)\Leftrightarrow b\equiv a\ (\operatorname{mod} n)$;
• 推移律: 任意の整数の組 $a,b,c$ に関して
  $a\equiv b\ (\operatorname{mod} n)$ かつ $b\equiv c\ (\operatorname{mod} n)$ ならば $a\equiv c\ (\operatorname{mod} n)$.

例えば

• $1\equiv 5\ (\operatorname{mod} 4),\quad 5\equiv 1\ (\operatorname{mod} 4)$
• $5\equiv 9\ (\operatorname{mod} 4),\quad 9\equiv 5\ (\operatorname{mod} 4)$
• $1\equiv 9\ (\operatorname{mod} 4),\quad 9\equiv 1\ (\operatorname{mod} 4)$

代数学的性質

• $a_1\equiv b_1\ (\operatorname{mod} n)$ かつ $a_2\equiv b_2\ (\operatorname{mod} n)$ ならば
  • $a_1+a_2\equiv b_1+b_2\ (\operatorname{mod} n)$
  • $a_1-a_2\equiv b_1-b_2\ (\operatorname{mod} n)$
  • $a_1\times a_2\equiv b_1\times b_2\ (\operatorname{mod} n)$
• $a\equiv b\ (\operatorname{mod} n)$ ならば
  • 任意の整数 $c$ に対して $a+c\equiv b+c\ (\operatorname{mod} n)$
  • 任意の整数 $c$ に対して $a-c\equiv b-c\ (\operatorname{mod} n)$
  • 任意の整数 $c$ に対して $a\times c\equiv b\times c\ (\operatorname{mod} n)$
  • 正整数 $q\gt 0$ に対して $a^q\equiv b^q\ (\operatorname{mod} n)$

単位的環

整数 $n\ge 2$ および 任意の $a, b, c$ において

• 加法の結合性: $(a+b)+c\equiv a+(b+c)\ (\operatorname{mod} n)$
• 加法の可換性: $a+b\equiv b+a\ (\operatorname{mod} n)$
• 加法単位元: $d\equiv 0\ (\operatorname{mod} n)$ ならば $d+a\equiv a+d\equiv a\ (\operatorname{mod} n)$
• 加法逆元: $a+(-a)\equiv(-a)+a\equiv 0\ (\operatorname{mod} n)$ となる $-a$ が存在する
• 乗法の結合性: $(a\times b)\times c\equiv a\times(b\times c)\ (\operatorname{mod} n)$
• 乗法単位元: $e\equiv 1\ (\operatorname{mod} n)$ ならば $e\times a\equiv a\times e\ (\operatorname{mod} n)$
• 左右分配性:
  • $a\times(b+c)\equiv(a\times b)+(a\times c)\ (\operatorname{mod} n)$
  • $(b+c)\times a\equiv(b\times a)+(c\times a)\ (\operatorname{mod} n)$

RSA暗号

桁数が大きい合成数の素因数分解が現実的な計算量では困難であることを安全性の前提とした公開鍵暗号の一つ。1977年に発明された。ja.wikipedia.org - RSA暗号

SSL公開鍵(WebのHTTPS通信など)、SSH公開鍵などで現在も使われている所が多いが、RSA暗号からは楕円曲線暗号(ECDSA, ED25519, ...)などのより強力な暗号方式への移行が徐々に進んでいる。例えば2031年以降、RSA 2048bit公開鍵は使用禁止扱いとなる見込み。 nict.go.jp - 暗号の安全性評価

en.wikipedia.org - RSA_Factoring_Challenge: RSA素因数分解チャレンジ（抜粋）

RSA_Number	10進桁数	2進桁数	素因数分解された日
RSA-100	100	330	1991-04-01
RSA-140	140	463	1999-02-02
RSA-576	174	576	2003-12-03
RSA-768	232	768	2009-12-12
RSA-240	240	795	2019-12-02
RSA-250	250	829	2020-02-28
RSA-260	260	862	未踏
RSA-270	270	895	未踏
RSA-896	270	896	未踏
RSA-1024	309	1024	未踏
RSA-1536	463	1536	未踏
RSA-2048	617	2048	未踏

さまざまな素数判定法

• 決定的素数判定法
  • 試し割り法・エラトステネスの篩 ( $\mathcal{O}(\sqrt{n})$ ) 遅い
  • Miller素数判定法
  • APR(Adleman–Pomerance–Rumely)素数判定法 ( $(\log n)^{\mathcal{O}(\log\log\log n)}$ ) 高速だが多項式時間ではない
  • 楕円曲線素数判定法 ( $\mathcal{O}((\log n)^{5+\varepsilon})$ ) 2022年6月、 $10^{50000}+65859$ の素数判定の記録
  • AKS(Agrawal–Kayal–Saxena)素数判定法 多項式時間で決定的なアルゴリズムとして2002年発表
• 確率的素数判定法 (今回はこちら、$2^{64}$未満の範囲で決定的に判定する)
  • Miller-Rabin 素数判定法 Miller素数判定法を確率的判定法にしたもの
  • Lucas 確率的素数判定法
  • Baillie–PSW 素数判定法 (Miller-Rabin + Lucas の複合)
• 特殊な形の数に固有の素数判定法 $2^{82589933}-1, 10223\cdot 2^{31172165}+1, 3\cdot 2^{18924988}-1$ など
  • メルセンヌ数( $N=2^n-1$ の形の自然数、$n$は自然数 )に対する素数判定法、Lucas–Lehmer test
  • プロス数( $N=k\cdot 2^n+1$ の形の自然数、$n$ は自然数、$k$ は奇数かつ $k\lt 2^n$ )に対する素数判定法
  • $N=k\cdot 2^n-1$ の形の自然数 に対する素数判定法、$n$ は自然数、$k$ は奇数かつ $k\lt 2^n$ Lucas–Lehmer–Riesel test

Miller-Rabin素数判定法

Miller 素数判定法 (決定的素数判定)

$n \ge 3$ の奇数に対して以下のように判定する。

1. $n - 1$ を $2$ のべき乗で割り、 $2^s\cdot d$の形式 (ただし $d$ は奇数) にする。
2. 以下を $[2, \min(n-1,\lfloor 2(\ln n)^2\rfloor)]$ の範囲の全ての base値 $a$ について繰り返す:
   • $a^d\underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}1$ かつ $[0,s-1]$ の範囲の全ての $r$ について $a^{2^rd}\underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}n-1$ なら、composite(合成数) を返す。
3. 2.で composite(合成数) と判定されなかった場合、 prime(素数) を返す。

• 拡張リーマン予想が真であることを前提とした決定的素数判定法のため、根拠がやや弱い。証明されていない予想に依存しない、他の決定的素数判定法もある
• $\lfloor 2(\ln n)^2\rfloor$ の値は $n\simeq 2^{63}$ だと $3813$ 、 $n\simeq 2^{64}$ だと $3935$ 。でも、 $n\lt 2^{64}$ の時に $4000$ 回弱もこの通りに調べる必要があるかと言うと…。

Miller-Rabin 素数判定法(1) (確率的素数判定)

$n \gt 3$ の奇数に対して以下のように判定する。 $k$ は判定の正確度のパラメータ。

1. $n - 1$ を $2$ のべき乗で割り、 $2^s\cdot d$の形式 (ただし $d$ は奇数) にする。
2. 以下を $k$ 回繰り返す:
   • $[2,n-2]$ の範囲から base値 $a$ をランダムに選ぶ。
   • $a^d\underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}1$ かつ $[0,s-1]$ の範囲の全ての $r$ について $a^{2^rd}\underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}n-1$ なら composite(合成数) を返す。
3. 2.で composite(合成数) と判定されなかった場合、 probably prime(おそらく素数) を返す。 (n が合成数なのに誤って probably prime(おそらく素数) を返す確率は $\displaystyle{\frac{1}{4^{k}}}$ 未満)

これは確率的な素数判定法だが、 $a$ の選び方を工夫すると…。

Miller-Rabin 素数判定法(2) (範囲限定:決定論的変形)

$n\lt 2^{64}$ の範囲では7回のループで判定可能な事が判明している。(Jim Sinclair, 2011年)

1. $n - 1$ を $2$ のべき乗で割り、 $2^s\cdot d$の形式 (ただし $d$ は奇数) にする。
2. base値 $a$ の値をそれぞれ $\{2,325,9375,28178,450775,9780504,1795265022\}$ として以下を行う:
   • $a \underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}0$ かつ $a^d\underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}1$ かつ $[0,s-1]$ の範囲の全ての $r$ について $a^{2^rd}\underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}n-1$ なら composite(合成数) を返す。
3. 2.で composite(合成数) と判定されなかった場合、 probably prime(おそらく素数) を返す。 (少なくとも $n\lt 2^{64}$ の範囲で probably prime(おそらく素数) を返した場合、 $n$ は必ず素数)

• $n\lt 2^{32}$ の範囲であれば、 $a$ は $\{2,7,61\}$ の3つで判定可能。 (Gerhard Jaeschke, 1993年)
• $a \ge n$ となる場合があるため、 $a\underset{\!\!\operatorname{mod}n\!\!}{\not\equiv}0$ の条件が追加されている点に注意。

Miller-Rabin 素数判定法(3) 実装

では、これを実装してみましょう

• PyPy3、C++17(gcc)、Rust それぞれでの実装例
• 63bit数1万個の素数判定の実行時間の例は以下の通り:
  • https://yukicoder.me/submissions/792430 : PyPy3: 426ms
  • https://yukicoder.me/submissions/794259 : C++17(gcc): 222ms
  • https://yukicoder.me/submissions/794257 : Rust: 208ms

# Python : nは64bitで表す事ができる自然数とし、1以下の整数もしくは合成数ならFalse、素数ならTrueを返す
def miller_rabin(n: int):
    if n == 2:
        return True
    if n < 2 or (n & 1) == 0:
        return False
    n1 = n - 1
    d = n1
    s = 0
    while (d & 1) == 0:
        d //= 2
        s += 1
    for a in [2,325,9375,28178,450775,9780504,1795265022]:
        if a % n == 0:
            continue
        t = pow(a, d, n) # = pow(a, d) % n
        if t == 1 or t == n1:
            continue
        for _ in range(s - 1):
            t = pow(t, 2, n) # = pow(t, 2) % n
            if t == n1:
                break # else節を実行せず次の for a in ... ループへ
        else: # breakでループを抜けず、普通に for _ in range() ... ループが終了した時
	    # for文に付くelse節の説明:
            # https://docs.python.org/ja/3/reference/compound_stmts.html#the-for-statement
            return False
    return True

// C++ : miller_rabin(n) : 1以下の整数もしくは合成数ならfalse、素数ならtrueを返す
#include <cstdbool>
#include <cstdint>
uint64_t modmul(uint64_t a, uint64_t b, uint64_t n) {
    return (uint64_t)(((__uint128_t)a) * ((__uint128_t)b) % ((__uint128_t)n));
}
uint64_t modpow(uint64_t a, uint64_t b, uint64_t n) {
    uint64_t t = ((b & 1) == 0) ? 1 : a;
    for (b >>= 1; b != 0; b >>= 1) {
        a = modmul(a, a, n);
	if ((b & 1) == 1) { t = modmul(t, a, n); } 
    }
    return t;
}
const uint64_t bases[] = {2,325,9375,28178,450775,9780504,1795265022};
bool miller_rabin(uint64_t n) {
    if (n == 2) { return true; }
    if (n < 2 || (n & 1) == 0) { return false; }
    uint64_t n1 = n - 1, d = n - 1;
    uint32_t s = 0;
    for (; (d & 1) == 0; d >>= 1) { s += 1; }
    for (const auto& base : bases) {
        uint64_t a = base;
        if (a >= n) {
	    a %= n;
	    if (a == 0) { continue; }
	}
        uint64_t t = modpow(a, d, n);
        if (t == 1) { continue; }
        for (uint32_t j = 1; t != n1; ++j) {
            if (j >= s) { return false; }
            t = modmul(t, t, n);
        }
    }
    return true;
}

// Rust : miller_rabin(n) : 1以下の整数もしくは合成数ならfalse、素数ならtrueを返す
pub fn modmul(a: u64, b: u64, n: u64) -> u64 {
    ((a as u128) * (b as u128) % (n as u128)) as u64
}
pub fn modpow(mut b: u64, mut p: u64, n: u64) -> u64 {
    let mut r = if (p & 1) == 0 { 1 } else { b };
    loop {
        p >>= 1;
	if p == 0 { return r; }
        b = modmul(b, b, n);
	if p & 1 != 0 { r = modmul(r, b, n); }
    }
}
pub fn miller_rabin(n: u64) -> bool {
    if n == 2 { return true; }
    if n < 2 || n & 1 == 0 { return false; }
    let n1 = n - 1;
    let s = n1.trailing_zeros();
    let d = n1 >> s;
    [2,325,9375,28178,450775,9780504,1795265022].iter().all(|&base| {
        let a = if base < n { base } else { base % n };
        if a == 0 { return true; }
        let mut t = modpow(a, d, n);
        if t == 1 || t == n1 { return true; }
        for _ in 1..s {
    	    t = modmul(t, t, n);
	    if t == n1 { return true; }
	}
        false
    })
}

もっとspeedupしてみよう

• 題意は既に満たしているが、 C++言語 や Rust言語 でさらに 10倍速 を目指す
• 128bit剰余算は遅いので多用を避けたい → モンゴメリ剰余乗算を使う
  • 128bit剰余算1回の代わりに、おおむね 64bit乗算2回+64bit加減算2回+上位64bit切り出し2回+下位64bit切り出し2回 を使う
  • Pythonではモンゴメリ剰余乗算での高速化は難しいのでこれは行わない
    • 例えばCPythonは整数型に$2^{30}$進数の内部表現を用いており、モンゴメリ剰余乗算を使うには実行コストが大きい
• C++ では128bit整数だけではなく、GCCのビルトイン関数(__builtin_uaddll_overflow, __builtin_usubll_overflow, __builtin_clzll, __builtin_ctzll)も使う
  • Rustでは整数型で overflowing_add, overflowing_sub, leading_zeros, trailing_zeros などが使える
• Miller-Rabin 素数判定法より速い方法は無いか → Baillie–PSW 素数判定法がある

乗法逆元

• ニュートン法 : $\operatorname{mod}2^{64}$ などの累乗数の法に対する乗法逆元の求め方の例
• 拡張ユークリッドの互除法 : 一般の法に対する乗法逆元の求め方の例

mod 2⁶⁴ における奇数 N の乗法逆元

与えられた奇数 $N$ に対し、整数 $N^{-1}$ を掛けて $2^{64}$ で割った余りが $1$ になる、整数 $N^{-1}$ の値を求める問題 $(NN^{-1}\underset{\!\!\operatorname{mod}2^{64}\!\!}\equiv 1)$ 実数の逆数を求めるニュートン法と同じ漸化式を用いる。

// C++
#include <cassert>
#include <cstdint>
uint64_t invmod_u64(uint64_t n) {
    assert(n % 2 == 1);
    uint64_t ni = n;
    for(int i = 0; i < 5; ++i) {
        ni = ni * (2 - n * ni);
    }
    assert(n * ni == 1);
    return ni;
}

// Rust
fn invmod_u64(n: u64) -> u64 {
    assert_eq!(n % 2, 1);
    let mut ni = n;
    for _ in 0..5 {
        ni = ni.wrapping_mul(
            2u64.wrapping_sub(n.wrapping_mul(ni))
        );
    }
    assert_eq!(n.wrapping_mul(ni), 1);
    ni
}

$$\begin{align*} &\textbf{algorithm | アルゴリズム:}\\ n'_0 &:= N \quad \begin{cases}N\text{ is positive odd number}\\ N\text{ は正の奇数}\end{cases} \\ n'_1 &:= n'_0\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}2^{64}\!\!}{-}\left(N\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}n'_0\right)\right) \\ n'_2 &:= n'_1\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}2^{64}\!\!}{-}\left(N\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}n'_1\right)\right) \\ n'_3 &:= n'_2\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}2^{64}\!\!}{-}\left(N\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}n'_2\right)\right) \\ n'_4 &:= n'_3\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}2^{64}\!\!}{-}\left(N\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}n'_3\right)\right) \\ n'_5 &:= n'_4\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}2^{64}\!\!}{-}\left(N\underset{\!\!\operatorname{mod}2^{64}\!\!}{\times}n'_4\right)\right) \\ &\text{return}\left(n'_5\right); \\ \end{align*}$$

$\textbf{proof | 証明:}$

• $x_j$ が、 $d\,x_j\mod m\equiv 1$ を満たす時、
  $x_{j+1}$ を $x_{j+1}=x_j(2-d\,x_j)$ と定義するなら、
  $dx_{j+1}\mod m^2\equiv 1$ を満たす。
  • $d\,x_j=1+km$ とおくと、
    $d\,x_{j+1}=d\,x_n(2-d\,x_n)=(1+km)(2-(1+km))$
    $=(1+km)(1-km)=(1-k^2m^2)\equiv 1\quad(\operatorname{mod}\ m^2)$
• 奇数$N$を $N=1+2k$ とおくと、 $\left(k(1+k)\right)$ は偶数で、 $\left(k(1+k)\right)/2$ は整数であるから
  • $n_0=N,\quad n_{j+1}=n_j(2-N\,n_j)$ とした時 $Nn_{j+1}=Nn_j(2-N\,n_j)$ より
  • $Nn_0=1+2^3\cdot \left(k(1+k)/2\right)\ \ \to\ \ Nn_0\mod 2^3\equiv 1$
  • $Nn_1=1-2^6\cdot \left(k(1+k)/2\right)^2\ \ \to\ \ Nn_1\mod 2^6\equiv 1$
  • $Nn_2=1-2^{12}\cdot \left(k(1+k)/2\right)^4\ \ \to\ \ Nn_2\mod 2^{12}\equiv 1$
  • $Nn_3=1-2^{24}\cdot \left(k(1+k)/2\right)^8\ \ \to\ \ Nn_3\mod 2^{24}\equiv 1$
  • $Nn_4=1-2^{48}\cdot \left(k(1+k)/2\right)^{16}\ \ \to\ \ Nn_4\mod 2^{48}\equiv 1$
  • $Nn_5=1-2^{96}\cdot \left(k(1+k)/2\right)^{32}\ \ \to\ \ Nn_5\mod 2^{96}\equiv 1$

mod 10⁹ における乗法逆元

$\mod 2^{64}$ の場合で触れた証明を踏まえると、$10$進数の場合でも同様に計算できる。有効桁数が1ループ毎に倍になる事の応用として、初期のループでは計算する桁数を抑え、徐々に桁数を増やして計算する例としても示す。

例題: $N=998244353, R=10^9$, $Nn'\mod R\equiv 1$ となる $n'$ を求める

$N$ と $R$ は互いに素 $(\gcd(N,R)=1)$ なので $\operatorname{mod}\ R$ における $N$ の乗法逆元 $n'$ は存在する

計算する桁数: 初期値:10進1桁 → 1ループ目:10進2桁 → 2ループ目:10進3桁 → 3ループ目:10進5桁 → 4ループ目:10進9桁

$$\begin{align*} &\textbf{algorithm | アルゴリズム:}\\ n'_0 &:= \begin{cases} 1\quad\text{ if }N\operatorname{mod}\ 10\equiv 1\\ 7\quad\text{ if }N\operatorname{mod}\ 10\equiv 3\\ 3\quad\text{ if }N\operatorname{mod}\ 10\equiv 7\\ 9\quad\text{ if }N\operatorname{mod}\ 10\equiv 9\\ \end{cases} \\ n'_1 &:= n'_0\underset{\!\!\operatorname{mod}10^2\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^2\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^2\!\!}{\times}n'_0\right)\right) \\ n'_2 &:= n'_1\underset{\!\!\operatorname{mod}10^3\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^3\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^3\!\!}{\times}n'_1\right)\right) \\ n'_3 &:= n'_2\underset{\!\!\operatorname{mod}10^5\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^5\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^5\!\!}{\times}n'_2\right)\right) \\ n'_4 &:= n'_3\underset{\!\!\operatorname{mod}10^9\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^9\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^9\!\!}{\times}n'_3\right)\right) \\ &\text{return}\left(n'_4\right); \\ \end{align*}$$

$\textbf{example | 例:}$

• $N=998244353, R=10^9$ として、 $Nn'\mod R\equiv 1$ を求める
• $N$ と $R$ は互いに素 $(\gcd(N,R)=1)$ なので乗法逆元は存在する
• $n'_0=7$
• $n'_1=n'_0\underset{\!\!\operatorname{mod}10^2\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^2\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^2\!\!}{\times}n'_0\right)\right)=7\underset{\!\!\operatorname{mod}10^2\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^2\!\!}{-}\left(53\underset{\!\!\operatorname{mod}10^2\!\!}{\times}7\right)\right)$
  $=7\underset{\!\!\operatorname{mod}10^2\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^2\!\!}{-}\left(371\mod 10^2\right)\right)=7\underset{\!\!\operatorname{mod}10^2\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^2\!\!}{-}71\right)$
  $=7\underset{\!\!\operatorname{mod}10^2\!\!}{\times}\left(-69\mod 10^2\right)=7\underset{\!\!\operatorname{mod}10^2\!\!}{\times}31=217\mod 10^2=17$
• $n'_2=n'_1\underset{\!\!\operatorname{mod}10^3\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^3\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^3\!\!}{\times}n'_1\right)\right)=17\underset{\!\!\operatorname{mod}10^3\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^3\!\!}{-}\left(353\underset{\!\!\operatorname{mod}10^3\!\!}{\times}17\right)\right)$
  $=17\underset{\!\!\operatorname{mod}10^3\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^3\!\!}{-}\left(6001\mod 10^3\right)\right)=17\underset{\!\!\operatorname{mod}10^3\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^3\!\!}{-}1\right)$
  $=17\underset{\!\!\operatorname{mod}10^3\!\!}{\times}1=17$
• $n'_3=n'_2\underset{\!\!\operatorname{mod}10^5\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^5\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^5\!\!}{\times}n'_2\right)\right)=17\underset{\!\!\operatorname{mod}10^5\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^5\!\!}{-}\left(44353\underset{\!\!\operatorname{mod}10^5\!\!}{\times}17\right)\right)$
  $=17\underset{\!\!\operatorname{mod}10^5\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^5\!\!}{-}\left(754001\mod 10^5\right)\right)=17\underset{\!\!\operatorname{mod}10^5\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^5\!\!}{-}54001\right)$
  $=17\underset{\!\!\operatorname{mod}10^5\!\!}{\times}\left(-53999\mod 10^5\right)=17\underset{\!\!\operatorname{mod}10^5\!\!}{\times}46001=782017\mod 10^5=82017$
• $n'_4=n'_3\underset{\!\!\operatorname{mod}10^9\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^9\!\!}{-}\left(N\underset{\!\!\operatorname{mod}10^9\!\!}{\times}n'_3\right)\right)=82017\underset{\!\!\operatorname{mod}10^9\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^9\!\!}{-}\left(998244353\underset{\!\!\operatorname{mod}10^9\!\!}{\times}82017\right)\right)$
  $=82017\underset{\!\!\operatorname{mod}10^9\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^9\!\!}{-}\left(81873007100001\mod 10^9\right)\right)=82017\underset{\!\!\operatorname{mod}10^9\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^9\!\!}{-}7100001\right)$
  $=82017\underset{\!\!\operatorname{mod}10^9\!\!}{\times}\left(2\underset{\!\!\operatorname{mod}10^9\!\!}{-}7100001\right)=82017\underset{\!\!\operatorname{mod}10^9\!\!}{\times}\left(-7099999\mod 10^9\right)$
  $=82017\underset{\!\!\operatorname{mod}10^9\!\!}{\times}992900001=81434679382017\mod 10^9=679382017$
• 検算: $Nn'\mod R=998244353\underset{\!\!\operatorname{mod}10^9\!\!}{\times}679382017=678189262000000001\mod 10^9=1$

mod R での乗法逆元を使った倍数判定

$\operatorname{mod}\ R$ における $D$ の乗法逆元を $d'$ として整数 $N$ に対し $(0\lt D\lt R,\ 0\le N\lt R,\ \gcd(R,D)=1,\ Dd'\mod R\equiv 1)$、

• $N$ が $D$ の倍数なら、 $(N/D)\mod R=(Nd')\mod R\le\lfloor(R-1)/D\rfloor$ である。
• $N$ が $D$ の倍数でないなら、 $(Nd')\mod R\gt\lfloor(R-1)/D\rfloor$ である。

#include <cstdbool>
#include <cstdint>
bool is_multiple_of_3_divmod(uint64_t n) { return n % 3 == 0; }
bool is_multiple_of_3_mulinvmod(uint64_t n) { // C++: n が 3 の倍数なら true、上の関数と同等
    return (uint64_t)(n * 0xAAAAAAAAAAAAAAABULL) <= (uint64_t)0x5555555555555555ULL; }

fn is_multiple_of_3_divmod(n: u64) -> bool { n % 3 == 0 }
fn is_multiple_of_3_mulinvmod(n: u64) -> bool { // Rust: n が 3 の倍数なら true、上の関数と同等
    n.wrapping_mul(0xAAAA_AAAA_AAAA_AAABu64) <= 0x5555_5555_5555_5555u64 }

いずれも同じアセンブリ出力になり、コンパイラも実は同様の最適化をしている事が分かる。

C++: https://godbolt.org/z/M8z1W1hqM
Rust: https://rust.godbolt.org/z/cvzKrvYTo

この例では除数が定数なのでコンパイル時の最適化が効いているが、変数の場合は必ずしもこのような最適化が行われるとは限らない。

拡張ユークリッドの互除法

$0$ でない整数 $a,b$ に対し、 $ax+by=\operatorname{gcd}(a,b)$ (ベズーの等式(Bézout's identity)) を満たす 整数 $x,y$ : $\operatorname{abs}(x)\le\operatorname{abs}(b/\operatorname{gcd}(a,b))$, $\operatorname{abs}(y)\le\operatorname{abs}(a/\operatorname{gcd}(a,b))$, $\max(\operatorname{abs}(x),\operatorname{abs}(y))\gt 0$ が存在する。特に、 $\operatorname{abs}(a) \ge 2, \operatorname{abs}(b) \ge 2$ かつ $a$ と $b$ が互いに素($a$ と $b$ の最大公約数が $1$)である時、 $ax+by=1$ を満たす整数 $x,y$ : $\operatorname{abs}(x)\lt \operatorname{abs}(b), \operatorname{abs}(y)\lt \operatorname{abs}(a)$ が存在する。ここで、$\operatorname{abs}(x)$ は絶対値関数、$\operatorname{gcd}(a,b)$ は最大公約数関数。

// Rust
pub fn extgcd_i64(a: std::num::NonZeroI64, b: std::num::NonZeroI64) -> (i64, i64, i64) {
    // a*x + b*y = d の解の一つを出力(x, y, d = gcd(a,b))
    let (mut wa, mut wb, mut x, mut y, mut nx, mut ny) = (a.get(), b.get(), 1, 0, 0, 1);
    while wb != 0 {
        let (q, r) = (wa / wb, wa % wb);
        let (tx, ty) = (x - (q * nx), y - (q * ny));
        wa = wb; wb = r; x = nx; y = ny; nx = tx; ny = ty;
    }
    assert_eq!(a.get().wrapping_mul(x).wrapping_add(b.get().wrapping_mul(y)), wa);
    (x, y, wa) // x, y, gcd(a,b)
}

$ax+by=\operatorname{gcd}(a,b)$ (ベズーの等式(Bézout's identity)) を変形して負数を取り扱わず計算する
$ax-by=\operatorname{gcd}(a,b)$ として $0\lt a, 0\lt b, 0\le x\le b/\gcd(a,b), 0\le y \le a/\gcd(a,b)$ の制約条件を付け解を求める

// Rust
pub fn extgcd_u64(a: std::num::NonZeroU64, b: std::num::NonZeroU64) -> (u64, u64, u64) {
    // a*x - b*y = d の解の一つを出力(x, y, d = gcd(a,b))
    let (mut wa, mut wb, mut x, mut y, mut nx, mut ny, mut f) = (a.get(), b.get(), 1, 0, 0, 1, false);
    while wb != 0 {
        let (q, r) = (wa / wb, wa % wb);
        let (tx, ty) = (x + (q * nx), y + (q * ny));
        wa = wb; wb = r; x = nx; y = ny; nx = tx; ny = ty; f = !f;
    }
    assert_eq!(nx * wa, b.get()); // nx * gcd(a,b) == b
    assert_eq!(ny * wa, a.get()); // ny * gcd(a,b) == a
    if f { x = nx - x; y = ny - y; }
    assert!(x <= nx); // x <= b / gcd(a,b)
    assert!(y <= ny); // y <= a / gcd(a,b)
    assert_eq!(a.get().wrapping_mul(x).wrapping_sub(b.get().wrapping_mul(y)), wa); // a*x-b*y == gcd(a,b)
    (x, y, wa) // x, y, gcd(a,b)
}

平方数の判定

$\operatorname{mod}8$ を取った時、平方数が取り得る値は $0,1,4$ しかない。逆に、
$\operatorname{mod}8$ を取った値が $2,3,5,6,7$ である場合、その数は平方数ではない事が分かる。
特に、整数 $k$ に対して $k(k+1)$ は偶数であるから、 $(2k+1)^2\equiv 4k(k+1)+1\equiv 1\ (\operatorname{mod}8)$ より、 奇数の平方数に $\operatorname{mod}8$ を取ると $1$ となる。

$$\begin{matrix} 0^2 \operatorname{mod}8\equiv 0 & 1^2 \operatorname{mod}8\equiv 1 \\ 2^2 \operatorname{mod}8\equiv 4 & 3^2 \operatorname{mod}8\equiv 1 \\ 4^2 \operatorname{mod}8\equiv 0 & 5^2 \operatorname{mod}8\equiv 1 \\ 6^2 \operatorname{mod}8\equiv 4 & 7^2 \operatorname{mod}8\equiv 1 \\ \end{matrix}$$

$\operatorname{mod}32$ を取った場合は 32個中7個(約22%)が平方数かもしれない値、
$\operatorname{mod}4095$ を取った場合は 4095個中336個(約8.2%)が平方数かもしれない値となり、
それでも平方数かもしれない数について、最後にニュートン法で平方根を計算している。

// 擬平方数判定 mod32 7/32 0.2188 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_mod32(x: u64) -> bool {
    (0x02030213u32 >> ((x as u32) & 31)) & 1 == 1
}
// 擬平方数判定 mod4095 336/4095 0.0821 (true: 平方数かもしれない, false: 平方数ではない)
fn issq_mod4095(x: u64) -> bool {
    const SQTABLE_MOD4095: [u64; 64] = [
        0x2001002010213,0x4200001008028001,0x20000010004,0x80200082010,
	0x1800008200044029,0x120080000010,0x2200000080410400,0x8100041000200800,
	0x800004000020100,0x402000400082201,0x9004000040,0x800002000880,
	0x18002000012000,0x801208,0x26100000804010,0x80000080000002,
	0x108040040101045,0x20c00004000102,0x400000100c0010,0x1300000040208,
        0x804000020010000,0x1008402002400080,0x201001000200040,0x4402000000806000,
	0x10402000000,0x1040008001200801,0x4080000000020400,0x10083080000002,
	0x8220140000040000,0x800084020100000,0x80010400010000,0x1200020108008060,
	0x180000000,0x400002400000018,0x4241000200,0x100800000000,
	0x10201008400483,0xc008000208201000,0x800420000100,0x2010002000410,
        0x28041000000,0x4010080000024,0x400480010010080,0x200040028000008,
	0x100810084020,0x20c0401000080000,0x1000240000220000,0x4000020800,
	0x410000000480000,0x8004008000804201,0x806020000104000,0x2080002000211000,
	0x1001008001000,0x20000010024000,0x480200002040000,0x48200044008000,
	0x100000000010080,0x80090400042,0x41040200800200,0x4000020100110,
        0x2000400082200010,0x1008200000000040,0x2004800002,0x2002010000080
    ];
    let p = (x % 4095) as usize;
    (SQTABLE_MOD4095[p >> 6] >> (p & 63)) & 1 == 1
}
// ニュートン法による整数平方根
fn sqrt_newton(x: u64) -> u64 {
    if x <= 1 { return x; }
    let k = 32 - ((x - 1).leading_zeros() >> 1);
    let mut s = (1 as u64) << k; // s = 2**k
    let mut t = (s + (x >> k)) >> 1; // t = (s + x/s)/2
    // whileループ回数=除算回数は u64:最大6回 u32:最大5回 を想定
    // s > floor(sqrt(x)) -> floor(sqrt(x)) <= t < s
    // s == floor(sqrt(x)) -> s == floor(sqrt(x)) <= t <= floor(sqrt(x)) + 1
    while t < s { s = t; t = (s + (x / s)) >> 1; }
    s
}
// 平方数判定 (true: 平方数である, false: 平方数ではない)
fn issq_newton(x: u64) -> bool { let sqrt = sqrt_newton(x); sqrt * sqrt == x }
// 平方数判定 (true: 平方数である, false: 平方数ではない)
fn issq(x: u64) -> bool {
    issq_mod32(x) && // 擬平方数判定 mod32 7/32 0.2188
    issq_mod4095(x) && // 擬平方数判定 mod4095 336/4095 0.0821
    issq_newton(x) // ニュートン法による整数平方根
}

(参考) 整数上でのニュートン法にて平方根が計算できる理由

• ニュートン・ラフソン法での平方根の近似値の初期値 $s_0$ が $s_0 \ge \lfloor\sqrt a\rfloor$ である場合を想定し（ $s_n \lt \lfloor\sqrt a\rfloor$ とした場合の挙動をここでは考慮していない事に注意 ）、 $s_{n+1}=\lfloor(s_n+\lfloor a/s_n\rfloor)/2\rfloor$ とする。 $s_n\in\mathbb{N},\ a\in\mathbb{N},\ s_n\ge 1,\ a\ge 1$ であるとする。
• (補題) $s_n$ は自然数であるから、
  $\displaystyle{s_{n+1}=\left\lfloor\left(s_n+\left\lfloor\frac{a}{s_n}\right\rfloor\right)/2\right\rfloor=\left\lfloor\left\lfloor s_n+\frac{a}{s_n}\right\rfloor/2\right\rfloor=\left\lfloor\left(s_n+\frac{a}{s_n}\right)/2\right\rfloor=\left\lfloor\frac{s_n^2+a}{2s_n}\right\rfloor}$
• (1) $s_n\gt\lfloor\sqrt{a}\rfloor$ のとき、 $\lfloor\sqrt{a}\rfloor\le s_{n+1}\lt s_n$ である。(この条件の間は単調減少する保証)
  • $s_n\gt\lfloor\sqrt{a}\rfloor$ より $s_n\gt\sqrt{a}$ であり、 $\varepsilon\gt 0$ を使って $s_n=(1+\varepsilon)\sqrt{a}$ とすると、
    
    $$\lfloor\sqrt{a}\rfloor\le\left\lfloor\left(1+\frac{\varepsilon^2}{2(1+\varepsilon)}\right)\sqrt{a}\right\rfloor=\left\lfloor\frac{2+2\varepsilon+\varepsilon^2}{2(1+\varepsilon)}\sqrt{a}\right\rfloor=\left\lfloor\frac{s_n^2+a}{2s_n}\right\rfloor=s_{n+1}\\\ \\s_{n+1}\le\frac{s_n^2+a}{2s_n}\lt\frac{s_n^2+s_n^2}{2s_n}=s_n$$
• (2) $s_n=\lfloor\sqrt{a}\rfloor$ のとき、 $\lfloor\sqrt{a}\rfloor\le s_{n+1}\le\lfloor\sqrt{a}\rfloor+1$ である。($\lfloor\sqrt{a}\rfloor$ 到達時の挙動の保証)
  • $s_n=\lfloor\sqrt{a}\rfloor,\quad\sqrt{a}-1\lt s_n\le\sqrt{a}$ より、 $s_n^2\le a\lt(s_n+1)^2$
  • $s_n$ は自然数であるから $\lfloor s_n\rfloor=s_n$ および $\displaystyle{\frac{1}{2s_n}\lt 1}$ より、
    
    $$\lfloor\sqrt{a}\rfloor=\left\lfloor\frac{s_n^2+s_n^2}{2s_n}\right\rfloor\le s_{n+1}\le\left\lfloor\frac{s_n^2+(s_n+1)^2}{2s_n}\right\rfloor=\left\lfloor s_n+1+\frac{1}{2s_n}\right\rfloor=\lfloor\sqrt{a}\rfloor+1$$

$a\lt 2^{64}$, 初期値$s_0$ が $\lfloor\sqrt{a}\rfloor\le s_0\le 2.9\sqrt{a}$ 程度の範囲であれば、 $s_{n+1}=\lfloor(s_n+\lfloor a/s_n\rfloor)/2\rfloor$ の操作 7回目までに $s_n=\lfloor\sqrt{a}\rfloor$ となる $0\le n\lt 7$ を見つける事ができる。

• $a\le 2^{4},\quad n\ge 1,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{10},\quad n\ge 2,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{23},\quad n\ge 3,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{32},\quad n\ge 4,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2.837\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{48},\quad n\ge 4,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{64},\quad n\ge 5,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2.916\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{99},\quad n\ge 5,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{128},\quad n\ge 6,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2.957\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$
• $a\le 2^{200},\quad n\ge 6,\quad \lfloor\sqrt{a}\rfloor\le s_0\le 2\sqrt{a} \quad\Rightarrow\quad \lfloor\sqrt{a}\rfloor\le s_n\le\lfloor\sqrt{a}\rfloor+1$

モンゴメリ剰余乗算

互いに素(最大公約数が $1$ )な自然数 $N, R\quad(N\lt R)$ を考える。
ここでは、 $N:\text{正の奇数}, R=2^{64}$ として考える。
そして、加減乗算・比較演算をする際、それぞれの変数にあらかじめ $\times R \mod N$ という変換(以下、モンゴメリ変換・モンゴメリ表現)をして、以下のように処理を置き換えていく。

$$\begin{align*} a+b\mod N=c\mod N & \quad\to\quad & (a\times R\mod N)+(b\times R\mod N)\mod N=c\times R\mod N \\ a-b\mod N= d\mod N & \quad\to\quad & (a\times R\mod N)-(b\times R\mod N)\mod N=d\times R\mod N \\ a\times b\mod N=e\mod N & \quad\to\quad & \{(a\times R\mod N)\times (b\times R\mod N)\}\times R^{-1}\mod N=e\times R\mod N \\ f\mod N=g\mod N & \quad\to\quad & f\times R\mod N=g\times R\mod N \\ \end{align*}$$

$R$ の値での除算・剰余算が効率的に行える場合、この モンゴメリ剰余乗算 $\{(a\times R\mod N)\times (b\times R \mod N)\}\times R^{-1}\mod N=e\times R\mod N$
の左辺が実は効率的に計算可能(モンゴメリリダクション: $\operatorname{MR}(T)=TR^{-1}\mod N$, $0\le T\lt NR$ )というのが肝になっている。

モンゴメリリダクション

$N\lt R,\ \gcd(N,R)=1$ であり、$N^{-1}$ は $NN^{-1}\underset{\!\!\operatorname{mod}R\!\!}{\equiv}1$ 、 $R^{-1}$ は $RR^{-1}\underset{\!\!\operatorname{mod}N\!\!}\equiv 1$ 、 $T$ は $0\le T\lt RN$ を満たすとする。
モンゴメリリダクション $\operatorname{MR}(T) = TR^{-1}\mod N$ は以下の通り。$\lfloor{x}\rfloor$ は床関数(実数に対してそれ以下の最大の整数を返す関数)。

$t \leftarrow \lfloor{T/R}\rfloor - \lfloor{(TN^{-1} \mod R)N/R}\rfloor;$
$\text{if }(t\lt 0)\text{ then \{ return }(t+N)\text{ \} else \{ return }(t)\text{ \}}$

証明:

• $T - (TN^{-1}\mod R)N\underset{\!\!\operatorname{mod}R\!\!}\equiv T+TN^{-1}N\underset{\!\!\operatorname{mod}R\!\!}\equiv T-T\underset{\!\!\operatorname{mod}R\!\!}\equiv 0$ より、 $(T - (TN^{-1}\mod R)N)$ は $R$ で割り切れる。
• $tR\underset{\!\!\operatorname{mod}N\!\!}\equiv T-(TN^{-1}\mod R)N\underset{\!\!\operatorname{mod}N\!\!}\equiv T$ より、 $t\underset{\!\!\operatorname{mod}N\!\!}\equiv TR^{-1}$ である。
• $T\lt RN,\ (TN^{-1}\mod R)N\lt RN$ より、 $-N\lt\ t\ =\ \lfloor{T/R}\rfloor - \lfloor{(TN^{-1} \mod R)N/R}\rfloor\ \lt N$ であるため、手続きが返す値は $N$ より小さい。

モンゴメリ変換・剰余乗算

あらかじめ $R_2=R^2\mod N$ を用意しておき、 $0\le a\lt N$ と $0\le b\lt N$ の乗算剰余 $c=ab\mod N$ を以下のようにして求める。$\operatorname{MR}(T) = TR^{-1}\mod N\quad(RR^{-1}\underset{\!\!\operatorname{mod}N\!\!}\equiv 1)$

$A\leftarrow\operatorname{MR}(aR_2)\ \{=aR\mod N\},\quad B\leftarrow\operatorname{MR}(bR_2)\ \{=bR\mod N\}$,$C\leftarrow\operatorname{MR}(AB)\ \{=abR\mod N\},\quad c\leftarrow\operatorname{MR}(C)\ \{=ab\mod N\}.$

ただし、$C$ の値をまた直ぐに次の演算に用いる場合、 $c\leftarrow\operatorname{MR}(C)$ の処理は不要となる。
乗算剰余を1回だけ求め、すぐにモンゴメリ表現から元に戻す場合は以下のように計算できる。

$c\leftarrow\operatorname{MR}(\operatorname{MR}(ab)R_2)=ab\mod N$

べき剰余 $a^k\mod N$ を求めたい場合、まず $A\leftarrow\operatorname{MR}(aR_2)=aR\mod N$ を計算し、 べき乗の計算は $C\leftarrow\operatorname{MR}(AB)=abR\mod N$ を繰り返し適用する。べき剰余の計算においても、通常のべき乗の計算と同様、バイナリ法($a^1,a^2,a^4,a^8,a^{16}\dots$を計算し、$k$ の$2$進数表現を基にそれらから必要なものを掛け合わせる)などの効率化手法が利用できる。

10進数でモンゴメリ剰余乗算

$N=998244353, R=10^9$ とする。 $N$ と $R$ は互いに素であり、 $N\lt R$ を満たし、 $R_2=R^2\mod N=716070898$ 、 $NN^{-1}\underset{\!\!\operatorname{mod}R\!\!}\equiv 1$ を満たす値は $N^{-1}=679382017$ である $(998244353\times 679382017=678189262000000001\underset{\!\!\operatorname{mod}R\!\!}\equiv 1)$。
$a=904894094, b=560163165$ とし、 $c\leftarrow\operatorname{MR}(\operatorname{MR}(ab)R_2)\quad\{=ab\mod N\}$ の手法で $ab\mod N$ を計算してみる。

• $N=998244353, R=10^9, R_2=716070898, N^{-1}=679382017$
• $t \leftarrow \lfloor{T/R}\rfloor - \lfloor{(TN^{-1} \mod R)N/R}\rfloor; \quad \text{if }(t\lt 0)\text{ then \{ return }(t+N)\text{ \} else \{ return }(t)\text{ \}}$
• $c\leftarrow\operatorname{MR}(\operatorname{MR}(ab)R_2)\quad\{=ab\mod N\}$
• $\operatorname{MR}(ab)=\operatorname{MR}(904894094 \times 560163165) = \operatorname{MR}(506888339684847510)$
  $=(\lfloor 506888339684847510/R\rfloor - \lfloor(684847510\times 679382017\mod R)\times 998244353/R\rfloor)\mod N$
  $=(\lfloor 506888339684847510/R\rfloor - \lfloor(465273082681227670\mod R)\times 998244353/R\rfloor)\mod N$
  $=(\lfloor 506888339684847510/R\rfloor - \lfloor 681227670\times 998244353/R\rfloor)\mod N$
  $=(\lfloor 506888339684847510/R\rfloor - \lfloor 680031674684847510/R\rfloor)\mod N$
  $=(506888339 - 680031674)\mod N=(-173143335)\mod N=N-173143335=825101018$
• $\operatorname{MR}(\operatorname{MR}(ab)R_2) = \operatorname{MR}(825101018\times 716070898) = \operatorname{MR}(590830826899974164)$
  $=(\lfloor 590830826899974164/R\rfloor - \lfloor(899974164\times 679382017\mod R)\times 998244353/R\rfloor)\mod N$
  $=(\lfloor 590830826899974164/R\rfloor - \lfloor(611426262786208788\mod R)\times 998244353/R\rfloor)\mod N$
  $=(\lfloor 590830826899974164/R\rfloor - \lfloor 786208788\times 998244353/R\rfloor)\mod N$
  $=(\lfloor 590830826899974164/R\rfloor - \lfloor 784828482899974164/R\rfloor)\mod N$
  $=(590830826 - 784828482)\mod N=(-193997656)\mod N=N-193997656=804246697$
• よって、 $ab\mod N=804246697$

モンゴメリ剰余乗算 まとめ

• 互いに素(最大公約数が $1$ )な自然数 $N, R\quad(N\lt R)$ を考え、かつ $N:\text{正の奇数}, R=2^{64}$ のように、 $R$ による除算・剰余算が容易な値を定める。
• モンゴメリ表現 (モンゴメリ変換 $\times R \mod N$ をした値) を用いると、加減乗算と $R$ による除算・剰余算 によって $\mod N$ における剰余乗算を計算できる。
• モンゴメリ剰余乗算は、同じ法 $N$ を繰り返し使う場合に有用。$N$ の変更が多いと、 その分 $R_2=R^2\mod N$ や $N^{-1}=1/N\mod R$ などの定数を求め直すコストが重くなる。
• モンゴメリ剰余乗算を用いて演算効率を上げるためには、モンゴメリ表現への変換（モンゴメリ変換）・逆変換（乗算に伴わないモンゴメリリダクション）は、特にループ処理の内側では極力行わない。モンゴメリ表現のままで演算や比較ができるよう、必要な定数はループのなるべく外側であらかじめ求めておく。
  • モンゴメリ表現された変数値とモンゴメリ表現されていない定数値の比較がループ内で必要な時も、あらかじめモンゴメリ変換した定数値を準備しておくと、ループ内でモンゴメリ表現された変数値をモンゴメリリダクションで逆変換する必要が無くなる
  • 例: $\operatorname{MR}(0)=0,\quad\operatorname{MR}(R\mod N)=1,\quad\operatorname{MR}(-R\mod N)=N-1$

Baillie-PSW素数判定法

https://en.wikipedia.org/wiki/Baillie–PSW_primality_test

Baillie-PSW素数判定法は3以上の奇数 $n$ に対して以下のように行う。

1. base値 を $2$ として Miller-Rabin素数判定を行う。失敗すれば composite(合成数) と出力して終了する。
   • 1.のみでは、以下の合成数は擬素数となる。(OEIS:A001262 Strong pseudoprimes to base 2.) $2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,\dots$
2. $D={5,-7,9,-11,\dots}$ の中で最初に $\left({D}|{n}\right)=-1$ となる $D$ を見つける。 $\left({a}|{n}\right)$ はヤコビ記号。 $n$ が平方数の場合 $D$ は永遠に見つからないため、適当なステップ数の経過時に平方数か確認し、平方数なら composite(合成数) と出力して終了する。
3. $P=1,Q=(1-D)/4$ とし、 $\gcd(n,Q)=1$ である事を確認する。 $U_k(P,Q),V_k(P,Q)$ のリュカ数列を用いて 強いリュカ擬素数(strong Lucas pseudoprime)テスト を行う。失敗すれば composite(合成数) と出力して終了する。
   • 2.3.のみでは、以下の合成数は擬素数となる。(OEIS:A217255 Strong Lucas pseudoprimes.) $5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,\dots$
4. 1.2.3.の何れでも終了しなかった場合、 probably prime(おそらく素数) と出力して終了する。これは少なくとも $n\lt 2^{64}$ の範囲では正しく素数と合成数を判別可能である。この判定法をすり抜けて擬素数となる合成数はおそらく無数に存在すると推測されるが、今のところその具体例は発見されていない。

クロネッカー記号・ヤコビ記号・ルジャンドル記号

https://ja.wikipedia.org/wiki/平方剰余
https://ja.wikipedia.org/wiki/ヤコビ記号

任意の整数 $a$ と $n$ に対し、クロネッカー記号(Kronecker symbol) $\left({a}|{n}\right)$ は $n$ の奇素因数 $p$ に対応する ルジャンドル記号(Legendre symbol) の積と、追加因数 $(a|0),(a|{\pm 1}),(a|2)$ によって定義される。 ヤコビ記号(Jacobi symbol) はそのうち、 $n$ が正の奇数の範囲における定義である。 $n=u\,p_1^{\alpha_1}p_2^{\alpha_2}\cdots p_k^{\alpha_k}$ と $n$ が素因数分解 $(u={\pm 1})$ されるとして

$$\left(a|n\right)=\left(a|u\right)\left(a|p_1\right)^{\alpha_1}\left(a|p_2\right)^{\alpha_2}\cdots\left(a|p_k\right)^{\alpha_k}$$

ルジャンドル記号 $\displaystyle{\left(a|p\right)}$ は 任意の整数 $a$ と全ての奇素数 $p$ に対して次のように定義される。

$$\left({a}|{p}\right)=\begin{cases} 0 & \text{: } a\equiv 0\ \ (\operatorname{mod} p), \\ 1 & \text{: } a\not\equiv 0\ \ (\operatorname{mod} p)\text{ and for some integer } x \text{ (平方剰余)}:\ a\equiv x^2\ \ (\operatorname{mod} p), \\ -1 & \text{: } a\not\equiv 0 \ \ (\operatorname{mod} p)\text{ and there is no such } x \text{ (平方非剰余)}. \\ \end{cases}$$

クロネッカー記号では、追加で以下の定義を行う。ヤコビ記号ではこのうち $(a|1)$ のみ使われる。

• $\left(a|1\right)=1$
• $\left(a|{-1}\right)=\begin{cases}1&\text{ if }a\ge 0,\\ -1&\text{ if }a\lt 0.\end{cases}$
• $\left(a|2\right)=\begin{cases}0&\text{ if }a\text{ is even},\\ 1&\text{ if }a\equiv 1,7\mod 8,\\ -1&\text{ if }a\equiv 3,5\mod 8.\end{cases}$
• $\left(a|0\right)=\begin{cases}1&\text{ if }a=\pm 1,\\ 0&\text{otherwise.}\end{cases}$

ヤコビ記号の性質

ヤコビ記号 ($n$ が正の奇数) の場合の性質を示す。必ずしもクロネッカー記号 ($n$ が整数) には当てはまらない。

• $a\equiv b\ \ (\operatorname{mod}\ n)$ の時、 $(a|n)=(b|n)$
• $(a|n)=\begin{cases}0&\text{if }\gcd(a,n)\ne 1\\\pm 1&\text{if }\gcd(a,n)=1\end{cases}$
• $(ab|n)=(a|n)(b|n)$
• $(-a|n)=(-1|n)(a|n)=(-1)^{\frac{n-1}{2}}(a|n)=\begin{cases}(a|n)&\text{if }n\not\equiv 3\ \ (\operatorname{mod}4)\\-(a|n)&\text{if }n\equiv 3\ \ (\operatorname{mod}4)\end{cases}$
• $(2a|n)=(2|n)(a|n)=(-1)^{\frac{n^2-1}{8}}(a|n)=\begin{cases}(a|n)&\text{if }n\equiv 1,7\ \ (\operatorname{mod}8)\\-(a|n)&\text{if }n\equiv 3,5\ \ (\operatorname{mod}8)\end{cases}$
• $m$ と $n$ が奇数で互いに素な正整数である場合、 $(m|n)(n|m)=(-1)^{\frac{m-1}{2}\frac{n-1}{2}}=\begin{cases}1&\text{if }n\equiv 1\ \ (\operatorname{mod}4)\text{ or }m\equiv 1\ \ (\operatorname{mod}4)\\-1&\text{if }n\equiv m\equiv 3\ \ (\operatorname{mod}4)\end{cases}$
• $(a|n)=-1$ の場合、 $a$ は $n$ を法として平方非剰余である。
• $a$ が $n$ を法として平方剰余であり、 $\gcd(a,n)=1$ の場合、 $(a|n)=1$ である。
• $(a|n)=1$ の場合、 $a$ は $n$ を法として平方剰余であるかもしれないし、ないかもしれない。奇素数 $p$ において $(a|p)=1$ の場合、 $a$ は $p$ を法として平方剰余である。

ヤコビ記号の計算

ユークリッドの互除法に似た方法を使い、 $\mathcal{O}(\log a\log b)$ 程度の処理回数で計算できる。

// C++ : ヤコビ記号 (a/n)
#include <cassert>
int jacobi(int a, int n) {
    assert(n > 0 && (n & 1) == 1); // n は正の奇数に限る
    // step 1 : n を法として a を減らす (ユークリッド除法)
    int t = 1; a = a % n; if (a < 0) { a = a + n; }
    // step 3 : a が 1 のとき、結果は 1。 a と n が互いに素でなければ、結果は 0。
    while (a != 0) {
        // step 2 : 偶数の a を抽出する
        while ((a & 1) == 0) { a = a >> 1; if ((n & 7) == 3 || (n & 7) == 5) { t = -t; } }
        // step 4 : 記号を反転し、 n を法として a を減らす
        if ((a & n & 3) == 3) { t = -t; }
        int r = n; n = a; a = r; a = a % n;
    }
    if (n == 1) { return t; } else { return 0; }
}

リュカ擬素数テスト

https://en.wikipedia.org/wiki/Lucas_pseudoprime

自然数 $P$ と 整数 $Q$ に対して $D=P^2-4Q$ とおき、 $U_k(P,Q),V_k(P,Q)$ をそれぞれ対応するリュカ数列とする。

$U_0(P,Q)=0, U_1(P,Q)=1, U_k(P,Q)=PU_{k-1}(P,Q)-QU_{k-2}(P,Q)\text{ for }k\gt 1$
$V_0(P,Q)=2, V_1(P,Q)=P, V_k(P,Q)=PV_{k-1}(P,Q)-QV_{k-2}(P,Q)\text{ for }k\gt 1$

$n$ を $\gcd(n,Q)=1$ である正の奇数とし、 $(D|n)$ をヤコビ記号とする。また、 $\delta(n)=n-(D|n)=d\cdot 2^s$ とおく。( $d$ は奇数)

（リュカ擬素数テスト）もし $n$ が素数であれば、 $U_{\delta(n)}=0\mod n$ が成り立つ。また、成り立たない場合 $n$ は素数ではない。

（強いリュカ擬素数テスト）もし $n$ が素数かつ $\gcd(n,D)=1$ である時、ある $0\le r\lt s$ において、 $U_d\equiv 0\mod n$ または $V_{d\cdot 2^r}\equiv 0\mod n$ のいずれかを満たす。

リュカ数列

https://ja.wikipedia.org/wiki/リュカ数列

$U_0(P,Q)=0, U_1(P,Q)=1, U_k(P,Q)=PU_{k-1}(P,Q)-QU_{k-2}(P,Q)\text{ for }k\gt 1$
$V_0(P,Q)=2, V_1(P,Q)=P, V_k(P,Q)=PV_{k-1}(P,Q)-QV_{k-2}(P,Q)\text{ for }k\gt 1$

$D=\{5,-7,9,-11,\cdots\},P=1,Q=(1-D)/4$ のリュカ数列は以下のように $\mathcal{O}(\log_2 n)$ で $(n+1)$ の$2$進数表現を最上位bitから見ながら $U_{n+1}\!\!\mod n,V_{n+1}\!\!\mod n$ を素早く計算する漸化式がある。

• $U_1=1,V_1=P=1$
• $U_{2k}\mod n=U_k\cdot V_k\mod n$
• $V_{2k}\mod n=V_k^2-2Q^k\mod n=(V_k^2+D\cdot U_k^2)/2\mod n=(V_k^2+D\cdot U_k^2+n)/2\mod n$
• $U_{2k+1}\mod n=(P\cdot U_{2k}+V_{2k})/2\mod n=(P\cdot U_{2k}+V_{2k}+n)/2\mod n$
• $V_{2k+1}\mod n=(D\cdot U_{2k}+P\cdot V_{2k})/2\mod n=(D\cdot U_{2k}+P\cdot V_{2k}+n)/2\mod n$

$/2\mod n$ の分子が奇数で割り切れない場合、右辺のように $n$ を足して偶数にしてから除算する。

モンゴメリ剰余乗算・BPSW素数判定法の実装結果

• モンゴメリ剰余乗算 導入前 (Miller-Rabin素数判定法+128bit剰余乗算)
  • https://yukicoder.me/submissions/792430 : PyPy3: 426ms
  • https://yukicoder.me/submissions/794259 : C++17(gcc): 222ms
  • https://yukicoder.me/submissions/794257 : Rust: 208ms
• Miller-Rabin素数判定法+モンゴメリ剰余乗算 導入後
  • https://yukicoder.me/submissions/793551 : C++17(gcc): 31ms
  • https://yukicoder.me/submissions/793560 : Rust: 27ms
• BPSW素数判定法+モンゴメリ剰余乗算 導入後
  • https://yukicoder.me/submissions/793997 : C++17(gcc): 22ms
  • https://yukicoder.me/submissions/794215 : Rust: 13ms

モンゴメリ剰余乗算実装前の10倍速をこれで実現できた。C言語・同様のアルゴリズムで 14ms を計測した提出例も存在したので、 C++言語でももう少し最適化の余地はあったかもしれない。