ピン留めされたアイテム

モジュラ逆数の計算

$0\lt n \lt m\le 2^{64}$ 辺りの範囲で $nn^{-1}\equiv 1\ (\operatorname{mod}m)$ が成り立つモジュラ逆数 $n^{-1}\mod m$ を求める方法を考えてみます。

(サンプルコードでは主に $0\lt n\lt m\lt 2^{32}$ を前提として書きます)

モジュラ逆数 $n^{-1}\mod m$ が存在するための必要十分条件は、 $n$ と $m$ の最大公約数 $\gcd(n,m)$ が $\gcd(n,m)=1$ を満たすことです。

ここでは、3つの方法でのモジュラ逆数の求め方を紹介します。

オイラーの定理を用いる方法 ( $m$ が素数である場合、 $m$ の素因数分解が既知である場合 )
変則的な拡張ユークリッド互除法を用いる方法 ( 一般の自然数 $m$ に対する場合 )
ニュートン・ラフソン法を用いる方法 ( $m$ が累乗数である場合 )

$m$ が素数、もしくは $m$ の素因数分解が既知である場合のモジュラ逆数 (オイラーの定理)

オイラーのトーシェント関数 $\varphi(m)$ を用いると、オイラーの定理より

$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)$

$n^{\varphi(m)-1}\ (\operatorname{mod}m)$ は、繰り返し二乗法（バイナリ法）などを使って計算できます。

$m$ が素数であれば、 $\varphi(m)=m-1$ であるので、 $\gcd(n,m)=1$ ならば $n^{-1}\equiv n^{m-2}\ (\operatorname{mod}m)$ と計算できます。

$m$ が素数でなければ、 $\varphi(m)$ の値は以下のように計算する必要があります。

オイラーのトーシェント関数

オイラーのトーシェント関数 $\varphi(m)$ は、正の整数 $m$ に対して $m$ と互いに素である $1$ 以上 $m$ 以下の自然数の個数です。

\displaystyle\varphi(m)=\sum_{1\le n\le m\atop \gcd(n,m)}1

$m$ が素数であれば、 $\varphi(m)=m-1$
$m$ が素数ではなく、 $m$ の素因数分解が以下のように ( $p_1,p_2,p_3,\dots,p_d$ は互いに異なる素数 ) 既知であれば、

\displaystyle m=\prod_{k=1}^d p_k^{\,e_k}=p_1^{\,e_1}\times p_2^{\,e_2}\times p_3^{\,e_3}\times\cdots\times p_d^{\,e_d}

\varphi(m)=m\prod_{k=1}^d\left(1-\frac{1}{p_k}\right)=m\left(1-\frac{1}{p_1}\right)\left(1-\frac{1}{p_2}\right)\left(1-\frac{1}{p_3}\right)\cdots\left(1-\frac{1}{p_d}\right)

サンプルコード (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$ に関するモジュラ逆数 (変則・拡張ユークリッド互除法)

拡張ユークリッドの互除法によるモジュラ逆数の計算は、多くの資料では符号付き整数型を用いた実装が多いですが、ここでは少し式を工夫して符号無し整数型を用いた実装を考えてみます。

符号付き整数型を用いた拡張ユークリッド互除法のコード例 (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$

目標

$\gcd(n,m)=1$ である事を前提に、 $i$ 回の操作後に $x_i=1$ もしくは $y_i=1$ を満たすように式変形が出来たとすると、

$m$ が奇数の場合、モジュラ逆数 $n^{-1}\mod m$ を以下のように導出します。

$\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)$

$m$ が偶数の場合、 $2^{s_i}$ の値が常に $1$ ( つまり、 $s_i=0$ ) となるように変形を行い、モジュラ逆数 $n^{-1}\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$

初期値

$\displaystyle\begin{pmatrix}s_0&a_0&c_0&x_0\\&b_0&d_0&y_0\end{pmatrix}=\begin{pmatrix}0&1&0&n\\&0&1&m\end{pmatrix}$

$(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. $x_i \ge t_iy_i$ : 1行目の式に2行目の式の $t_i$ 倍を加える操作

$\gcd(n,m)=\gcd(x_i,y_i)=\gcd(x_i-t_iy_i,y_i)=\gcd(x_{i+1},y_{i+1})$

$(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$

\begin{pmatrix}s_{i+1}&a_{i+1}&c_{i+1}&x_{i+1}\\&b_{i+1}&d_{i+1}&y_{i+1}\end{pmatrix}= \begin{pmatrix}s_i&a_i+b_it_i&c_i+d_it_i&x_i-t_iy_i\\&b_i&d_i&y_i\end{pmatrix}

B. $t_ix_i \le y_i$ : 2行目の式に1行目の式の $t_i$ 倍を加える操作

$\gcd(n,m)=\gcd(x_i,y_i)=\gcd(x_i,y_i-t_ix_i)=\gcd(x_{i+1},y_{i+1})$

$(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$

\begin{pmatrix}s_{i+1}&a_{i+1}&c_{i+1}&x_{i+1}\\&b_{i+1}&d_{i+1}&y_{i+1}\end{pmatrix}= \begin{pmatrix}s_i&a_i&c_i&x_i\\&a_it_i+b_i&c_it_i+d_i&y_i-t_ix_i\end{pmatrix}

C. $m\equiv 1\ (\operatorname{mod}2)$ かつ $x_i$ が $2^{t_i}$ の倍数 : $x_i$ の値を $1/2^{t_i}$ する操作

$x_i$ が $2^{t_i}$ の倍数、 $y_i$ が奇数であれば $\gcd(n,m)=\gcd(x_i,y_i)=\gcd(x_i/2^{t_i},y_i)=\gcd(x_{i+1},y_{i+1})$

$(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$

\begin{pmatrix}s_{i+1}&a_{i+1}&c_{i+1}&x_{i+1}\\&b_{i+1}&d_{i+1}&y_{i+1}\end{pmatrix}= \begin{pmatrix}s_i+t_i&a_i&c_i&x_i/2^{t_i}\\&2^{t_i}b_i&2^{t_i}d_i&y_i\end{pmatrix}

D. $m\equiv 1\ (\operatorname{mod}2)$ かつ $y_i$ が $2^{t_i}$ の倍数 : $y_i$ の値を $1/2^{t_i}$ する操作

$x_i$ が奇数、 $y_i$ が $2^{t_i}$ の倍数であれば $\gcd(n,m)=\gcd(x_i,y_i)=\gcd(x_i,y_i/2^{t_i})=\gcd(x_{i+1},y_{i+1})$

$(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}$

\begin{pmatrix}s_{i+1}&a_{i+1}&c_{i+1}&x_{i+1}\\&b_{i+1}&d_{i+1}&y_{i+1}\end{pmatrix}= \begin{pmatrix}s_i+t_i&2^{t_i}a_i&2^{t_i}c_i&x_i\\&b_i&d_i&y_i/2^{t_i}\end{pmatrix}

式変形手順の具体例 ( $m$ が $3$ 以上の奇数の場合 )

同じ mod でのモジュラ逆数の計算を多数行う想定で調整した例です。実行環境、言語環境、問題の条件などによって最適な手順ではない場合もあるかと思います。

$\lbrace 2^0, 2^{-1}, 2^{-2}, ..., 2^{-63}\rbrace$ の $\operatorname{mod}m$ 事前計算を省いた手順は後述の「 $m$ が $2$ 以上の偶数の場合」を参照してください。

例えば $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)=y\ge 3$ のためモジュラ逆数が存在しない
$x\ge 2$ の偶数であれば $x$ を奇数にするような自然数 $t$ にて式変形C. を実行
$x=1$ であれば正常終了 : $\displaystyle n^{-1}\equiv 2^{-s}a\quad(\operatorname{mod}m)$ を計算して出力
$x=y$ であれば異常終了 : $\gcd(n,m)=x=y\ge 3$ のためモジュラ逆数が存在しない
$x\lt 8y$ であれば $t=\lfloor y/x\rfloor$ として式変形B. を実行
$x\lt y$ であれば $t=1$ として式変形B. を実行
$y=0$ であれば異常終了 : $\gcd(n,m)=x\ge 3$ のためモジュラ逆数が存在しない
$y\ge 2$ の偶数であれば $y$ を奇数にするような自然数 $t$ にて式変形D. を実行
$y=1$ であれば正常終了 : $\displaystyle n^{-1}\equiv -2^{-t}b\quad(\operatorname{mod}m)$ を計算して出力
$8x\gt y$ であれば $t=\lfloor x/y\rfloor$ として式変形A. を実行
$x\gt y$ であれば $t=1$ として式変形A. を実行
手順 3. に戻る

補題: $m$ が $3$ 以上の奇数の場合の $2$ の累乗数に対するモジュラ逆数の事前計算

$m$ が奇数 : $m\equiv 1\ (\operatorname{mod}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}$

例:

\begin{align*} 1 \times 2^{0} &\equiv 1 \mod 998244353\\ (1+998244353)/2 \times 2^{1} \equiv 499122177 \times 2^{1} &\equiv 1 \mod 998244353\\ (499122177+998244353)/2 \times 2^{2} \equiv 748683265 \times 2^{2} &\equiv 1 \mod 998244353\\ (748683265+998244353)/2 \times 2^{3} \equiv 873463809 \times 2^{3} &\equiv 1 \mod 998244353\\ (873463809+998244353)/2 \times 2^{4} \equiv 935854081 \times 2^{4} &\equiv 1 \mod 998244353\\ (935854081+998244353)/2 \times 2^{5} \equiv 967049217 \times 2^{5} &\equiv 1 \mod 998244353\\ (967049217+998244353)/2 \times 2^{6} \equiv 982646785 \times 2^{6} &\equiv 1 \mod 998244353\\ &\vdots\\ (436731897+998244353)/2 \times 2^{28} \equiv 717488125 \times 2^{28} &\equiv 1 \mod 998244353\\ (717488125+998244353)/2 \times 2^{29} \equiv 857866239 \times 2^{29} &\equiv 1 \mod 998244353\\ (857866239+998244353)/2 \times 2^{30} \equiv 928055296 \times 2^{30} &\equiv 1 \mod 998244353\\ (928055296)/2 \times 2^{31} \equiv 464027648 \times 2^{31} &\equiv 1 \mod 998244353\\ (464027648)/2 \times 2^{32} \equiv 232013824 \times 2^{32} &\equiv 1 \mod 998244353\\ (232013824)/2 \times 2^{33} \equiv 116006912 \times 2^{33} &\equiv 1 \mod 998244353\\ &\vdots\\ \end{align*}

式変形手順の具体例 ( $m$ が $2$ 以上の偶数の場合 )

$(a,b,x,y)=(1,0,n,m)$ と初期化 ( $c,d,s$ は使わない )
$x=0$ であれば異常終了 : $\gcd(n,m)=y\ge 2$ のためモジュラ逆数が存在しない
$x=1$ であれば正常終了 : $\displaystyle n^{-1}\equiv a\quad(\operatorname{mod}m)$ を計算して出力
$x\le y$ であれば $t=\lfloor y/x\rfloor$ として式変形B. を実行
$y=0$ であれば異常終了 : $\gcd(n,m)=x\ge 2$ のためモジュラ逆数が存在しない
$y=1$ であれば正常終了 : $\displaystyle n^{-1}\equiv -b\quad(\operatorname{mod}m)$ を計算して出力
$x\ge y$ であれば $t=\lfloor x/y\rfloor$ として式変形A. を実行
手順 2. に戻る

$m=1$ の場合

$nn^{-1} \equiv 0 \equiv 1\ (\operatorname{mod}1)$ より、 $n^{-1}\equiv 0\ (\operatorname{mod}1)$ を出力

$m=0$ の場合

未定義動作

サンプルコード (Rust)

// モジュラ演算用定数構造体
pub struct ModBin {
    modulus: u32,
    im: u64,
    ib: [u32; 64],
}
impl ModBin {
    pub fn new(modulus: u32) -> Self {
        assert_ne!(modulus, 0);
        // `im` $= \lceil (2^{64} / M) \rceil = \lfloor (2^{64}-1) / M\rfloor + 1$
        let im = (!0u64) / u64::from(modulus) + 1;
        let mut ib = [0u32; 64];
        ib[0] = 1;
        if modulus % 2 != 0 {
            // `modulus` が奇数の時は、 $\{2^0,2^{-1},...,2^{-63}\}$ を事前計算
            for i in 1..64 {
                let prev = ib[i - 1];
                if prev % 2 == 0 {
                    ib[i] = prev / 2;
                } else {
                    ib[i] = (prev / 2) + (modulus / 2) + 1; 
                }
            }
        }
        Self { modulus, im, ib }
    }
    pub fn add(&self, a: u32, b: u32) -> u32 {
        let v = u64::from(a) + u64::from(b);
        match v.overflowing_sub(u64::from(self.modulus)) {
            (_, true) => v as u32,
            (w, false) => w as u32,
        }
    }
    pub fn sub(&self, a: u32, b: u32) -> u32 {
        match a.overflowing_sub(b) {
            (v, true) => v.wrapping_add(self.modulus),
            (v, false) => v,
        }
    }
    pub fn neg(&self, a: u32) -> u32 {
        self.sub(0, a)
    }
    pub fn mul(&self, a: u32, b: u32) -> u32 {
        // Barrett Reduction
        let z = u64::from(a) * u64::from(b);
        let x = ((u128::from(z) * u128::from(self.im)) >> 64) as u64;
        match z.overflowing_sub(x.wrapping_mul(u64::from(self.modulus))) {
            (v, true) => (v as u32).wrapping_add(self.modulus),
            (v, false) => v as u32,
        }
    }
    pub fn div(&self, a: u32, b: u32) -> Option<u32> {
        self.inv_bingcd(b).map(|ib| self.mul(a, ib))
    }
    // モジュラ逆数の計算(バイナリGCD法無しの実装例)
    pub fn inv_plain(&self, n: u32) -> Option<u32> {
        let m = self.modulus;
        assert_ne!(m, 0);
        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) = 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) = x \ge 2$
                0 => return None,
                // $n^{-1} \equiv -b$
                1 => return Some(m - b),
                _ => {},
            }
            a += b * (x / y);
            x %= y;
        }
    }
    // モジュラ逆数の計算(modが奇数の場合、バイナリGCD法のハイブリッド)
    pub fn inv_bingcd(&self, n: u32) -> Option<u32> {
        let m = self.modulus;
        assert_ne!(m, 0);
        let (mut a, mut b, mut x, mut y) = (1, 0, n, m);
        if m == 1 {
            return Some(0);
        } else if m % 2 == 0 {
            // m が偶数: 拡張ユークリッド互除法
            loop {
                match x {
                    // $\gcd(n,m) = 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) = x \ge 2$
                    0 => return None,
                    // $n^{-1} \equiv -b$
                    1 => return Some(m - b),
                    _ => {},
                }
                a += b * (x / y);
                x %= y;
            }
        } else {
            // m が奇数: バイナリGCD + 拡張ユークリッド互除法
            if n == 0 {
                // $\gcd(n,m) = m \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}$
                    return Some(self.mul(a, self.ib[s as usize]));
                }
                if x == y {
                    // $\gcd(n,m) = x = y \ge 3$
                    return None;
                }
                if x < (y >> 3) {
                    // 式変形B: $t = \lfloor y / x \rfloor$
                    b += a * (y / x);
                    y %= x;
                    if y == 0 {
                        // $\gcd(n,m) = 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}$
                    return Some(self.mul(m - b, self.ib[s as usize]));
                }
                if (x >> 3) > y {
                    // 式変形A: $t = \lfloor x / y \rfloor$
                    a += b * (x / y);
                    x %= y;
                    if x == 0 {
                        // $\gcd(n,m) = 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;
                    }
                }
            }
        }
    }
}

$m$ が累乗数の場合のモジュラ逆数 (ニュートン・ラフソン法)

自然数 $p,q$ を用いて $m=p^q$ とできる場合、自然数 $r$ とニュートン・ラフソン法における、実数の逆数近似計算の漸化式 $x_{i+1}=x_i(2-nx_i)$ を用いて

$nx_i\equiv 1\ (\operatorname{mod}p^r),\quad x_{i+1}=x_i(2-nx_i) \quad\Longrightarrow\quad nx_{i+1}\equiv 1\ (\operatorname{mod}p^{2r})$

と計算する事ができます。また、 $q\le 2r$ であれば、

$q\le 2r$ かつ $nx_i\equiv 1\ (\operatorname{mod}p^r) \quad\Longrightarrow\quad nx_{i+1}\equiv 1\ (\operatorname{mod}p^q)$

となります。

証明:

$nx_i=1+kp^r$ とおくと、

$nx_{i+1}=nx_i(2-nx_i)=(1+kp^r)(2-(1+kp^r))$

$=(1+kp^r)(1-kp^r)=(1-k^2p^{2r})\equiv 1\ (\operatorname{mod}p^{2r})$

また、 1回の反復で有効桁数を2倍ではなく3倍にする方法もあります。

$nx_i\equiv 1\ (\operatorname{mod}p^r),\quad x_{i+1}=x_i(nx_i(nx_i-3)+3)\quad\Longrightarrow\quad nx_{i+1}\equiv 1\ (\operatorname{mod}p^{3r})$

証明:

$nx_i=1+kp^r$ とおくと、
$nx_{i+1}=nx_i(nx_i(nx_i-3)+3)$
$=1+(nx_i-1)^3=1+k^3p^{3r}\equiv 1\quad(\operatorname{mod}\ p^{3r})$

例: $m$ が $2$ の累乗数の場合

$n$ が奇数 : $n\ \operatorname{mod}2\equiv 1$ であれば、 $n^2\equiv 1\ (\operatorname{mod}8)$ となります。 $m=2^q$ とすると、

$x_0=n\quad\Longrightarrow\quad nx_0\equiv 1\ (\operatorname{mod}2^3)$
$x_1=x_0(2-nx_0)\mod 2^q\quad\Longrightarrow\quad nx_1\equiv 1\ (\operatorname{mod}2^{\min(6,q)})$
$x_2=x_1(2-nx_1)\mod 2^q\quad\Longrightarrow\quad nx_2\equiv 1\ (\operatorname{mod}2^{\min(12,q)})$
$x_3=x_2(2-nx_2)\mod 2^q\quad\Longrightarrow\quad nx_3\equiv 1\ (\operatorname{mod}2^{\min(24,q)})$
$x_4=x_3(2-nx_3)\mod 2^q\quad\Longrightarrow\quad nx_4\equiv 1\ (\operatorname{mod}2^{\min(48,q)})$
$x_5=x_4(2-nx_4)\mod 2^q\quad\Longrightarrow\quad nx_5\equiv 1\ (\operatorname{mod}2^{\min(96,q)})$

のように反復計算を必要な回数 ( $q \le 3\times 2^i$ ) 行えば、 $n^{-1}\equiv x_i\ (\operatorname{mod}m)$ を計算することができます。

また、以下のように計算する桁数を反復の初期でケチる事もできます。

$x_0=n\quad\Longrightarrow\quad nx_0\equiv 1\ (\operatorname{mod}2^3)$
$x_1=x_0(2-nx_0)\mod 2^4\quad\Longrightarrow\quad nx_1\equiv 1\ (\operatorname{mod}2^4)$
$x_2=x_1(2-nx_1)\mod 2^8\quad\Longrightarrow\quad nx_2\equiv 1\ (\operatorname{mod}2^8)$
$x_3=x_2(2-nx_2)\mod 2^{16}\quad\Longrightarrow\quad nx_3\equiv 1\ (\operatorname{mod}2^{16})$
$x_4=x_3(2-nx_3)\mod 2^{32}\quad\Longrightarrow\quad nx_4\equiv 1\ (\operatorname{mod}2^{32})$
$x_5=x_4(2-nx_4)\mod 2^{64}\quad\Longrightarrow\quad nx_5\equiv 1\ (\operatorname{mod}2^{64})$

サンプルコード (Rust)

// モジュラ逆数 $n^{-1} \mod 2^{32}$ の計算
pub fn invmod_u32(n: u32) -> Option<u32> {
    // n が偶数だとモジュラ逆数は存在しない
    if n % 2 == 0 {
        return None;
    }
    // 初期値は n 自身
    let mut x = n;
    for _ in 0..4 {
        // $x := x * (2 - (n * x)) \mod 2^{32}$ を4回繰り返す
        x = x.wrapping_mul(2u32.wrapping_sub(n.wrapping_mul(x)));
    }
    // できあがり。
    Some(x)
}
// モジュラ逆数 $n^{-1} \mod 2^{64}$ の計算
pub fn invmod_u64(n: u64) -> Option<u64> {
    // n が偶数だとモジュラ逆数は存在しない
    if n % 2 == 0 {
        return None;
    }
    // 初期値は n 自身
    let mut x = n;
    for _ in 0..5 {
        // $x := x * (2 - (n * x)) \mod 2^{64}$ を5回繰り返す
        x = x.wrapping_mul(2u64.wrapping_sub(n.wrapping_mul(x)));
    }
    // できあがり。
    Some(x)
}

Mizar/みざー 2023/01/03に更新

極力、（記述するプログラムコードの上では）符号なし整数型を用いた演算による拡張ユークリッドの互除法により、モジュラ逆数の計算を目指す。

互いに素な正の整数 $a,\ b$ について、 $ax-by=1$ となる整数 $0\lt x \lt b,\ 0\lt y \lt a$ の組を求める。特に、 $\gcd(a,b)=1$ かつ $0\lt a \lt b$ である時の、 $a\times x\operatorname{mod}b=1$ を満たす整数 $x$ （モジュラ逆数）を高速に求める事を目指す。

例として

$100^{-1}\operatorname{\ mod\ }998244353$ を求める
$100x-998244353y=1$ を満たす整数 $x,y$ の組において、 $x$ が最小の正の整数となるものを求める

拡張ユークリッドの互除法とベズーの等式と一次不定方程式:

https://manabitimes.jp/math/674
https://manabitimes.jp/math/672

この項目では、拡張ユークリッドの互除法における係数の操作の説明として、行列の積による表現を用いている。

https://ja.wikipedia.org/wiki/行列
https://manabitimes.jp/math/1023

除算・剰余算による拡張ユークリッドの互除法

演算回数の最悪ケース: $O(\log_\tau \min(a,b))$ , $\tau = \frac{1+\sqrt{5}}{2}$ (黄金数), $(5\lfloor\log_{10} 10\min(a,b)\rfloor)$ 回以下(ラメの定理)

フィボナッチ数列 (OEIS000045, $F_{0}$ ～ $F_{2000}$ の数表) の隣接項、例えば $F_{43}=433494437, F_{44}=701408733$ を係数に用いた、 $433494437x-701408733y=1$ の拡張ユークリッド互除法

https://oeis.org/A000045

https://oeis.org/A000045/b000045.txt

$ax-by=\gcd(a,b)$ での拡張ユークリッドの互除法では、縦ベクトル $\begin{pmatrix}a\\-b\end{pmatrix}$ に対して以下のような操作を行う。

縦ベクトルの絶対値が小さな方が0ではない間、絶対値の大きな方を、小さな方で割り算し、割り算の商を用いた操作によって割り算の余りに置き換える操作を繰り返す。
縦ベクトルの絶対値が小さな方が0になったら、必要な補正を行った上で出力を行う。縦ベクトルの絶対値が0でない方 $=\gcd(a,b)$

\begin{align*} \text{初期値}\\ \begin{pmatrix}1&0\\0&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}100\\-998244353\end{pmatrix}\\ \text{2行目を1行目で割る}&\text{商:9982443 余り:53}\\ \begin{pmatrix}1&0\\9982443&1\end{pmatrix} \begin{pmatrix}1&0\\0&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&0\\9982443&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}\\ \begin{pmatrix}1&0\\9982443&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}100\\-53\end{pmatrix}\\ \text{1行目を2行目で割る}&\text{商:1 余り:47}\\ \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}1&0\\9982443&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}100\\-53\end{pmatrix}\\ \begin{pmatrix}9982444&1\\9982443&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}47\\-53\end{pmatrix}\\ \text{2行目を1行目で割る}&\text{商:1 余り:6}\\ \begin{pmatrix}1&0\\1&1\end{pmatrix} \begin{pmatrix}9982444&1\\9982443&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&0\\1&1\end{pmatrix} \begin{pmatrix}47\\-53\end{pmatrix}\\ \begin{pmatrix}9982444&1\\19964887&2\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}47\\-6\end{pmatrix}\\ \text{1行目を2行目で割る}&\text{商:7 余り:5}\\ \begin{pmatrix}1&7\\0&1\end{pmatrix} \begin{pmatrix}9982444&1\\19964887&2\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&7\\0&1\end{pmatrix} \begin{pmatrix}47\\-6\end{pmatrix}\\ \begin{pmatrix}149736653&15\\19964887&2\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}5\\-6\end{pmatrix}\\ \text{2行目を1行目で割る}&\text{商:1 余り:1}\\ \begin{pmatrix}1&0\\1&1\end{pmatrix} \begin{pmatrix}149736653&15\\19964887&2\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&0\\1&1\end{pmatrix} \begin{pmatrix}5\\-6\end{pmatrix}\\ \begin{pmatrix}149736653&15\\169701540&17\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}5\\-1\end{pmatrix}\\ \text{1行目を2行目で割る}&\text{商:5 余り:0}\\ \begin{pmatrix}1&5\\0&1\end{pmatrix} \begin{pmatrix}149736653&15\\169701540&17\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&5\\0&1\end{pmatrix} \begin{pmatrix}5\\-1\end{pmatrix}\\ \begin{pmatrix}998244353&100\\169701540&17\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}0\\-1\end{pmatrix}\\ \text{右辺で}{-1}\text{の方が残ってしまったので、}\\ \text{右辺が}{1}\text{の時の値を求めるために}\\ \text{1行目から2行目を引くと}\\ \begin{pmatrix}1&-1\\0&1\end{pmatrix} \begin{pmatrix}998244353&100\\169701540&17\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&-1\\0&1\end{pmatrix} \begin{pmatrix}0\\-1\end{pmatrix}\\ \begin{pmatrix}828542813&83\\169701540&17\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1\\-1\end{pmatrix}\\ \text{よって}\\ 828542813 \times 100 - 83 \times 998244353 &= 1\\ 169701540 \times 100 - 17 \times 998244353 &= -1\\ \end{align*}

$x=828542813, y=83$ の時、 $100x-998244353y=1$ を満たす。

Mizar/みざー 2022/12/13に更新

除算・剰余算による拡張ユークリッドの互除法のPythonコード例

def extgcd(a, b):
  x, y, z, w = 1, 0, 0, 1
  while True:
    if a == 0: return x - y, z - w, b
    y, w, b = y + (b // a) * x, w + (b // a) * z, b % a
    if b == 0: return x, z, a
    x, z, a = x + (a // b) * y, z + (a // b) * w, a % b

import math
m = 998244353
for i in range(105):
  # i * x - m * y = g; 0 <= x <= m/g; -1 <= y <= i/g; g = gcd(i, m)
  # if i == 0 then y == -1, else y >= 0
  x, y, g = extgcd(i, m)
  print(i, x, y, g, i * x % m)

0 1 -1 998244353 0
1 1 0 1 1
2 499122177 1 1 1
3 332748118 1 1 1
4 748683265 3 1 1
5 598946612 3 1 1
6 166374059 1 1 1
7 855638017 6 1 1
8 873463809 7 1 1
9 443664157 4 1 1
10 299473306 3 1 1
11 272248460 3 1 1
12 582309206 7 1 1
13 460728163 6 1 1
14 926941185 13 1 1
15 865145106 13 1 1
16 935854081 15 1 1
17 939524097 16 1 1
18 720954255 13 1 1
19 105078353 2 1 1
20 149736653 3 1 1
21 617960790 13 1 1
22 136124230 3 1 1
23 217009642 5 1 1
24 291154603 7 1 1
25 319438193 8 1 1
26 729486258 19 1 1
27 480636170 13 1 1
28 962592769 27 1 1
29 481911067 14 1 1
30 432572553 13 1 1
31 128805723 4 1 1
32 967049217 31 1 1
33 756245722 25 1 1
34 968884225 33 1 1
35 370776474 13 1 1
36 859599304 31 1 1
37 242816194 9 1 1
38 551661353 21 1 1
39 486324172 19 1 1
40 573990503 23 1 1
41 97389693 4 1 1
42 308980395 13 1 1
43 510729669 22 1 1
44 68062115 3 1 1
45 288381702 13 1 1
46 108504821 5 1 1
47 191153174 9 1 1
48 644699478 31 1 1
49 692659347 34 1 1
50 658841273 33 1 1
51 313174699 16 1 1
52 364743129 19 1 1
53 828731161 44 1 1
54 240318085 13 1 1
55 54449692 3 1 1
56 980418561 55 1 1
57 700522353 40 1 1
58 740077710 43 1 1
59 727534020 43 1 1
60 715408453 43 1 1
61 654586461 40 1 1
62 563525038 35 1 1
63 205986930 13 1 1
64 982646785 63 1 1
65 890741115 58 1 1
66 378122861 25 1 1
67 432075914 29 1 1
68 983564289 67 1 1
69 405084665 28 1 1
70 185388237 13 1 1
71 702988981 50 1 1
72 429799652 31 1 1
73 382888245 28 1 1
74 121408097 9 1 1
75 771975633 58 1 1
76 774952853 59 1 1
77 894530654 69 1 1
78 243162086 19 1 1
79 404352143 32 1 1
80 786117428 63 1 1
81 825708292 67 1 1
82 547817023 45 1 1
83 60135202 5 1 1
84 653612374 55 1 1
85 387553690 33 1 1
86 754487011 65 1 1
87 493385140 43 1 1
88 533153234 47 1 1
89 392568004 35 1 1
90 144190851 13 1 1
91 65818309 6 1 1
92 553374587 51 1 1
93 42935241 4 1 1
94 95576587 9 1 1
95 819611153 78 1 1
96 322349739 31 1 1
97 699800165 68 1 1
98 845451850 83 1 1
99 584830025 58 1 1
100 828542813 83 1 1
101 889524671 90 1 1
102 655709526 67 1 1
103 38766771 4 1 1
104 681493741 71 1 1

$a=0$ の時、この実装でのextgcd関数が負の値を返している（モジュラ逆数と最大公約数だけを求めたい場合は計算を省略しても良い部分）。符号なし整数型でこの手法による拡張ユークリッドの互除法を実装する際には、この $a=0$ のケースについて実装上の注意が必要。

Mizar/みざー 2022/12/16

拡張ユークリッドの互除法による乗法逆元（除算・剰余算）

from typing import Tuple, Optional
def mod_inv_extgcd(n: int, m: int) -> Tuple[Optional[int], int]:
    # extended euclidean algorithm
    # return (n^{-1} mod m) or None
    x, y, a, b = 1, 0, n, m
    assert a >= 0 and b >= 0
    while True:
        if a == 0:
            if b == 1: # b = gcd(n,m) = 1
                return x - y, b
            else: # b = gcd(n,m) != 1
                return None, b
        y, b = y + (b // a) * x, b % a
        if b == 0:
            if a == 1: # a = gcd(n,m) = 1
                return x, a
            else: # a = gcd(n,m) != 1
                return None, a
        x, a = x + (a // b) * y, a % b

m = 998244353
for i in range(105):
  inv, g = mod_inv_extgcd(i, m)
  if inv is None:
    print(i, g, inv)
  else:
    print(i, g, inv, i * inv % m)

拡張ユークリッドの互除法による乗法逆元（バイナリ法1）

from typing import Tuple, Optional
def mod_inv_binextgcd(n: int, m: int) -> Tuple[Optional[int], int]:
    x, y, a, b, s = 1, 0, n, m, 0
    # n * x mod m = 1 となる 整数x を返す
    # mが偶数である場合、mが負である場合、nが負である場合は動作未定義
    assert (m & 1) == 1 and m > 0 and n >= 0
    if a == 0:
        return None, b
    while (a & 1) == 0:
        a, s = (a>>1), (s+1)
    while True:
        if a == 1: # a = gcd(n,m) = 1
            # n^{-1} mod m = x * 2^{-s} mod m
            for _ in range(s):
                if (x & 1) == 0:
                    x = (x >> 1)
                else:
                    x = ((x + m) >> 1)
            return x, a
        if a == b: # a = b = gcd(n,m) != 1
            return None, a
        while a < b:
            y, b = (x+y), (b-a)
            while (b & 1) == 0:
                x, b, s = (x<<1), (b>>1), (s+1)
        if b == 1: # b = gcd(n,m) = 1
            # n^{-1} mod m = -y * 2^{-s} mod m
            for _ in range(s):
                if (y & 1) == 0:
                    y = (y >> 1)
                else:
                    y = ((y + m) >> 1)
            return m - y, b
        while a > b:
            x, a = (x+y), (a-b)
            while (a & 1) == 0:
                y, a, s = (y<<1), (a>>1), (s+1)

m = 998244353
for i in range(105):
  inv, g = mod_inv_binextgcd(i)
  if inv is None:
    print(i, g, inv)
  else:
    print(i, g, inv, i * inv % m)

拡張ユークリッドの互除法による乗法逆元（バイナリ法2）

$\left[ 2^0, 2^{-1}, 2^{-2}, \dots \right] \mod m$ を事前計算する版

from typing import Tuple, Optional
class OddModulo:
    def __init__(self, m: int) -> None:
        self.m = m
        # must be odd modulo
        assert m & 1 == 1 and m > 0
        # b = [2^{0} mod m, 2^{-1} mod m, ...]
        b = [1]
        for _ in range(m.bit_length() * 2):
            if (b[-1] & 1) == 0:
                b.append(b[-1] // 2)
            else:
                b.append((b[-1] + m) // 2)
        self.b = b
    def mod_inv_binextgcd(self, n: int) -> Tuple[Optional[int], int]:
        # binary extended euclidean algorithm
        x, y, a, b, s = 1, 0, n, self.m, 0
        assert a >= 0
        if a == 0:
            return None, b
        while (a & 1) == 0:
            a, s = (a>>1), (s+1)
        if a == 1: # a = gcd(n,m) = 1
            # n^{-1} mod m = x * 2**-s mod m
            return x * self.b[s] % self.m, a
        if a < (b >> 9):
            # b が a に比べてかなり大きい場合は
            # 除算・剰余算を使うメリットのが大きいので使っておく
            y, b = (y + x * (b // a)), (b % a)
            if b == 0:
                return None, a
            while (b & 1) == 0:
                x, b, s = (x << 1), (b >> 1), (s + 1)
        while True:
            if a == 1: # a = gcd(n,m) = 1
                # n^{-1} mod m = x * 2**-s mod m
                return x * self.b[s] % self.m, a
            if a == b: # a = b = gcd(n,m) != 1
                # modular inverse does not exist.
                return None, a
            while a < b:
                y, b = (y + x), (b - a)
                while (b & 1) == 0:
                    x, b, s = (x << 1), (b >> 1), (s + 1)
            if b == 1: # b = gcd(n,m) = 1
                # n^{-1} mod m = -y * 2^{-s} mod m
                return self.m - y * self.b[s] % self.m, b
            while a > b:
                x, a = (x + y), (a - b)
                while (a & 1) == 0:
                    y, a, s = (y << 1), (a >> 1), (s + 1)

m = 998244353
modulo = OddModulo(m)
for i in range(105):
  inv, g = modulo.mod_inv_binextgcd(i)
  if inv is None:
    print(i, g, inv)
  else:
    print(i, g, inv, i * inv % m)

Mizar/みざー 2023/01/11に更新

Python 3.8 以降の場合、モジュラ逆数の計算は組み込み関数 pow(base, exp[, mod]) で出来ます。
(AtCoderで使われている PyPy3 7.3.0 だと Python 3.6 相当のためまだ対応していない)

Python (3.8.2) import sys;print(sys.version)

3.8.2 (default, Feb 26 2020, 02:56:10) 
[GCC 7.4.0]

PyPy3 (7.3.0) import sys;print(sys.version)

3.6.9 (7.3.0+dfsg-1~ppa1~ubuntu18.04, Dec 24 2019, 08:12:19)
[PyPy 7.3.0 with GCC 7.4.0]

# Python 3.8以降
n = 100
m = 998244353
inv = pow(n, -1, m)
print(n, m, inv, n * inv % m)

100 998244353 828542813 1

Mizar/みざー 2022/12/17に更新

フェルマーの小定理

https://manabitimes.jp/math/680

$n$ が整数で $m$ が素数の場合、 $n^{m} \equiv n \mod m$
$n$ が整数で $m$ が素数で $\gcd(n,m)=1$ の場合、 $n^{m-1} \equiv 1 \mod m$
$n$ が整数で $m$ が素数で $\gcd(n,m)=1$ の場合、 $n^{m-2} \equiv n^{-1} \mod m$

# Python 3.8以降
n = 100
m = 998244353
inv = pow(n, m-2, m)
print(n, m, inv, n * inv % m)

100 998244353 828542813 1

Mizar/みざー 2023/01/20に更新

$A=\left(\begin{array}{c:c}a&b\\\hdashline c&d\end{array}\right)$ , $B=\left(\begin{array}{c:c}p&q\\\hdashline r&s\end{array}\right)$ , $\vec{x}=\left(\begin{array}{c}t\\\hdashline u\end{array}\right)$ の時の、
行列の積の結合法則 $(AB)\vec{x}=A(B\vec{x})$ とそれぞれの積の計算の仕方の説明

行列・縦ベクトルの外側にある矢印は、行列の積の計算時に要素を見る順番の方向を表したもの（行列の積の説明用であり、通常はこの矢印は書かない）
行列の内側にある点線は、行列の要素の区切りを示したもの（説明用であり、通常はあまりこの区切り線は書かない）
「1行2列の行ベクトル $A$ 」と「2行2列の行列 $B$ 」の積 → 「1行2列の行ベクトル」（「 $j$ 列目の結果」は「 $A$ の列ベクトル」と「 $B$ の＄j$列目の列ベクトル」との内積）
「2行2列の行列 $A$ 」と「2行1列の列ベクトル $B$ 」の積 → 「2行1列の列ベクトル」（「 $i$ 行目の結果」は「 $A$ の $i$ 行目の行ベクトル」と「 $B$ の列ベクトル」との内積）
「2行2列の行列 $A$ 」と「2行2列の行列 $B$ 」の積 → 「2行2列の行列」（「 $i$ 行 $j$ 列目の結果」は「 $A$ の $i$ 行目の行ベクトル」と「 $B$ の $j$ 列目の列ベクトル」との内積）

\begin{align*} (AB)\vec{x}&=A(B\vec{x})\\ \underrightarrow{\overrightarrow{\left(\begin{array}{c:c}a&b\\\hdashline c&d\end{array}\right)}} \left\downarrow\left(\begin{array}{c:c}p&q\\\hdashline r&s\end{array}\right)\right\downarrow \left(\begin{array}{c}t\\\hdashline u\end{array}\right)&= \left(\begin{array}{c:c}a&b\\\hdashline c&d\end{array}\right) \underrightarrow{\overrightarrow{\left(\begin{array}{c:c}p&q\\\hdashline r&s\end{array}\right)}} \left\downarrow\left(\begin{array}{c}t\\\hdashline u\end{array}\right)\right\downarrow\\ \underrightarrow{\overrightarrow{\left(\begin{array}{c:c}ap+br&aq+bs\\\hdashline cp+dr&cq+ds\end{array}\right)}} \left\downarrow\left(\begin{array}{c}t\\\hdashline u\end{array}\right)\right\downarrow&= \underrightarrow{\overrightarrow{\left(\begin{array}{c:c}a&b\\\hdashline c&d\end{array}\right)}} \left\downarrow\left(\begin{array}{c}pt+qu\\\hdashline rt+su\end{array}\right)\right\downarrow\\ \left(\begin{array}{c}(ap+br)t+(aq+bs)u\\\hdashline(cp+dr)t+(cq+ds)u\end{array}\right) &= \left(\begin{array}{c}a(pt+qu)+b(rt+su)\\\hdashline c(pt+qu)+d(rt+su)\end{array}\right) \\ \left(\begin{array}{c}apt+aqu+brt+bsu\\\hdashline cpt+cqu+drt+dsu\end{array}\right)&=\left(\begin{array}{c}apt+aqu+brt+bsu\\\hdashline cpt+cqu+drt+dsu\end{array}\right) \end{align*}

Mizar/みざー 2023/01/03に更新

バイナリGCDアルゴリズム（主に加減算・シフト演算）によるモジュラ逆数

演算回数の最悪ケース?例: $O(2\log_2 ab)$

$12297829382473034410x-12297829382473034411y=1$ のバイナリGCDアルゴリズム
$(2^{64}-2^{63})x-(2^{64}-1)y=1$ のバイナリGCDアルゴリズム
$(2^{64}-2^{62})x-(2^{64}-3)y=1$ のバイナリGCDアルゴリズム
もっと回数が多いケースが存在するかもしれないが、行列操作で $2\log_2 ab$ 回を上回ることはない（絶対値を減らすように行う減算と2の累乗での除算を交互に行うので、この回数を上回る事ができない）

$ax-by=\gcd(a,b)$ のバイナリGCDアルゴリズムでは、縦ベクトル $\begin{pmatrix}a\\-b\end{pmatrix}$ に対して、以下のような操作・動作を行う。

片方、もしくは両方が $0$ であれば、結果を出力する
片方の絶対値が $1$ になるか、双方が同じ絶対値になるまでの間、次の2つの操作（偶数を奇数にする操作・引き算する操作）を繰り返す
a. 偶数があれば、それを $2$ の累乗で割って奇数にする操作
b. ともに奇数であれば、絶対値の大きな方をもう一方の数で引き算する操作
片方の絶対値が $1$ になる ( $\gcd(a,b)/\gcd(a,b,2^s)=1$ ) か、双方が同じ絶対値 ( $=\gcd(a,b)/\gcd(a,b,2^s)$ ) になった場合は、 $2$ の累乗で割った操作の補正を必要に応じて行った上で、結果を出力する

非負整数 $a$ および正奇数 $m$ において、 $2h\equiv a\ (\operatorname{mod}m)$ を満たす非負整数 $h$ はこのように求められる:
$h=\begin{cases}a/2&\left(a\equiv 0 \left(\operatorname{mod}2\right)\right)\\(a+m)/2&\left(a\equiv 1 \left(\operatorname{mod}2\right)\right)\\\end{cases}$
これを繰り返し適用する事で、 $116006912 \times 2^{33} \equiv 1 \mod 998244353$ を求める事ができる。

1/2^n を事前計算する場合の最終結果の求め方

\begin{align*} 1 \times 2^{0} &\equiv 1 \mod 998244353\\ (1+998244353)/2 \times 2^{1} \equiv 499122177 \times 2^{1} &\equiv 1 \mod 998244353\\ (499122177+998244353)/2 \times 2^{2} \equiv 748683265 \times 2^{2} &\equiv 1 \mod 998244353\\ (748683265+998244353)/2 \times 2^{3} \equiv 873463809 \times 2^{3} &\equiv 1 \mod 998244353\\ (873463809+998244353)/2 \times 2^{4} \equiv 935854081 \times 2^{4} &\equiv 1 \mod 998244353\\ (935854081+998244353)/2 \times 2^{5} \equiv 967049217 \times 2^{5} &\equiv 1 \mod 998244353\\ (967049217+998244353)/2 \times 2^{6} \equiv 982646785 \times 2^{6} &\equiv 1 \mod 998244353\\ &\vdots\\ (436731897+998244353)/2 \times 2^{28} \equiv 717488125 \times 2^{28} &\equiv 1 \mod 998244353\\ (717488125+998244353)/2 \times 2^{29} \equiv 857866239 \times 2^{29} &\equiv 1 \mod 998244353\\ (857866239+998244353)/2 \times 2^{30} \equiv 928055296 \times 2^{30} &\equiv 1 \mod 998244353\\ (928055296)/2 \times 2^{31} \equiv 464027648 \times 2^{31} &\equiv 1 \mod 998244353\\ (464027648)/2 \times 2^{32} \equiv 232013824 \times 2^{32} &\equiv 1 \mod 998244353\\ (232013824)/2 \times 2^{33} \equiv 116006912 \times 2^{33} &\equiv 1 \mod 998244353\\ \end{align*}

\begin{align*} \text{さらに}\\ 445267313 \times 100 &\equiv 2^{33} \mod 998244353\\ 116006912 \times 2^{33} &\equiv 1 \mod 998244353\text{ より}\\ 116006912 \times 445267313 \times 100 &\equiv 1 \mod 998244353\\ 116006912 \times 445267313 &\equiv 828542813 \mod 998244353\text{ より}\\ 828542813 \times 100 &\equiv 1 \mod 998244353\\ \end{align*}

$116006912 \times 445267313 \equiv 828542813 \mod 998244353$ より、 $x=828542813$ の時、 $100x-998244353y=1$ を満たす整数 $y$ が存在する。

x/2^n を事後計算する場合の最終結果の求め方

\begin{align*} 445267313 \times 100 &\equiv 2^{33} \mod 998244353\\ (445267313 + 998244353)/2 \times 100 \equiv 721755833 \times 100 &\equiv 2^{32} \mod 998244353\\ (721755833 + 998244353)/2 \times 100 \equiv 860000093 \times 100 &\equiv 2^{31} \mod 998244353\\ (860000093 + 998244353)/2 \times 100 \equiv 929122223 \times 100 &\equiv 2^{30} \mod 998244353\\ (929122223 + 998244353)/2 \times 100 \equiv 963683288 \times 100 &\equiv 2^{29} \mod 998244353\\ &\vdots\\ (559016838)/2 \times 100 \equiv 279508419 \times 100 &\equiv 2^{4} \mod 998244353\\ (279508419 + 998244353)/2 \times 100 \equiv 638876386 \times 100 &\equiv 2^{3} \mod 998244353\\ (638876386)/2 \times 100 \equiv 319438193 \times 100 &\equiv 2^{2} \mod 998244353\\ (319438193 + 998244353)/2 \times 100 \equiv 658841273 \times 100 &\equiv 2^{1} \mod 998244353\\ (658841273 + 998244353)/2 \times 100 \equiv 828542813 \times 100 &\equiv 2^{0} \equiv 1 \mod 998244353\\ \end{align*}

Mizar/みざー 2023/01/03に更新

除算による拡張ユークリッド互除法とBinaryGCD法の混合によるモジュラ逆数の計算

\begin{align*} \text{初期値}\\ \begin{pmatrix}1&0\\0&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}100\\-998244353\end{pmatrix}\\ \text{1行目を }2^2\text{ で割る}\\ \frac{1}{2^2} \begin{pmatrix}1&0\\0&2^2\end{pmatrix} \begin{pmatrix}1&0\\0&1\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \frac{1}{2^2} \begin{pmatrix}1&0\\0&2^2\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}\\ \frac{1}{2^2} \begin{pmatrix}1&0\\0&4\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}25\\-998244353\end{pmatrix}\\ \text{2行目を1行目で割る}&\text{商:39929774 余り:3}\\ \frac{1}{2^2} \begin{pmatrix}1&0\\39929774&1\end{pmatrix} \begin{pmatrix}1&0\\0&4\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&0\\39929774&1\end{pmatrix} \begin{pmatrix}25\\-998244353\end{pmatrix}\\ \frac{1}{2^2} \begin{pmatrix}1&0\\39929774&4\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}25\\-3\end{pmatrix}\\ \text{1行目に2行目の値を足す}\\ \frac{1}{2^2} \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}1&0\\39929774&4\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}25\\-3\end{pmatrix}\\ \frac{1}{2^2} \begin{pmatrix}39929775&4\\39929774&4\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}22\\-3\end{pmatrix}\\ \text{1行目を }2\text{ で割る}\\ \frac{1}{2^2} \frac{1}{2} \begin{pmatrix}1&0\\0&2\end{pmatrix} \begin{pmatrix}39929775&4\\39929774&4\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \frac{1}{2} \begin{pmatrix}1&0\\0&2\end{pmatrix} \begin{pmatrix}22\\-3\end{pmatrix}\\ \frac{1}{2^3} \begin{pmatrix}39929775&4\\79859548&8\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}11\\-3\end{pmatrix}\\ \text{1行目に2行目の値を足す}\\ \frac{1}{2^3} \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}39929775&4\\79859548&8\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}11\\-3\end{pmatrix}\\ \frac{1}{2^3} \begin{pmatrix}119789323&12\\79859548&8\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}8\\-3\end{pmatrix}\\ \text{1行目を }2^3\text{ で割る}\\ \frac{1}{2^3} \frac{1}{2^3} \begin{pmatrix}1&0\\0&2^3\end{pmatrix} \begin{pmatrix}119789323&12\\79859548&8\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \frac{1}{2^3} \begin{pmatrix}1&0\\0&2^3\end{pmatrix} \begin{pmatrix}8\\-3\end{pmatrix}\\ \frac{1}{2^6} \begin{pmatrix}119789323&12\\638876384&64\end{pmatrix} \begin{pmatrix}100\\-998244353\end{pmatrix}&= \begin{pmatrix}1\\-3\end{pmatrix}\\ \text{よって}\\ 119789323 \times 100 - 12 \times 998244353 &= 2^{6}\\ 638876384 \times 100 - 64 \times 998244353 &= -3 \times 2^{6}\\ \end{align*}

1/2^n を事前計算する場合の最終結果の求め方

\begin{align*} 1 \times 2^0 &\equiv 1 \mod 998244353\\ (1+998244353)/2 \times 2^1 \equiv 499122177 \times 2^1 &\equiv 1 \mod 998244353\\ (499122177+998244353)/2 \times 2^2 \equiv 748683265 \times 2^2 &\equiv 1 \mod 998244353\\ (748683265+998244353)/2 \times 2^3 \equiv 873463809 \times 2^3 &\equiv 1 \mod 998244353\\ (873463809+998244353)/2 \times 2^4 \equiv 935854081 \times 2^4 &\equiv 1 \mod 998244353\\ (935854081+998244353)/2 \times 2^5 \equiv 967049217 \times 2^5 &\equiv 1 \mod 998244353\\ (967049217+998244353)/2 \times 2^6 \equiv 982646785 \times 2^6 &\equiv 1 \mod 998244353\\ \end{align*}

\begin{align*} \text{さらに}\\ 119789323 \times 100 &\equiv 2^{6} \mod 998244353\\ 982646785 \times 2^{6} &\equiv 1 \mod 998244353\text{ より}\\ 982646785 \times 119789323 \times 100 &\equiv 1 \mod 998244353\\ 982646785 \times 119789323 &\equiv 828542813 \mod 998244353\text{ より}\\ 828542813 \times 100 &\equiv 1 \mod 998244353\\ \end{align*}

x/2^n を事後計算する場合の最終結果の求め方

\begin{align*} 119789323 \times 100 &\equiv 2^6 \mod 998244353\\ (119789323 +998244353)/2 \times 100 \equiv 559016838 \times 100 &\equiv 2^5 \mod 998244353\\ (559016838 +998244353)/2 \times 100 \equiv 279508419 \times 100 &\equiv 2^4 \mod 998244353\\ (279508419 +998244353)/2 \times 100 \equiv 638876386 \times 100 &\equiv 2^3 \mod 998244353\\ (638876386)/2 \times 100 \equiv 319438193 \times 100 &\equiv 2^2 \mod 998244353\\ (319438193 +998244353)/2 \times 100 \equiv 658841273 \times 100 &\equiv 2^1 \mod 998244353\\ (658841273 +998244353)/2 \times 100 \equiv 828542813 \times 100 &\equiv 2^0 \equiv 1 \mod 998244353\\ \end{align*}

Mizar/みざー 2023/01/03

バイナリGCDアルゴリズムの変形（基数3）

$ax-by=\gcd(a,b)$ のバイナリGCDアルゴリズム（基数3）では、縦ベクトル $\begin{pmatrix}a\\-b\end{pmatrix}$ に対して、以下のような操作・動作を行う。

片方、もしくは両方が $0$ であれば、結果を出力する
片方の絶対値が $1$ になるか、双方が同じ絶対値になるまでの間、次の2つの操作（偶数を奇数にする操作・引き算する操作）を繰り返す
a. 3の倍数があれば、それを $3$ の累乗で割って 3の倍数でない数にする操作
b. ともに3の倍数でなければ、絶対値の大きな方をもう一方の数で足し算もしくは引き算して、3の倍数にする操作
片方の絶対値が $1$ になる ( $\gcd(a,b)/\gcd(a,b,3^s)=1$ ) か、双方が同じ絶対値 ( $=\gcd(a,b)/\gcd(a,b,3^s)$ ) になった場合は、 $3$ の累乗で割った操作の補正を必要に応じて行った上で、結果を出力する

モジュラ逆数の計算

m が素数、もしくは m の素因数分解が既知である場合のモジュラ逆数 (オイラーの定理)

オイラーのトーシェント関数

サンプルコード (Rust)

一般の m に関するモジュラ逆数 (変則・拡張ユークリッド互除法)

対象

目標

制約

初期値

式変形パターン

式変形手順の具体例 ( m が 3 以上の奇数 の場合 )

補題: m が 3 以上の奇数の場合の 2 の累乗数に対するモジュラ逆数の事前計算

式変形手順の具体例 ( m が 2 以上の偶数 の場合 )

m=1 の場合

m=0 の場合

サンプルコード (Rust)

m が累乗数の場合のモジュラ逆数 (ニュートン・ラフソン法)

例: m が 2 の累乗数の場合

サンプルコード (Rust)

除算・剰余算による拡張ユークリッドの互除法

除算・剰余算による拡張ユークリッドの互除法のPythonコード例

拡張ユークリッドの互除法による乗法逆元（除算・剰余算）

拡張ユークリッドの互除法による乗法逆元（バイナリ法1）

拡張ユークリッドの互除法による乗法逆元（バイナリ法2）

フェルマーの小定理

バイナリGCDアルゴリズム（主に加減算・シフト演算）によるモジュラ逆数

1/2^n を事前計算する場合の最終結果の求め方

x/2^n を事後計算する場合の最終結果の求め方

除算による拡張ユークリッド互除法とBinaryGCD法の混合によるモジュラ逆数の計算

1/2^n を事前計算する場合の最終結果の求め方

x/2^n を事後計算する場合の最終結果の求め方

バイナリGCDアルゴリズムの変形（基数3）

$m$ が素数、もしくは $m$ の素因数分解が既知である場合のモジュラ逆数 (オイラーの定理)

一般の $m$ に関するモジュラ逆数 (変則・拡張ユークリッド互除法)

式変形手順の具体例 ( $m$ が $3$ 以上の奇数の場合 )

補題: $m$ が $3$ 以上の奇数の場合の $2$ の累乗数に対するモジュラ逆数の事前計算

式変形手順の具体例 ( $m$ が $2$ 以上の偶数の場合 )

$m=1$ の場合

$m=0$ の場合

$m$ が累乗数の場合のモジュラ逆数 (ニュートン・ラフソン法)

例: $m$ が $2$ の累乗数の場合