background

同余方程学习笔记

解各种同余方程是同余问题的一个重要部分,本文介绍各种同余方程的求解方法。

Table of Contents

二元线性不定方程

形式

不定方程的范围很广泛,只要没有确定的解的方程都可以叫不定方程。本文主要讨论数论中的不定方程。

ax+by=c (a,b,c,x,y∈Z) ax+by=c\ (a,b,c,x,y\in \mathrm{Z})

形如上面这种形式,满足各项的系数和解均为整数的方程就是数论中的不定方程。

求解

Theorem 1:∀a,b∈Z\forall a, b \in \mathrm{Z},一定存在 x,y∈Zx,y\in \mathrm{Z},使得 ax+by=gcd⁡(a,b)ax+by=\gcd(a,b)。

这个定理又叫做裴蜀定理,根据这个定理,可以得出不定方程 ax+by=cax+by=c 有解的充要条件就是 gcd⁡(a,b)∣c\gcd(a,b)\mid c。记 d=gcd⁡(a,b)d=\gcd(a,b),则我们可以先求出 ax′+by′=dax'+by'=d 的解,再令 x=cdx′,y=cdy′x=\frac{c}{d}x',y=\frac{c}{d}y',就得出原方程的一组解。

求 x′,y′x',y' 的过程可以用扩展欧几里得算法,即 exgcd 解决。

Theorem 2:若 gcd⁡(a,b)=1\gcd(a,b)=1,且 x0,y0x_0,y_0 是方程 ax+by=cax+by=c 的一组解,则该方程的通解为 x=x0+kb,y=y0−kax=x_0+kb,y=y_0-ka,其中 k∈Zk\in \mathrm{Z}。

换句话说,x,yx,y 这两个解是具有周期性的。根据这个定理,我们知道一定有一个 x∈[0,b)x\in [0,b),而对于方程 ax+by=cax+by=c,则一定有一个 x∈[0,bgcd⁡(a,b))x\in [0,\frac{b}{\gcd(a,b)}),由此可以求一个不定方程的一个元的最小非负整数解,具体地,假设知道了一个解 x0x_0,则最小非负整数 x=(x0 mod m+m) mod mx=(x_0 \bmod m + m) \bmod m,其中 m=bgcd⁡(a,b)m=\frac{b}{\gcd(a,b)}。

总结一下,求解二元线性不定方程 ax+by=cax+by=c 的流程如下:

  1. 求出 ax+by=gcd⁡(a,b)ax+by=\gcd(a,b) 的一组解 x′,y′x',y',并根据 gcd⁡(a,b)∣c\gcd(a,b)\mid c 判断方程有无解
  2. 令 x=cgcd⁡(a,b)x′,y=cgcd⁡(a,b)y′x=\frac{c}{\gcd(a,b)}x',y=\frac{c}{\gcd(a,b)}y',x,yx,y 为 ax+by=cax+by=c 的一组解
  3. 若需要求 xx 的最小非负整数解,则将 x=(x0 mod m+m) mod mx=(x_0 \bmod m + m) \bmod m 作为最终的答案,其中 m=bgcd⁡(a,b)m=\frac{b}{\gcd(a,b)}

线性同余方程

形式

ax≡c(modp) ax \equiv c \pmod{p}

形如这样的方程,就被称线性同余方程。

求解

任意的线性同余方程都可以与一个二元线性不定方程互相转化,比如上面的方程就与二元线性不定方程有如下关系:

ax≡c(modp)  ⟺  ax+py=c ax \equiv c \pmod{p} \iff ax+py=c

所以可以把线性同余方程转化为一个二元线性不定方程,求出 xx 即可,包括有无解的判定条件都是一样的。

代码

template <typename T> T exgcd(T a, T b, T & x, T & y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    } else {
        T d = exgcd(b, a % b, x, y);
        T t = x;
        x = y;
        y = t - a / b * y;
        return d;
    }
}

template <typename T> optional<T> solveEquation(T a, T c, T p) {
    a = (a % p + p) % p;
    c = (c % p + p) % p;
    T x, y;
    T d = exgcd(a, p, x, y);

    if (c % d != 0)
        return nullopt;
    
    x = (1ll * (c / d) * x % (p / d) + (p / d)) % (p / d);
    return x;
}

中国剩余定理

形式

{x≡a1(modm1)x≡a2(modm2)⋯x≡an(modmn)∀i≠j,gcd⁡(mi,mj)=1 \begin{cases} x \equiv a_1 \pmod{m_1}\\\\ x \equiv a_2 \pmod{m_2}\\\\ \cdots\\\\ x \equiv a_n \pmod{m_n}\\\\ \end{cases} \\\\ \forall i\ne j,\gcd(m_i,m_j)=1

形如上面这样的线性同余方程组且满足每个方程的模数两两互质,就可以用中国剩余定理解决。中国剩余定理也叫做 CRT (Chinese Remainder Theorem)。

求解

  1. 令 M=∏i=1nmiM=\prod_{i=1}^{n} m_i,MM 不对任何数取模
  2. 对于 mim_i,令 Mi=∏i≠jmj=MmiM_i=\prod_{i\ne j}m_j=\frac{M}{m_i},Mi−1M_i^{-1} 为 MiM_i 对 mim_i 的逆元,ci=MiMi−1c_i=M_iM_i^{-1},cic_i 不能对 mim_i 取模
  3. 令 x=∑i=1naici mod Mx=\sum_{i=1}^{n}a_ic_i \bmod M,则 xx 为模 MM 意义下的唯一解

算法正确性证明

对于 ∀i≠j\forall i\ne j,mi∣Mjm_i\mid M_j,故 mi∣ajcjm_i\mid a_jc_j,所以在模 mim_i 情况下,x≡aici(modmi)x \equiv a_ic_i \pmod{m_i}。

因为 Mi−1M_i^{-1} 为 MiM_i 对 mim_i 的逆元,故在模 mim_i 的同余式中 ci≡1(modmi)c_i\equiv 1 \pmod{m_i},aici≡ai(modmi)a_ic_i\equiv a_i \pmod{m_i},即 x≡ai(modmi)x \equiv a_i \pmod{m_i},满足第 ii 个方程。

解在模 MM 意义下的唯一性证明

设 x,yx,y 不相等,x,y<Mx,y<M,但都满足方程组,则可以有下面的方程组:

{x−y≡0(modm1)x−y≡0(modm2)⋯x−y≡0(modmn) \begin{cases} x - y\equiv 0 \pmod{m_1}\\\\ x - y\equiv 0 \pmod{m_2}\\\\ \cdots\\\\ x - y\equiv 0 \pmod{m_n}\\\\ \end{cases}

则 ∀i,mi∣∣x−y∣\forall i,m_i\mid |x-y|,所以 M∣∣x−y∣M \mid |x-y|,即 ∣x−y∣=kM|x-y|=kM,因为 ∣x−y∣<M|x-y|<M,故 k=0k=0,推出 x=yx=y,与假设矛盾。

代码

int exgcd(int a, int b, int & x, int & y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    } else {
        int d = exgcd(b, a % b, x, y);
        int t = x;
        x = y;
        y = t - (a / b) * y;
        return d;
    }
}

int solveEquation(int a, int c, int p) {
    int x, y;
    exgcd(a, p, x, y);
    x = (1ll * x * c % p + p) % p;
    return x;
}

inline int inv(int x, int p) {
    return solveEquation(x, 1, p);
}

int chineseRemainderTheorem(int n, int a[], int m[]) {
    int mod = 1, ans = 0;

    for (int i = 1; i <= n; ++i)
        mod *= m[i];
    
    for (int i = 1; i <= n; ++i) {
        int modi = mod / m[i];
        ans = (ans + 1ll * a[i] * modi * inv(modi)) % mod; // 用 exgcd 求逆元
    }

    return ans;
}

扩展中国剩余定理

形式

{a1x≡c1(modm1)a2x≡c2(modm2)⋯anx≡cn(modmn) \begin{cases} a_1 x \equiv c_1 \pmod{m_1}\\\\ a_2 x \equiv c_2 \pmod{m_2}\\\\ \cdots\\\\ a_n x \equiv c_n \pmod{m_n}\\\\ \end{cases}

相比普通中国剩余定理求解的方程组,少了模数互质的限制。

求解

考虑合并方程。设当前要合并前 ii 个方程,显然可以转化为合并前 i−1i-1 个方程与第 ii 个方程。

设 xix_i 为合并前 ii 个方程的解,di=lcm⁡(m1,m2,…,mi)d_i=\operatorname{lcm}(m_1,m_2,\dots,m_i)。构造 xi=sxi−1+tx_i=sx_{i-1}+t,我们要使 xix_i 能够满足前 i−1i-1 个方程,就要满足以下条件:

  1. ∀j<i\forall j < i,aj(sxi−1+t)≡ajxi−1(modmj)a_j(sx_{i-1}+t)\equiv a_j x_{i-1} \pmod{m_j}
  2. ai(sxi−1+t)≡ci(modmi)a_i(sx_{i-1}+t) \equiv c_i \pmod{m_i}

最简单的办法就是令 s=1,t=di−1yis=1,t=d_{i-1}y_i,并求一个合适的 yiy_i,也就是要求出以下方程的解 yiy_i:

aidi−1yi≡ci−aixi−1(modmi) a_i d_{i-1}y_i \equiv c_i - a_i x_{i-1} \pmod{m_i}

用 exgcd 解决,若该线性同余方程无解,则方程组无解,否则求出前 ii 个方程的解 xix_i,进行下一次合并。

假设各个 mim_i 同阶,时间复杂度 O(nlog⁡2mi)O(n \log_2 m_i)。

代码

int gcd(int a, int b) {
    return b == 0 ? a : gcd(b, a % b);
}

int exgcd(int a, int b, int & x, int & y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    } else {
        int d = exgcd(b, a % b, x, y);
        int t = x;
        x = y;
        y = t - (a / b) * y;
        return d;
    }
}

optional<int> solveEquation(int a, int c, int p) {
    a = (a % p + p) % p;
    c = (c % p + p) % p;
    int x, y, d = exgcd(a, p, x, y);
    
    if (c % d)
        return nullopt;
    
    x *= (c / d);
    int m = p / d;
    x = (x % m + m) % m;
    return x;
}

optional<int> exChineseRemainderTheorem(int n, int a[], int c[], int m[]) {
    int x = 0, d = 1;
    
    for (int i = 1; i <= n; ++i) {
        auto oy = solveEquation(a[i] * d, c[i] - a[i] * x, m[i]);

        if (!oy.has_value())
            return nullopt;
        
        int y = oy.value();
        int g = gcd(d, m[i]);
        d = d / g * m[i];
        x = (x + 1ll * (d / m[i] * g) * y) % d;
    }

    return x;
}

BSGS

离散对数

与离散对数相对的是我们熟悉的连续对数,也就是最经常接触的对数。而离散对数是在同余中的一个概念。

若在模 mm 意义下,xx 满足以下等式:

ax≡b(modm) a^x \equiv b \pmod{m}

则 xx 为以 aa 为底的 bb 的离散对数。注意离散对数不同于连续对数,一个数的离散对数有很多个,并不是唯一的。

如果 gcd⁡(a,m)=1\gcd(a,m)=1,则我们可以使用 BSGS 算法求解一个数的最小的离散对数。

求解

因为 gcd⁡(a,m)=1\gcd(a,m)=1,则 aka^k 对 mm 的逆元一定存在。所以我们可以把 axa^x 写成 aps−qa^{ps-q} 并直接将 a−qa^{-q} 提取出来,即:

aps≡baq(modm) a^{ps} \equiv ba^q \pmod{m}

考虑暴力的思路,可以枚举两边的 p,qp,q,如果有同余的两个 aps,baqa^{ps},ba^q,就更新答案。

考虑到枚举会计算许多重复的东西,可以用某种数据结构维护映射表,存下每个 (baq,q)(ba^q,q),再枚举左边找相应的同余的项。我们先假设 x=ps−q∈[0,m)x=ps-q \in [0,m),则当 s=⌈m⌉s=\lceil\sqrt{m}\rceil 时,时间复杂度最优,为 Θ(m)\Theta(\sqrt{m}) 或 Θ(mlog⁡2m)\Theta(\sqrt{m}\log_2 \sqrt{m}),具体为哪种复杂度,取决于是用哈希表还是平衡树实现。

可解性证明

至于为什么一定有 x∈[0,m)x\in [0,m),使用抽屉原理证明:

对于任意的 xx,ax mod ma^x \bmod m 最多只有 mm 种取值,而在 x∈[0,m]x\in [0,m] 时,一共有 m+1m+1 种指数,则根据抽屉原理得出必定存在一对 x,yx,y 满足 ax≡ay(modm)a^x \equiv a^y \pmod{m}。设 t=∣x−y∣t=|x-y|,则 axak=ax+taka^x a^k=a^{x+t}a^k,因为已经保证 gcd⁡(a,m)=1\gcd(a,m)=1,则逆元一定存在,所以可以规定 t∈Zt\in\mathrm{Z},所以指数具有周期性,一定可以得到 x∈[0,m)x\in [0,m)。

代码

int power(int x, int y, int m) {
    int res = 1;

    for (; y; y /= 2) {
        if (y % 2)
            res = 1ll * res * x % m;
        
        x = 1ll * x * x % m;
    }

    return res;
}

optional<int> bsgs(int a, int b, int m) {
    a %= m;
    b %= m;
    
    unordered_map<int, int> buc;
    int s = ceil(sqrt(m)), prod = b, base = power(a, s, m);
    buc[b] = 0;

    for (int i = 0; i < s; ++i) {
        prod = 1ll * prod * a % m;
        buc[prod] = i;
    }

    prod = 1;

    for (int i = 0; i <= s; ++i) {
        auto it = buc.find(prod);

        if (it != buc.end() && i * s - it->second >= 0)
            return i * s - it->second;
            
        prod = 1ll * prod * base % m;
    }

    return nullopt;
}

扩展 BSGS

形式

ax≡b(modm) a^x \equiv b \pmod{m}

其中 a,ma,m 不一定互质。这种情况下,逆元不一定存在,也就不可以把指数拆成两个部分来优化。

考虑先把式子化为 gcd⁡(a,m)=1\gcd(a,m)=1 的形式。方法是不断提取出 aa 与 mm 互质的部分,任何整个式子同除以这个数。

a≡c(modm)  ⟺  ad≡cd(modmd) (d∣a∧d∣c∧d∣m) a \equiv c \pmod{m} \iff \frac{a}{d} \equiv \frac{c}{d} \pmod{\frac{m}{d}}\ (d\mid a \wedge d\mid c \wedge d\mid m)

这个式子是同余式的基本性质。所以按照这样的方法提取:

ax≡b(modm)⇓ad1ax−1≡bd1(modmd1)⇓a2d1d2ax−2≡bd1d2(modmd1d2)⇓⋮⇓akd1d2⋯dkax−k≡bd1d2⋯dk(modmd1d2⋯dk) a^x \equiv b \pmod{m}\\\\ \Downarrow\\\\ \frac{a}{d_1}a^{x-1} \equiv \frac{b}{d_1} \pmod{\frac{m}{d_1}}\\\\ \Downarrow\\\\ \frac{a^2}{d_1 d_2}a^{x-2} \equiv \frac{b}{d_1 d_2} \pmod{\frac{m}{d_1 d_2}}\\\\ \Downarrow\\\\ \vdots\\\\ \Downarrow\\\\ \frac{a^k}{d_1 d_2\cdots d_k}a^{x-k} \equiv \frac{b}{d_1 d_2\cdots d_k} \pmod{\frac{m}{d_1 d_2\cdots d_k}}\\\\

假设现在 aa 与 md1d2⋯dk\frac{m}{d_1 d_2\cdots d_k} 互质,就可以用 BSGS 求出最后一个方程的解 x−kx-k,再加上 kk 就是 xx 了。

需要注意两点:

  1. 若在提取过程中,出现 d1d2⋯di∤bd_1 d_2 \cdots d_i \nmid b,则方程无法再继续提取,方程无解
  2. 有可能答案在 [0,k)[0,k) 内,只要 O(k)O(k) 枚举检验即可。
int gcd(int a, int b) {
    int t;

    while (b) {
        t = a;
        a = b;
        b = a % b;
    }
    
    return a;
}

int power(int x, int y, int m) {
    int res = 1;

    for (; y; y /= 2) {
        if (y % 2)
            res = 1ll * res * x % m;
        
        x = 1ll * x * x % m;
    }

    return res;
}

int exgcd(int a, int b, int & x, int & y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    } else {
        int d = exgcd(b, a % b, x, y);
        int t = x;
        x = y;
        y = t - a / b * y;
        return d;
    }
}

int inverse(int a, int p) {
    int x, y;
    exgcd(a, p, x, y);
    x = (x % p + p) % p;
    return x;
}

optional<int> bsgs(int a, int b, int m) {
    a %= m;
    b %= m;
    
    unordered_map<int, int> buc;
    int s = ceil(sqrt(m)), prod = b, base = power(a, s, m);
    buc[b] = 0;

    for (int i = 0; i < s; ++i) {
        prod = 1ll * prod * a % m;
        buc[prod] = i;
    }

    prod = 1;

    for (int i = 0; i <= s; ++i) {
        auto it = buc.find(prod);

        if (it != buc.end() && i * s - it->second >= 0)
            return i * s - it->second;
            
        prod = 1ll * prod * base % m;
    }

    return nullopt;
}

optional<int> exBsgs(int a, int b, int m) {
    a = (a % m + m) % m;
    b = (b % m + m) % m;
    
    if (m == 1 || b == 1)
        return 0;
        
    int prod = 1, d, k = 0;

    while (true) {
        d = gcd(a, m);
        
        if (d == 1)
            break;
        
        if (b % d)
            return nullopt;
        
        b /= d;
        m /= d;
        ++k;
        prod = prod * (a / d) % m;

        if (prod == b)
            return k;
    }

    auto res = bsgs(a, 1ll * b * inverse(prod, m) % p, m);

    if (!res.has_value())
        return nullopt;
    
    return res.value() + k;
}

阶与原根

阶

对于一对互质的整数 a,ma,m,若 rr 满足 r≠0r \ne 0 且 ar≡1(modm)a^r \equiv 1 \pmod{m} 且不存在任何小于 rr 的正整数满足该同余式,则称 rr 为 aa 对 mm 的阶,记作 r=δm(a)r = \delta_m(a)。

Theorem 1:0<δm(a)≤φ(m)0 < \delta_m(a) \le \varphi(m)。

若 gcd⁡(a,m)=1\gcd(a,m)=1,则欧拉定理一定成立,即 aφ(m)≡1(modm)a^{\varphi(m)} \equiv 1 \pmod{m},显然 0<δm(a)≤φ(m)0 < \delta_m(a) \le \varphi(m)。

Theorem 2:若 ak≡1(modm)a^k \equiv 1 \pmod{m},则 δm(a)∣k\delta_m(a) \mid k。

假设满足 δm(a)∤k\delta_m(a) \nmid k 但 ak≡1(modm)a^k \equiv 1 \pmod{m},则 k=sδm(a)+tk=s\delta_m(a)+t,其中 t∈(0,δm(a))t \in (0,\delta_m(a)),所以:

ak≡(aδm(a))sat≡at(modm) a^k \equiv (a^{\delta_m(a)})^sa^t \equiv a^t \pmod{m}

根据阶的定义,at≢1(modm)a^t \not\equiv 1 \pmod{m},所以矛盾,故 δm(a)∣k\delta_m(a) \mid k。

原根

若 aa 满足 δm(a)=φ(m)\delta_m(a)=\varphi(m),则 aa 为 mm 的原根。

Theorem 1:aa 为 mm 的原根的充要条件为 gcd⁡(a,m)=1\gcd(a,m)=1 且对于 φ(m)\varphi(m) 的每一个质因数,aφ(m)p≢1(modm)a^{\frac{\varphi(m)}{p}} \not\equiv 1 \pmod{m}。

根据这个定理,可以轻松判断一个数是否是另一个数的原根。

int gcd(int a, int b) {
    int t;

    while (b) {
        t = a;
        a = b;
        b = t % b;
    }

    return a;
}

int phi(int x) {
    int res = x;

    for (int i = 2; i * i <= x; ++i) {
        if (x % i)
            continue;
        
        res = res / i * (i - 1);

        while (x % i)
            x /= i;
    }

    if (x > 1)
        res = res / x * (x - 1);
    
    return res;
}

int power(int x, int y, int m) {
    int res = 1;

    for (; y; y /= 2) {
        if (y % 2)
            res = 1ll * res * x % m;
        
        x = 1ll * x * x % m;
    }

    return res;
}

bool isPrimitiveRoot(int a, int m) {
    if (gcd(a, m) != 1)
        return false;
    
    int pm = phi(m), tmp = pm;
    a %= m;

    for (int i = 2; i * i <= tmp; ++i) {
        if (tmp % i)
            continue;
        
        if (power(a, pm / i, m) == 1)
            return false;
        
        while (tmp % i == 0)
            tmp /= i;
    }

    if (tmp > 1)
        if (power(a, pm / tmp, m) == 1)
            return false;

    return true;
}

Theorem 2:一个数 mm 有原根,当且仅当 mm 为 2,4,pk,2pk2,4,p^k,2p^k,其中 k∈N∗k\in \mathrm{N^*},pp 为奇质数。

Theorem 3:若一个数 mm 有原根,则其原根个数为 φ(φ(m))\varphi(\varphi(m))。

Theorem 4:若一个数 mm 有原根,则其最小原根的数量级为 m14m^{\frac{1}{4}}。

根据这个性质,只需要暴力枚举每个数,求一个数的最小原根的时间复杂度是可以接受的。

int gcd(int a, int b) {
    int t;

    while (b) {
        t = a;
        a = b;
        b = t % b;
    }

    return a;
}

int phi(int x) {
    int res = x;

    for (int i = 2; i * i <= x; ++i) {
        if (x % i)
            continue;
        
        res = res / i * (i - 1);

        while (x % i)
            x /= i;
    }

    if (x > 1)
        res = res / x * (x - 1);
    
    return res;
}

int power(int x, int y, int m) {
    int res = 1;

    for (; y; y /= 2) {
        if (y % 2)
            res = 1ll * res * x % m;
        
        x = 1ll * x * x % m;
    }

    return res;
}

int primitiveRoot(int m) {
    int pm = phi(m), tmp = pm;
    vector<int> factor;

    for (int i = 2; i * i <= tmp; ++i) {
        if (tmp % i)
            continue;
        
        factor.push_back(i);
        
        while (tmp % i == 0)
            tmp /= i;
    }

    if (tmp > 1)
        factor.push_back(tmp);
    
    for (int i = 2; i <= m; ++i) {
        if (gcd(i, m) != 1)
            continue;

        bool found = true;
        
        for (int p : factor) {
            if (power(i, pm / p, m) == 1) {
                found = false;
                break;
            }
        }

        if (found)
            return i;
    }

    return -1;
}

Theorem 5:若 aa 为 mm 的原根,则 a,a2,a3,…,aφ(m)a,a^2,a^3,\dots,a^{\varphi(m)} 两两模 mm 不同余。

Theorem 6:若 aa 为 mm 的原根,则对于任意的 x∈[1,φ(m)]x \in [1,\varphi(m)],一定有一个 kk 与其一一对应,满足 x≡ak(modm)x \equiv a^k \pmod{m} 且 k∈[1,φ(m)]k\in [1,\varphi(m)],kk 被称为 xx 的指标,记作 I(x)I(x)。

比如 m=13m=13,其最小原根 a=2a=2,则把 xx 与 I(x)I(x) 写成表格:

xx112233445566778899101011111212
I(x)I(x)121211442299551111338810107766

这个定理可以用于把底数转换为原根的幂,可以利用这个性质解模数为质数的 NN 次剩余问题。

N 次剩余

形式

xn≡a(modm) x^n \equiv a \pmod{m}

其中 mm 为奇质数,求 x∈[0,m)x\in [0,m)。由于 mm 不是奇质数的方法较复杂,这里不考虑。

求解

设 gg 为其原根,则可以利用原根和指标的性质将 xx 转化为 gyg^y,即:

(gy)n≡(gn)y≡a(modm) (g^y)^n \equiv (g^n)^y \equiv a \pmod{m}

这样就可以用 BSGS 求出 yy,再用快速幂求得 xx 即可。

也可以再将 aa 化为 gzg^z,然后方程化为:

gny≡gz(modm) g^{ny} \equiv g^z \pmod{m}

于是就可以解出线性同余方程的解 yy:

ny≡z(modm) ny \equiv z \pmod{m}

代码就不写了,上面组合一下就可以了。