【数学】同余方程解法合集

更新于 分类于 ， 阅读次数： 本文字数： 6.2k 阅读时长 ≈ 6 分钟

包括（一元、多元）线性同余方程（组）、高次同余方程的解法模板。

线性同余方程

一元线性同余方程

形式：ax≡b(mod  m)ax\equiv b(\mod m) ，求解xx ，引入逆元后可转化为x≡ba−1(mod  m)x\equiv ba^{-1}(\mod m)

等价式ax+my=b,y∈Nax+my=b,y\in N

有解：gcd⁡(a,m)=d\gcd(a,m)=d ，d∣bd\mid b 则有dd 个模mm 不同余的解，否则无解

解法：转化为ax+my=dax+my=d ，则乘b/db/d 即可得到原式的解。

1
2
3
4
5
6
7
8
9
10
11
// 求解ax+my=b的x，即ax==b(mod m)，答案存在ans[0...d-1]中，返回值d为答案个数
int mod_equation(int a, int b, int m)
{
    exgcd(a, m, d, x, y);
    if (b % d) // 无解
        return -1;
    x = x * (b / d) % m;
    for (int i = 0; i < d; i++)
        ans[i] = (x + i * m / d) % m;
    return d;
}

一元线性同余方程组

形如：

{x≡a1(mod  m1)x≡a2(mod  m2)⋮x≡an(mod  mn)\begin{cases} x\equiv a_1&(\mod m_1)\\ x\equiv a_2&(\mod m_2)\\ &\vdots \\ x\equiv a_n&(\mod m_n)\\ \end{cases}

方法一：两两相消

先看两个方程：

{x≡a1(mod  m1)x≡a2(mod  m2)\begin{cases} x\equiv a_1(\mod m_1)\\ x\equiv a_2(\mod m_2)\\ \end{cases}

转化为

\begin{numcases}{} x+m_1k_1=a_1\tag{1}\\ x+m_2k_2=a_2\\ \end{numcases}

则得到a1−m1k1=a2−m2k2a_1-m_1k_1=a_2-m_2k_2 ，即m1k1−m2k2=a1−a2m_1k_1-m_2k_2=a_1-a_2 ，则有解的等价条件为gcd⁡(m1,m2)∣(a1−a2)\gcd(m_1,m_2)\mid(a_1-a_2) 。用上一节的一元线性同余方程解法解出k1k_1 的一个特解k′k' ，于是通解k1=k′+cm2gcd⁡(m1,m2),c∈Nk_1=k'+\frac{cm_2}{\gcd(m_1,m_2)},c\in N 。代入(1)(1) 式可得，x=a1−m1k′+cm1m2gcd⁡(m1,m2)x=a_1-m_1k'+\frac{cm_1m_2}{\gcd(m_1,m_2)} 。于是转换为了新的同余方程：x≡a1−m1k′(mod  lcm(m1,m2))x\equiv a_1-m_1k'(\mod \text{lcm}(m_1,m_2)) 。依次合并即可解得最终答案。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
// 求解 x==a[] mod m[]，数组从0开始编号，共n个式子
int mod_equations(int a[], int m[], int n)
{
    int a1 = a[0], a2, m1 = m[0], m2, k1, k, d;
    for (int i = 1; i < n; i++)
    {
        a2 = a[i], m2 = m[i];
        int k2, tmp;
        exgcd(m1, m2, d, k1, k2);
        if ((a2 - a1) % d)
            return -1;
        tmp = m2 / d;
        k1 = (k1 * (a2 - a1) / d % tmp + tmp) % tmp; // 特解k'，保证最小正数
        a1 += m1 * k1;
        m1 *= m2 / d;
        a1 = (a1 + m1) % m1;
    }
    return a1; // 通解是a1+cm1,c是整数
}

方法二：中国剩余定理

条件：原方程组m1,m2⋯mnm_1,m_2\cdots m_n 互素

定理：有模M=m1⋅m2⋯mnM=m_1\cdot m_2\cdots m_n 的唯一解

推导：令Mi=M/miM_i=M/m_i ，由互素知gcd⁡(Mi,mi)=1\gcd(M_i,m_i)=1 ，又Mi⋅Mi−1≡1(mod  mi)M_i\cdot M_i^{-1}\equiv 1(\mod m_i) ，则通解为cM+∑i=1naiMiMi−1cM+\sum_{i=1}^{n}a_iM_iM_i^{-1} 。

1
2
3
4
5
6
7
8
9
10
11
12
13
// 中国剩余定理解一元线性同余方程组，数组从0标号
ll CRT(int a[], int m[], int n)
{
    ll M = 1ll, ans = 0, tmp;
    for (int i = 0; i < n; i++)
        M *= m[i];
    for (int i = 0; i < n; i++)
    {
        tmp = M / m[i];
        ans = (ans + a[i] * tmp % M * inv(tmp) % M) % M;
    }
    return ans;
}

多元线性同余方程

形如：a1x1+a2x2+⋯+anxn≡b(mod  m)a_1x_1+a_2x_2+\cdots+a_nx_n\equiv b(\mod m)

由a1x+a2y=gcd⁡(a1,a2)a_1x+a_2y=\gcd(a_1,a_2) 有解可知左式能表达cgcd⁡(a1,a2),c∈Nc\gcd(a_1,a_2),c\in N ；引入a3a_3 ，则gcd(a1,a2)x+a3ygcd(a_1,a_2)x+a_3y 可以表示cgcd⁡(a1,a2,a3),c∈Nc\gcd(a_1,a_2,a_3),c\in N 。以此类推，a1x1+a2x2+⋯+anxna_1x_1+a_2x_2+\cdots+a_nx_n

可以表达cgcd⁡(a1,a2,⋯ ,an),c∈Nc\gcd(a_1,a_2,\cdots,a_n),c\in N ，因此有解的条件为gcd⁡(a1,a2,⋯ ,an,m)∣b\gcd(a_1,a_2,\cdots, a_n,m)\mid b 。

一个很好理解的推导：我们将mm 添加到aa 数列中，仍设数量为nn 。设d=gcd⁡(a1,a2,⋯ ,an),d1=gcd⁡(a1,a2,⋯ ,an−1)d=\gcd(a_1,a_2,\cdots, a_n),d_1=\gcd(a_1,a_2,\cdots, a_{n-1}) ，于是有gcd⁡(d1,an)=d\gcd(d_1,a_n)=d ，则原式转化为anxn+d1k=ba_nx_n+d_1k=b ，可以通过前面一元线性同余方程的解法解出xn,kx_n,k ，则问题转化为a1x1+a2x2+⋯+an−1xn−1=d1∗ka_1x_1+a_2x_2+\cdots+a_{n-1}x_{n-1}=d_1*k，仍采取相同的解法。

但O(n2)O(n^2) 求每次的gcd⁡\gcd 显然不现实。设di=gcd⁡(a1,a2,⋯ ,ai)d_i=\gcd(a_1,a_2,\cdots,a_{i}) ，x[1...n]x[1...n] 表示解（从11 标号，a[n+1]a[n+1] 是mm ）则可以求解diXi+ai+1yi=di+1d_iX_i+a_{i+1}y_i=d_{i+1} 的Xi,yiX_i,y_i 。观察最后一个式子：dnXn+an+1yn=dn+1d_{n}X_{n}+a_{n+1}y_{n}=d_{n+1} ，要求解答案，等式两边乘上c=bdn+1c=\frac{b}{d_{n+1}} ，所以xn+1=yncx_{n+1}=y_{n}c ，b−an+1xn+1=dnXncb-a_{n+1}x_{n+1}=d_{n}X_{n}c ；则前一个式子dn−1Xn−1+anyn−1=dnd_{n-1}X_{n-1}+a_{n}y_{n-1}=d_{n} ，要求解答案，等式两边乘上c′=cXnc'=cX_{n} ，依次求解即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 求解a1x1+a2x2+...+anxn==b mod m，a从1开始标号
int X[N], y[N];
bool mod_Xn_equation(int a[], int b, int m, int n, int x[])
{
    a[n + 1] = m;
    int d = a[1];
    for (int i = 2; i <= n + 1; i++)
        d = exgcd(d, a[i], X[i - 1], y[i - 1]); // 这个exgcd返回求解的ax+by=d的
    if (b % d)
        return false;
    int c = b / d;
    x[n + 1] = y[n] * c;
    y[0] = 1;
    for (int i = n; i >= 1; i--)
    {
        c = (c * X[i] % m + m) % m;
        x[i] = (y[i - 1] * c % m + m) % m;
    }
    return true;
}

高次同余方程

只考虑：ax≡b(mod  m)a^x\equiv b(\mod m)

为保持统一，放一下手写哈希

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// (x,y)的HASH
#define MOD 76543
int hs[MOD], head[MOD], nxt[MOD], id[MOD], top;
void insert(int x, int y)
{
    int k = x % MOD;
    hs[top] = x, id[top] = y, nxt[top] = head[k], head[k] = top++;
}
int find(int x)
{
    int k = x % MOD;
    for (int i = head[k]; i != -1; i = nxt[i])
        if (hs[i] == x)
            return id[i];
    return -1;
}

BSGS算法

条件：gcd⁡(a,m)=1\gcd(a,m)=1

令x=ks−t,s=⌈m⌉,0≤t<s,0<k≤s+1x=ks-t,s=\lceil \sqrt m\rceil,0\leq t<s,0<k\leq s+1 （ss 的取值有严格证明这样复杂度最低），则原式变为aks≡bat(mod  m)a^{ks}\equiv ba^{t}(\mod m) ，枚举可能的tt ，存下batba^t 的哈希值，再枚举左边的kk 算aksa^{ks} 查表。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 普通BSGS: a^x==b mod m, gcd(a,m)=1
int BSGS(int a, int b, int m)
{
    if(b == 1 || m == 1) // 注意这里的特判
        return 0;
    int s = sqrt(m) + 1;
    memset(head,-1,sizeof(head));
    top=0;
    for (int i = 0, tmp = b; i < s; i++, tmp = tmp * a % m)
        insert(tmp, i);
    int as = qpow(a, s, m); // 快速幂算a^s
    for (int i = 1, tmp = as; i <= s + 1; i++, tmp = tmp * as % m)
    {
        int qry = find(tmp);
        if (qry == -1)
            continue;
        return i * s - qry;
    }
    return -1;
}

扩展BSGS算法

条件：gcd⁡(a,m)=d≠1\gcd(a,m)=d\neq 1

等价变换为ax+my=ba^x+my=b ，有解的条件是d∣bd\mid b 。两边同除dd ，则变为adax−1+ymd=bd\frac{a}{d}a^{x-1}+y\frac{m}{d}=\frac{b}{d} ，不断检查gcd⁡(md,a)\gcd(\frac{m}{d},a) 并除掉，直到为11 ，记每次除的数乘积为DD ，则原式变为akDax−k+ymD=bD\frac{a^k}{D}a^{x-k}+y\frac{m}{D}=\frac{b}{D} ，即ax−k≡bDinv(akD)(mod  mD)a^{x-k}\equiv \frac{b}{D}inv(\frac{a^k}{D})(\mod \frac{m}{D}) ，用BSGS求解即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
int exBSGS(int a, int b, int m)
{
    if (b == 1 || m == 1)
        return 0;
    int k = 0, atmp = 1;
    while (1)
    {
        int tmp = __gcd(m, a);
        if (tmp == 1)
            break;
        if(b % tmp)
            return -1;
        b /= tmp, m /= tmp, k++;
        atmp = atmp * (a / tmp) % m;
        if (b == atmp)
            return k;
    } // 此时b=b/D, m=m/D, atmp=a^k / D
    int ans = BSGS(a, b * inv(atmp, m) % m, m);
    return ans == -1 ? ans  + k;
}