モジュラ逆数の計算法
モジュラ逆数の計算
(サンプルコードでは主に
モジュラ逆数
ここでは、3つの方法でのモジュラ逆数の求め方を紹介します。
- オイラーの定理を用いる方法 (
が素数である場合、m の素因数分解が既知である場合 )m - 変則的な拡張ユークリッド互除法を用いる方法 ( 一般の自然数
に対する場合 )m - ニュートン・ラフソン法を用いる方法 (
が累乗数である場合 )m
m が素数、もしくは m の素因数分解が既知である場合のモジュラ逆数 (オイラーの定理)
オイラーのトーシェント関数
\gcd(n,m)=1\quad\Longrightarrow\quad n^{\varphi(m)}\equiv 1\ (\operatorname{mod}m) \gcd(n,m)=1\quad\Longrightarrow\quad n^{-1}\equiv n^{\varphi(m)-1}\ (\operatorname{mod}m)
オイラーのトーシェント関数
オイラーのトーシェント関数
-
が素数であれば、m \varphi(m)=m-1 -
が素数ではなく、m の素因数分解が以下のように (m は互いに異なる素数 ) 既知であれば、p_1,p_2,p_3,\dots,p_d
サンプルコード (Rust)
// モジュラ逆数 $n^{-1} \mod 998244353$ をオイラーの定理から計算 (998244353 は素数)
pub fn invmod_euler_998244353(n: u32) -> Option<u32> {
const MOD: u32 = 998_244_353;
let nm = n % MOD;
if nm == 0 {
// n が m の倍数なら gcd(n,m) != 1 なので、モジュラ逆数は存在しない
return None;
}
// 内部計算には64bit整数を使用(変数r,変数s)
let (mut r, mut s, mut t) = (1u64, u64::from(nm), MOD - 2);
// 繰り返し二乗法(バイナリ法)にて $n^{M-2}$ を計算
while t > 0 {
if t % 2 == 1 {
r = (r * s) % u64::from(MOD);
}
s = (s * s) % u64::from(MOD);
t /= 2;
}
// 32ビット整数に変換して出力
Some(r as u32)
}
m が累乗数の場合のモジュラ逆数 (ニュートン・ラフソン法)
自然数
と計算する事ができます。また、
となります。
証明:
また、 1回の反復で有効桁数を2倍ではなく3倍にする方法もあります。
証明:
m が 2 の累乗数 ( m=2^w ) の場合
例:
- 奇数
を 整数 n を用いて k とすると、 n=2k+1 \displaystyle n^2=8\times\left(\frac{k(k+1)}{2}\right)+1 は偶数なので、 k(k+1) は整数 \displaystyle \left(\frac{k(k+1)}{2}\right) - よって、
より n^2=8\times(\text{整数})+1 n^2\equiv 1\ (\operatorname{mod}8)
また、
は排他的論理和によるビット演算関数。 \operatorname{XOR} \displaystyle\operatorname{XOR}(3n,2)=\begin{cases}3n-2&(n\equiv 1\ (\operatorname{mod}4)\ \Leftrightarrow\ 3n\equiv 3\ (\operatorname{mod}4))\\3n+2&(n\equiv 3\ (\operatorname{mod}4)\ \Leftrightarrow\ 3n\equiv 1\ (\operatorname{mod}4))\end{cases} のとき、 n\equiv 1\ (\operatorname{mod}4) とすると n=4k+1
\begin{align*}&\quad n\times\operatorname{XOR}(3n, 2)\\&=(4k+1)(3(4k+1)-2)\\ &=32\left(k^2+\frac{k(k+1)}{2}\right)+1\\ &\equiv 1\ (\operatorname{mod}32)\end{align*}
は偶数なので、 k(k+1) は整数 \displaystyle\left(k^2+\frac{k(k+1)}{2}\right) のとき、 n\equiv 3\ (\operatorname{mod}4) とすると n=4k-1
\begin{align*}&\quad n\times\operatorname{XOR}(3n, 2)\\&=(4k-1)(3(4k-1)+2)\\ &=32\left(k^2+\frac{k(k-1)}{2}\right)+1\\ &\equiv 1\ (\operatorname{mod}32)\end{align*}
は偶数なので、 k(k-1) は整数 \displaystyle\left(k^2+\frac{k(k-1)}{2}\right) - よって、
が奇数 : n であれば、 n\equiv 1\ (\operatorname{mod}2) (n\times \operatorname{XOR}(3n,2))\equiv 1\ (\operatorname{mod}32)
x_0=n\quad\Longrightarrow\quad nx_0\equiv 1\ (\operatorname{mod}2^{\min(3,w)}) x_1=x_0(2-nx_0)\mod 2^w\quad\Longrightarrow\quad nx_1\equiv 1\ (\operatorname{mod}2^{\min(6,w)}) x_2=x_1(2-nx_1)\mod 2^w\quad\Longrightarrow\quad nx_2\equiv 1\ (\operatorname{mod}2^{\min(12,w)}) x_3=x_2(2-nx_2)\mod 2^w\quad\Longrightarrow\quad nx_3\equiv 1\ (\operatorname{mod}2^{\min(24,w)}) x_4=x_3(2-nx_3)\mod 2^w\quad\Longrightarrow\quad nx_4\equiv 1\ (\operatorname{mod}2^{\min(48,w)}) x_5=x_4(2-nx_4)\mod 2^w\quad\Longrightarrow\quad nx_5\equiv 1\ (\operatorname{mod}2^{\min(96,w)})
x_0=\operatorname{XOR}(3n,2)\mod 2^w\quad\Longrightarrow\quad nx_0\equiv 1\ (\operatorname{mod}2^{\min(5,w)}) x_1=x_0(2-nx_0)\mod 2^w\quad\Longrightarrow\quad nx_1\equiv 1\ (\operatorname{mod}2^{\min(10,w)}) x_2=x_1(2-nx_1)\mod 2^w\quad\Longrightarrow\quad nx_2\equiv 1\ (\operatorname{mod}2^{\min(20,w)}) x_3=x_2(2-nx_2)\mod 2^w\quad\Longrightarrow\quad nx_3\equiv 1\ (\operatorname{mod}2^{\min(40,w)}) x_4=x_3(2-nx_3)\mod 2^w\quad\Longrightarrow\quad nx_4\equiv 1\ (\operatorname{mod}2^{\min(80,w)})
のように
tmp = n * 3
の計算を乗算を使わず tmp = n + (n << 1)
のように変形して1サイクル以下で実行できる ) の計算に 1サイクル、乗算に 3サイクル の実行時間が掛かるとした実行時間の例 (30サイクル)。
0. `tmp0 := n * 3 mod 2^{64}` : 1cycle
1. `x0 := xor(tmp0, 2)` : 1cycle
2. `tmp1 := n * x0 mod 2^{64}` : 3cycle
3.
4.
5. `tmp2 := 2 - tmp1 mod 2^{64}` : 1cycle
6. `x1 := x0 * tmp2 mod 2^{64}` : 3cycle
7.
8.
9. `tmp3 := n * x1 mod 2^{64}` : 3cycle
10.
11.
12. `tmp4 := 2 - tmp3 mod 2^{64}` : 1cycle
13. `x2 := x1 * tmp4 mod 2^{64}` : 3cycle
14.
15.
16. `tmp5 := n * x2 mod 2^{64}` : 3cycle
17.
18.
19. `tmp6 := 2 - tmp5 mod 2^{64}` : 1cycle
20. `x3 := x1 * tmp6 mod 2^{64}` : 3cycle
21.
22.
23. `tmp7 := n * x3 mod 2^{64}` : 3cycle
24.
25.
26. `tmp8 := 2 - tmp7 mod 2^{64}` : 1cycle
27. `x4 := x1 * tmp8 mod 2^{64}` : 3cycle
28.
29.
30. `return x4`
Jeffrey Hurchalla’s method
さらにこの漸化式
x_0=\operatorname{XOR}(3n,2)\mod 2^w,\quad y_0=1-nx_0\mod 2^w\quad\Longrightarrow\quad nx_0\equiv 1\ (\operatorname{mod}2^{\min(5,w)}) x_1=x_0(1+y_0)\mod 2^w,\quad y_1=(y_0)^2\mod 2^w\quad\Longrightarrow\quad nx_1\equiv 1\ (\operatorname{mod}2^{\min(10,w)}) x_2=x_1(1+y_1)\mod 2^w,\quad y_2=(y_1)^2\mod 2^w\quad\Longrightarrow\quad nx_2\equiv 1\ (\operatorname{mod}2^{\min(20,w)}) x_3=x_2(1+y_2)\mod 2^w,\quad y_3=(y_2)^2\mod 2^w\quad\Longrightarrow\quad nx_3\equiv 1\ (\operatorname{mod}2^{\min(40,w)}) x_4=x_3(1+y_3)\mod 2^w,\quad y_4=(y_3)^2\mod 2^w\quad\Longrightarrow\quad nx_4\equiv 1\ (\operatorname{mod}2^{\min(80,w)})
このように ニュートン・ラフソン法の漸化式
nx_0\equiv(1-y_0)\quad(\operatorname{mod}2^w) x_1\equiv x_0(2-nx_0)\equiv x_0(1+y_0)\quad(\operatorname{mod}2^w) nx_1\equiv nx_0(1+y_0)\equiv(1-y_0)(1+y_0)\equiv(1-(y_0)^2)\quad(\operatorname{mod}2^w) x_2\equiv x_1(2-nx_1)\equiv x_1(1+(y_0)^2)\quad(\operatorname{mod}2^w) nx_2\equiv nx_1(1+(y_0)^2)\equiv(1-(y_0)^2)(1+(y_0)^2)\equiv(1-(y_0)^4)\quad(\operatorname{mod}2^w) x_3\equiv x_2(2-nx_2)\equiv x_2(1+(y_0)^4)\quad(\operatorname{mod}2^w) nx_3\equiv nx_2(1+(y_0)^4)\equiv(1-(y_0)^4)(1+(y_0)^4)\equiv(1-(y_0)^8)\quad(\operatorname{mod}2^w) x_4\equiv x_3(2-nx_3)\equiv x_3(1+(y_0)^8)\quad(\operatorname{mod}2^w) nx_4\equiv nx_3(1+(y_0)^8)\equiv(1-(y_0)^8)(1+(y_0)^8)\equiv(1-(y_0)^{16})\quad(\operatorname{mod}2^w)
tmp = n * 3
の計算を乗算を使わず tmp = n + (n << 1)
のように変形して1サイクル以下で実行できる ) の計算に 1サイクル、乗算に 3サイクル の実行時間が掛かるとして、以下のように 19サイクル で
0. `tmp0 := n * 3 mod 2^{64}` : 1cycle
1. `x0 := xor(tmp0, 2)` : 1cycle
2. `tmp1 := n * x0 mod 2^{64}` : 3cycle
3.
4.
5. `y0 := 1 - tmp1 mod 2^{64}` : 1cycle
6. `tmp2 := y0 + 1 mod 2^{64}` : 1cycle, `y1 := y0 * y0 mod 2^{64}`: 3cycle
7. `x1 := x0 * tmp2 mod 2^{64}` : 3cycle,
8.
9. `tmp3 := y1 + 1 mod 2^{64}` : 1cycle, `y2 := y1 * y1 mod 2^{64}` : 3cycle
10. `x2 := x1 * tmp3 mod 2^{64}` : 3cycle
11.
12. `tmp4 := y2 + 1 mod 2^{64}` : 1cycle, `y3 := y2 * y2 mod 2^{64}` : 3cycle
13. `x3 := x2 * tmp4 mod 2^{64}` : 3cycle
14.
15. `tmp5 := y3 + 1 mod 2^{64}` : 1cycle
16. `x4 := x3 * tmp5 mod 2^{64}` : 3cycle
17.
18.
19. `return x4`
サンプルコード (Rust)
// モジュラ逆数 $n^{-1} \mod 2^{32}$ の計算
pub fn invmod_u32(n: u32) -> Option<u32> {
// n が偶数だとモジュラ逆数は存在しない
if n % 2 == 0 {
return None;
}
// Jeffrey Hurchalla’s method: https://arxiv.org/abs/2204.04342
let x0 = n.wrapping_mul(3) ^ 2; // x0 := BITWISEXOR(n * 3, 2) mod 2^{32}
let y0 = 1u32.wrapping_sub(n.wrapping_mul(x0)); // y0 := 1 - n * x0 mod 2^{32}
let x1 = x0.wrapping_mul(y0.wrapping_add(1)); // x1 := x0 * (y0 + 1) mod 2^{32}
let y1 = y0.wrapping_mul(y0); // y1 := y0 * y0 mod 2^{32}
let x2 = x1.wrapping_mul(y1.wrapping_add(1)); // x2 := x1 * (y1 + 1) mod 2^{32}
let y2 = y1.wrapping_mul(y1); // y2 := y1 * y1 mod 2^{32}
let x3 = x2.wrapping_mul(y2.wrapping_add(1)); // x3 := x2 * (y2 + 1) mod 2^{32}
Some(x3) // return x3
}
// モジュラ逆数 $n^{-1} \mod 2^{64}$ の計算
pub fn invmod_u64(n: u64) -> Option<u64> {
// n が偶数だとモジュラ逆数は存在しない
if n % 2 == 0 {
return None;
}
// Jeffrey Hurchalla’s method: https://arxiv.org/abs/2204.04342
let x0 = n.wrapping_mul(3) ^ 2; // x0 := BITWISEXOR(n * 3, 2) mod 2^{64}
let y0 = 1u64.wrapping_sub(n.wrapping_mul(x0)); // y0 := 1 - n * x0 mod 2^{64}
let x1 = x0.wrapping_mul(y0.wrapping_add(1)); // x1 := x0 * (y0 + 1) mod 2^{64}
let y1 = y0.wrapping_mul(y0); // y1 := y0 * y0 mod 2^{64}
let x2 = x1.wrapping_mul(y1.wrapping_add(1)); // x2 := x1 * (y1 + 1) mod 2^{64}
let y2 = y1.wrapping_mul(y1); // y2 := y1 * y1 mod 2^{64}
let x3 = x2.wrapping_mul(y2.wrapping_add(1)); // x3 := x2 * (y2 + 1) mod 2^{64}
let y3 = y2.wrapping_mul(y2); // y3 := y2 * y2 mod 2^{64}
let x4 = x3.wrapping_mul(y3.wrapping_add(1)); // x4 := x3 * (y3 + 1) mod 2^{64}
Some(x4) // return x4
}
m に関するモジュラ逆数 (変則・拡張ユークリッド互除法)
一般の 拡張ユークリッドの互除法によるモジュラ逆数の計算は、多くの資料では符号付き整数型を用いた実装が多いですが、ここでは少し式を工夫して符号無し整数型を用いた実装を考えてみます。
符号付き整数型を用いた拡張ユークリッド互除法のコード例 (Rust):
// extgcd(a, b) -> (g, x, y)
// ax + by = gcd(a,b) = g を満たす整数の組 (g, x, y) をひとつ返す
pub fn extgcd(a: i64, b: i64) -> (i64, i64, i64) {
match b {
0 => (a, 1, 0),
_ => {
let (d, y, x) = extgcd(b, a % b);
(d, x, y - (a / b) * x)
},
}
}
対象
このような形の連立方程式を変形する操作を考えてみます。
(a_in-c_im)/2^{s_i}=x_i (b_in-d_im)/2^{s_i}=-y_i
目標
\displaystyle x_i=1\quad\Longrightarrow\quad a_in\equiv 2^{s_i}\quad(\operatorname{mod}m) \quad\Longrightarrow\quad n^{-1}\equiv a_i/2^{s_i}\quad(\operatorname{mod}m) \displaystyle y_i=1\quad\Longrightarrow\quad b_in\equiv -2^{s_i}\quad(\operatorname{mod}m) \quad\Longrightarrow\quad n^{-1}\equiv -b_i/2^{s_i}\quad(\operatorname{mod}m)
\displaystyle x_i=1\quad\Longrightarrow\quad a_in\equiv 1\quad(\operatorname{mod}m) \quad\Longrightarrow\quad n^{-1}\equiv a_i\quad(\operatorname{mod}m) \displaystyle y_i=1\quad\Longrightarrow\quad b_in\equiv -1\quad(\operatorname{mod}m) \quad\Longrightarrow\quad n^{-1}\equiv -b_i\quad(\operatorname{mod}m)
制約
\{n,m,a_i,b_i,c_i,d_i,s_i,t_i,x_i,y_i\} は非負整数 \gcd(n,m)=1
初期値
(a_0n-c_0m)/2^{s_0}=x_0=(1\cdot n-0\cdot m)/1=n (b_0n-d_0m)/2^{s_0}=-y_0=(0\cdot n-1\cdot m)/1=-m
式変形パターン
A.
(a_{i+1}n-c_{i+1}m)/2^{s_{i+1}}=x_{i+1}=((a_i+b_it_i)n-(c_i+d_it_i)m)/2^{s_i}=x_i-t_iy_i (b_{i+1}n-d_{i+1}m)/2^{s_{i+1}}=-y_{i+1}=(b_in-d_im)/2^{s_i}=-y_i
B.
(a_{i+1}n-c_{i+1}m)/2^{s_{i+1}}=x_{i+1}=(a_in-c_im)/2^{s_i}=x_i (b_{i+1}n-d_{i+1}m)/2^{s_{i+1}}=-y_{i+1}=((a_it_i+b_i)n-(c_it_i+d_i)m)/2^{s_i}=t_ix_i-y_i
C.
(a_{i+1}n-c_{i+1}m)/2^{s_{i+1}}=x_{i+1}=(a_in-c_im)/2^{s_i+t_i}=x_i/2^{t_i} (b_{i+1}n-d_{i+1}m)/2^{s_{i+1}}=-y_{i+1}=(2^{t_i}b_in-2^{t_i}d_im)/2^{s_i+t_i}=-y_i
D.
(a_{i+1}n-c_{i+1}m)/2^{s_{i+1}}=x_{i+1}=(2^{t_i}a_in-2^{t_i}c_im)/2^{s_i+t_i}=x_i (b_{i+1}n-d_{i+1}m)/2^{s_{i+1}}=-y_{i+1}=(b_in-d_im)/2^{s_i+t_i}=-y_i/2^{t_i}
m について、 A.B. のみの式変形を使う場合 )
式変形手順の例 ( 一般の A.B. のみの式変形を使う事で、
-
と初期化 ((a,b,x,y)=(1,0,n,m) は使わない )c,d,s -
であれば異常終了 :x=0 のためモジュラ逆数が存在しない\gcd(n,m)=\gcd(x,y)=y\ge 2 -
であれば正常終了 :x=1 を計算して出力\displaystyle n^{-1}\equiv a\quad(\operatorname{mod}m) -
であればx\le y として 式変形B. を実行t=\lfloor y/x\rfloor -
であれば異常終了 :y=0 のためモジュラ逆数が存在しない\gcd(n,m)=\gcd(x,y)=x\ge 2 -
であれば正常終了 :y=1 を計算して出力\displaystyle n^{-1}\equiv -b\quad(\operatorname{mod}m) -
であればx\ge y として 式変形A. を実行t=\lfloor x/y\rfloor - 手順 2. に戻る
サンプルコード (Rust):
// 拡張ユークリッド互除法によるモジュラ逆数 $n^{-1} \mod m$
// 符号なし整数型を用いたバリエーション
pub fn modinv_plain(n: u32, m: u32) -> Option<u32> {
let (mut a, mut b, mut x, mut y) = (1, 0, n, m);
if m == 1 {
return Some(0);
}
loop {
match x {
// $\gcd(n,m) = \gcd(x,y) = y \ge 2$
0 => return None,
// $n^{-1} \equiv a$
1 => return Some(a),
_ => {}
}
b += a * (y / x);
y %= x;
match y {
// $\gcd(n,m) = \gcd(x,y) = x \ge 2$
0 => return None,
// $n^{-1} \equiv -b$
1 => return Some(m - b),
_ => {}
}
a += b * (x / y);
x %= y;
}
}
m が 3 以上の奇数 の場合、バイナリGCD法の複合 )
式変形手順の具体例 (
- 例えば
の場合、0\le n\lt m\lt 2^{32} の\lbrace 2^0, 2^{-1}, 2^{-2}, ..., 2^{-63}\rbrace をそれぞれ事前計算しておく ( 計算方法は後述の補題参照、\operatorname{mod}m が取りうる値くらいの数だけ用意 )\lfloor\log_2 mn\rfloor -
と初期化 ((a,b,s,x,y)=(1,0,0,n,m) は使わない )c,d -
であれば異常終了 :x=0 のためモジュラ逆数が存在しない\gcd(n,m)=\gcd(x,y)=y\ge 3 -
の偶数であればx\ge 2 を奇数にするような 自然数x にて 式変形C. を実行t -
であれば正常終了 :x=1 を計算して出力\displaystyle n^{-1}\equiv 2^{-s}a\quad(\operatorname{mod}m) -
であれば異常終了 :x=y のためモジュラ逆数が存在しない\gcd(n,m)=\gcd(x,y)=x=y\ge 3 -
であればx\lt 8y として 式変形B. を実行t=\lfloor y/x\rfloor -
であればx\lt y として式変形B. を実行t=1 -
であれば異常終了 :y=0 のためモジュラ逆数が存在しない\gcd(n,m)=\gcd(x,y)=x\ge 3 -
の偶数であればy\ge 2 を奇数にするような 自然数y にて 式変形D. を実行t -
であれば正常終了 :y=1 を計算して出力\displaystyle n^{-1}\equiv -2^{-t}b\quad(\operatorname{mod}m) -
であれば8x\gt y として 式変形A. を実行t=\lfloor x/y\rfloor -
であればx\gt y として式変形A. を実行t=1 - 手順 3. に戻る
サンプルコード (Rust): 記事末尾のサンプルコード fn inv_mx
関数の部分を参照。
m が 3 以上の奇数の場合の 2 の累乗数に対するモジュラ逆数の事前計算
補題: 1/2^0\ \operatorname{mod}m=1 \displaystyle 1/2^{q+1}\ \operatorname{mod}m=\begin{cases}(1/2^{q}\ \operatorname{mod}m)/2&((1/2^{q}\ \operatorname{mod}m)\equiv 0\ \operatorname{mod}2)\\ ((1/2^{q}\ \operatorname{mod}m)+m)/2&((1/2^{q}\ \operatorname{mod}m)\equiv 1\ \operatorname{mod}2)\end{cases}
例:
m=1 の場合
-
より、nn^{-1} \equiv 0 \equiv 1\ (\operatorname{mod}1) を出力n^{-1}\equiv 0\ (\operatorname{mod}1)
m=0 の場合
- 未定義動作とする。異常終了しても良いし、
やm=2^{32} などとみなして処理しても良いし、鼻から悪魔が出ても良い。m=2^{64}
m の場合、主に m が 2 以上の偶数 の場合)
ニュートン・ラフソン法、中国剰余定理 を組み合わせる (任意の -
とする。 ここで、m=m_xm_y,\ m_x=d,\ m_y=2^r はr,d の形に分解した時の値 (m=2^r\times d は非負整数、r は正の奇数 ) 。d は正の奇数、m_x はm_y の累乗数 ( ただし2 が奇数であれば、m ) であり、m_y=1 を満たす。また、\gcd(m_x, m_y)=1 であれば、\gcd(n,m)=1 を満たす。\gcd(n,m_x)=1,\ \gcd(n,m_y)=1 -
を満たす、m_x\overline{m_x}\equiv 1\quad(\operatorname{mod}m_y) における\operatorname{mod}m_y のモジュラ逆数m_x の値を、前述の「ニュートン・ラフソン法」で求める。\overline{m_x} の時は、m_y=1 とする。\overline{m_x}=0 -
における\operatorname{mod}m_x のモジュラ逆数n の値を、前述の「変則的拡張ユークリッドの互除法・式変形手順の具体例 (x=n^{-1}\ \operatorname{mod}m_x がm 以上の奇数 の場合 )」で求める。3 の時は、m_x=1 とする。x=0 であれば、\gcd(n,m_x)\ne 1 であるので、モジュラ逆数\gcd(n,m)\ne 1 の値が存在しないとして異常終了する。n^{-1}\ \operatorname{mod}m -
における\operatorname{mod}m_y のモジュラ逆数n の値を、前述の「ニュートン・ラフソン法」で求める。y=n^{-1}\ \operatorname{mod}m_y の時は、m_y=1 とする。y=0 であれば、\gcd(n,m_y)\ne 1 であるので、モジュラ逆数\gcd(n,m)\ne 1 の値が存在しないとして異常終了する。n^{-1}\ \operatorname{mod}m -
であるので、m=m_xm_y,\ \gcd(m_x,m_y)=1 をそれぞれ満たすような 整数z\equiv x\ (\operatorname{mod}m_x),\ z\equiv y\ (\operatorname{mod}m_y) の値を、後述の「中国剰余定理 (Garner のアルゴリズム)」を用いて、0\le z\lt m もしくはz=(x-m_x((\overline{m_x}(x-y))\operatorname{mod}m_y))\operatorname{mod}m のように求める。z=(x+m_x((\overline{m_x}(y-x))\operatorname{mod}m_y))\operatorname{mod}m -
よりm=m_xm_y,\ \gcd(m_x,m_y)=1, nz\equiv 1\ (\operatorname{mod}m_x),\ nz\equiv 1\ (\operatorname{mod}m_y) であるので、nz\equiv 1\ (\operatorname{mod}m) を モジュラ逆数z の値として正常終了する。n^{-1}\ \operatorname{mod}m
補題: 中国剰余定理 (Garner のアルゴリズム)
x,y\in\text{整数}\mathbb{Z} M_x,M_y,\overline{M_x}\in\text{自然数}\mathbb{N} \gcd(M_x,M_y)=1 \overline{M_x}M_x\equiv 1\quad(\operatorname{mod}M_y) t=M_x((\overline{M_x}(x-y))\operatorname{mod}M_y) z=\begin{cases}x-t&(x \ge t)\\x-t+M_xM_y&(x \lt t)\\\end{cases}
とすると
0\le t\le M_x(M_y-1) x-M_x(M_y-1)\le x-t\le x+M_xM_y z\equiv x\quad(\operatorname{mod}M_x)\quad\because\quad t\equiv 0\quad(\operatorname{mod}M_x) z\equiv y\quad(\operatorname{mod}M_y)\quad\because\quad t\equiv x-y\quad(\operatorname{mod}M_y) - 特に、
であれば0\le x\lt M_xM_y 0\le z\lt M_xM_y
サンプルコード (Rust)
拡張ユークリッド互除法 + バイナリGCD法 + ニュートン・ラフソン法 + Garnerのアルゴリズム を組み合わせてモジュラ逆元を計算するコードの例です。
// Modulus値分類
pub enum ModIntDynType {
Zero,
Odd,
P0998244353,
P1000000007,
Power2,
Even,
}
// モジュラ演算用定数構造体
pub struct ModIntDyn {
ty: ModIntDynType,
m: u32,
tz: u32,
mx: u32,
mimx: u32,
mym: u32,
dim: u64,
dimx: u64,
ib: [u32; 64],
}
impl ModIntDyn {
pub fn new_modulus(m: u32) -> Self {
if m == 0 {
// $m = 0$ であれば $m = 2^{32}$ とみなして処理
return Self {
m,
ty: ModIntDynType::Zero,
tz: 32,
mx: 1,
mimx: 0,
mym: !0u32,
dim: 1,
dimx: 0,
ib: [1u32; 64],
};
}
// tz = tzcnt(m), my = 2^tz
let tz = m.trailing_zeros();
// mx = modulus / my
let mx = m >> tz;
// mx * imx % 2^{32} = 1 => mx * imx % my = 1 % my
let mimx = {
let mut t = mx.wrapping_mul(3) ^ 2;
let mut u = 1u32.wrapping_sub(t.wrapping_mul(mx));
t = t.wrapping_mul(u.wrapping_add(1));
u = u.wrapping_mul(u);
t = t.wrapping_mul(u.wrapping_add(1));
u = u.wrapping_mul(u);
t = t.wrapping_mul(u.wrapping_add(1));
t
};
// mym = my - 1, x % my == x & mym
let mym = u32::from(tz != 0).wrapping_neg() >> ((32 - tz) & 31);
// `dim` $\ = \lceil (2^{64} / m) \rceil = \lfloor (2^{64}-1) / m \rfloor + 1$
let dim = ((!0u64) / u64::from(m)).wrapping_add(1);
// `dimx` $\ = \lceil (2^{64} / m_x) \rceil = \lfloor (2^{64}-1) / m_x \rfloor + 1$
let dimx = ((!0u64) / u64::from(mx)).wrapping_add(1);
// ib: $\{2^0,2^{-1},\dots,2^{-63}\} \mod m_x$ を事前計算
let mut ib = [0u32; 64];
let mut t = 1;
for e in ib.iter_mut() {
*e = t;
t = t / 2 + ((t % 2).wrapping_neg() & (mx / 2 + 1));
}
// オブジェクト生成
Self {
m,
ty: {
if m == 998_244_353 {
ModIntDynType::P0998244353
} else if m == 1_000_000_007 {
ModIntDynType::P1000000007
} else if tz == 0 {
ModIntDynType::Odd
} else if mx == 1 {
ModIntDynType::Power2
} else {
ModIntDynType::Even
}
},
tz,
mx,
mimx,
mym,
dim,
dimx,
ib,
}
}
// $a \mod m$ (normalize)
#[inline]
pub fn norm(&self, a: u32) -> u32 {
match self.ty {
// $m = 2^{32}$
ModIntDynType::Zero => a,
// $m = m_y$
ModIntDynType::Power2 => a & self.mym,
// others
_ => a % self.m,
}
}
// $a + b \mod m$
#[inline]
pub fn add(&self, a: u32, b: u32) -> u32 {
debug_assert!(a < self.m);
debug_assert!(b < self.m);
let v = u64::from(a) + u64::from(b);
match self.ty {
// $m = 2^{32}
ModIntDynType::Zero => v as u32,
// others
_ => match v.overflowing_sub(u64::from(self.m)) {
(_, true) => v as u32,
(w, false) => w as u32,
},
}
}
// $a - b \mod m$
#[inline]
pub fn sub(&self, a: u32, b: u32) -> u32 {
debug_assert!(a < self.m);
debug_assert!(b < self.m);
let (v, f) = a.overflowing_sub(b);
match self.ty {
// $m = 2^{32}
ModIntDynType::Zero => v,
// others
_ => v.wrapping_add(u32::from(f).wrapping_neg() & self.m),
}
}
// $-a \mod m$
#[inline]
pub fn neg(&self, a: u32) -> u32 {
debug_assert!(a < self.m);
self.sub(0, a)
}
// $a * b \mod m$
#[inline]
pub fn mul(&self, a: u32, b: u32) -> u32 {
debug_assert!(a < self.m);
debug_assert!(b < self.m);
let z = u64::from(a) * u64::from(b);
match self.ty {
// $m = m_y$
ModIntDynType::Power2 => (z as u32) & self.mym,
// $m = 998244353$
ModIntDynType::P0998244353 => {
(z - (((u128::from(z) * 0x89AE_4087_5DE0_CC3F) >> 93) as u64) * 998_244_353) as u32
}
// $m = 1000000007$
ModIntDynType::P1000000007 => {
(z - (((u128::from(z) * 0x8970_5F31_12A2_8FE5) >> 93) as u64) * 1_000_000_007)
as u32
}
// $m = 2^{32}$
ModIntDynType::Zero => z as u32,
// others
ModIntDynType::Odd | ModIntDynType::Even => {
// Barrett Reduction
let x = ((u128::from(z) * u128::from(self.dim)) >> 64) as u64;
let (v, f) = z.overflowing_sub(x.wrapping_mul(u64::from(self.m)));
(v as u32).wrapping_add(u32::from(f).wrapping_neg() & self.m)
}
}
}
// $a * b \mod m_x$
#[inline]
fn mulmx(&self, a: u32, b: u32) -> u32 {
debug_assert!(a < self.mx);
debug_assert!(b < self.mx);
let z = u64::from(a) * u64::from(b);
match self.ty {
// $m = 998244353$
ModIntDynType::P0998244353 => {
(z - (((u128::from(z) * 0x89AE_4087_5DE0_CC3F) >> 93) as u64) * 998_244_353) as u32
}
// $m = 1000000007$
ModIntDynType::P1000000007 => {
(z - (((u128::from(z) * 0x8970_5F31_12A2_8FE5) >> 93) as u64) * 1_000_000_007)
as u32
}
// $m_x = 0$
ModIntDynType::Power2 | ModIntDynType::Zero => 0,
// others
ModIntDynType::Odd | ModIntDynType::Even => {
// Barrett Reduction
let x = ((u128::from(z) * u128::from(self.dimx)) >> 64) as u64;
let (v, f) = z.overflowing_sub(x.wrapping_mul(u64::from(self.mx)));
(v as u32).wrapping_add(u32::from(f).wrapping_neg() & self.mx)
}
}
}
// $a^b \mod m$
#[inline]
pub fn pow(&self, a: u32, mut b: u32) -> u32 {
debug_assert!(a < self.m);
if b == 0 {
return (self.m > 1) as u32;
}
let mut t = a;
while b % 2 == 0 {
t = self.mul(t, t);
b /= 2;
}
let mut r = t;
b /= 2;
while b > 0 {
t = self.mul(t, t);
if b % 2 != 0 {
r = self.mul(r, t);
}
b /= 2;
}
r
}
// $a / b \mod m$
#[inline]
pub fn div(&self, a: u32, b: u32) -> Option<u32> {
debug_assert!(a < self.m);
debug_assert!(b < self.m);
self.inv(b).map(|ib| self.mul(a, ib))
}
// モジュラ逆数の計算(非負整数による拡張ユークリッド互除法の実装例)
#[inline]
pub fn inv_plain(&self, n: u32) -> Option<u32> {
let m = self.m;
let (mut a, mut b, mut x, mut y) = (1, 0, n, m);
if m == 1 {
return Some(0);
}
loop {
match x {
// $\gcd(n,m) = \gcd(x,y) = y \ge 2$
0 => return None,
// $n^{-1} \equiv a$
1 => return Some(a),
_ => {}
}
b += a * (y / x);
y %= x;
match y {
// $\gcd(n,m) = \gcd(x,y) = x \ge 2$
0 => return None,
// $n^{-1} \equiv -b$
1 => return Some(m - b),
_ => {}
}
a += b * (x / y);
x %= y;
}
}
// inv の子関数: 拡張ユークリッド互除法 + バイナリGCD法 による モジュラ逆数の計算 ($\mod m_x$) ($m_x$ は 正の奇数)
#[inline]
fn inv_mx(&self, n: u32) -> Option<u32> {
debug_assert_ne!(self.mx % 2, 0);
if self.mx % 2 == 0 {
unsafe { core::hint::unreachable_unchecked() }
}
if self.mx == 1 {
return Some(0);
}
let (mut a, mut b, mut x, mut y) = (1, 0, n, self.mx);
if x == 0 {
// $\gcd(n,m_x) = \gcd(x,y) = m_x \ge 3$
return None;
}
// 式変形C: $t = \operatorname{tzcnt}(x)
let mut s = x.trailing_zeros();
x >>= s;
loop {
if x == 1 {
// $n^{-1} \equiv a / 2^s \mod m_x$
return Some(self.mulmx(a, self.ib[s as usize]));
}
if x == y {
// $\gcd(n,m_x) = \gcd(x,y) = x = y \ge 3$
return None;
}
if x <= (y >> 3) {
// 式変形B: $t = \lfloor y / x \rfloor$
if x == 0 {
unsafe { core::hint::unreachable_unchecked() }
}
b += a * (y / x);
y %= x;
if y == 0 {
// $\gcd(n,m_x) = \gcd(x,y) = x \ge 3$
return None;
}
// 式変形D: $t = \operatorname{tzcnt}(y)$
let t = y.trailing_zeros();
y >>= t;
a <<= t;
s += t;
} else {
while x < y {
// 式変形B: $t = 1$
b += a;
y -= x;
// 式変形D: $t = \operatorname{tzcnt}(y)$
let t = y.trailing_zeros();
y >>= t;
a <<= t;
s += t;
}
}
if y == 1 {
// $n^{-1} \equiv -b / 2^s \mod m_x$
return Some(self.mulmx(self.mx - b, self.ib[s as usize]));
}
if (x >> 3) >= y {
// 式変形A: $t = \lfloor x / y \rfloor$
if y == 0 {
unsafe { core::hint::unreachable_unchecked() }
}
a += b * (x / y);
x %= y;
if x == 0 {
// $\gcd(n,m_x) = \gcd(x,y) = y \ge 3$
return None;
}
// 式変形C: $t = \operatorname{tzcnt}(x)$
let t = x.trailing_zeros();
x >>= t;
b <<= t;
s += t;
} else {
while x > y {
// 式変形A: $t = 1$
a += b;
x -= y;
// 式変形C: $t = \operatorname{tzcnt}(x)$
let t = x.trailing_zeros();
x >>= t;
b <<= t;
s += t;
}
}
}
}
// inv の子関数: ニュートン・ラフソン法 による モジュラ逆数 の計算 ($\mod m_y$) ($m_y$ は 1 もしくは 2 の累乗数)
#[inline]
fn inv_my(&self, n: u32) -> Option<u32> {
if self.mym == 0 {
// 0bit, $m_y == 1$
return Some(0);
}
if n % 2 == 0 {
// $\gcd(n,m_y) \ge 2$
return None;
}
Some(
match self.tz {
// 21bit ~ 32bit
21..=0xffff_ffff => {
let mut t = n.wrapping_mul(3) ^ 2;
let mut u = 1u32.wrapping_sub(t.wrapping_mul(n));
t = t.wrapping_mul(u.wrapping_add(1));
u = u.wrapping_mul(u);
t = t.wrapping_mul(u.wrapping_add(1));
u = u.wrapping_mul(u);
t = t.wrapping_mul(u.wrapping_add(1));
t
}
// 11bit ~ 20bit
11..=20 => {
let mut t = n.wrapping_mul(3) ^ 2;
let mut u = 1u32.wrapping_sub(t.wrapping_mul(n));
t = t.wrapping_mul(u.wrapping_add(1));
u = u.wrapping_mul(u);
t = t.wrapping_mul(u.wrapping_add(1));
t
}
// 6bit ~ 10bit
6..=10 => {
let t = n.wrapping_mul(3) ^ 2;
t.wrapping_mul(2u32.wrapping_sub(n.wrapping_mul(t)))
}
// 4bit ~ 5bit
4..=5 => n.wrapping_mul(3) ^ 2,
// 1bit ~ 3bit
0..=3 => n,
} & self.mym,
)
}
// モジュラ逆数の計算(拡張ユークリッドの互除法、バイナリGCD法、ニュートン・ラフソン法、Garnerのアルゴリズムのハイブリッド)
#[inline]
pub fn inv(&self, n: u32) -> Option<u32> {
match self.ty {
ModIntDynType::Power2 => self.inv_my(n),
ModIntDynType::P0998244353 | ModIntDynType::P1000000007 | ModIntDynType::Odd => {
self.inv_mx(n)
}
ModIntDynType::Even => {
if let Some(y) = self.inv_my(n) {
// 拡張ユークリッドの互除法+バイナリGCD法 にて $n^{-1} \mod m_x$ を計算
if let Some(x) = self.inv_mx(n) {
// 中国剰余定理・Garner のアルゴリズム による結果の合成
let (t, f) = x.overflowing_sub(
self.mx * (self.mimx.wrapping_mul(x.wrapping_sub(y)) & self.mym),
);
return Some(t.wrapping_add(u32::from(f).wrapping_neg() & self.m));
}
}
None
}
ModIntDynType::Zero => {
if n % 2 == 0 {
return None;
}
let mut t = n.wrapping_mul(3) ^ 2;
let mut u = 1u32.wrapping_sub(t.wrapping_mul(n));
t = t.wrapping_mul(u.wrapping_add(1));
u = u.wrapping_mul(u);
t = t.wrapping_mul(u.wrapping_add(1));
u = u.wrapping_mul(u);
t = t.wrapping_mul(u.wrapping_add(1));
Some(t)
}
}
}
}
// 1 / values[0] / values[1] / ... (mod m)
pub fn moddyn_div_test(modulus: u32, values: &[u32]) -> u32 {
let moddyn = ModIntDyn::new_modulus(modulus);
values.iter().fold(1u32, |p, &c| moddyn.div(p, c).unwrap())
}
Discussion