[笔记]同余方程
本文最后更新于 245 天前,其中的信息可能已经有所发展或是发生改变。

拓展欧几里得算法、中国剩余定理、拓展中国剩余定理、BSGS、exBSGS

Before the Beginning

转载请将本段放在文章开头显眼处,如有二次创作请标明。

原文链接:https://www.codein.icu/congruence/

拓展欧几里得算法

解不定方程 $ax + by = c$,可以使用拓展欧几里得算法。

首先解 $ax + by = \gcd (a,b)$.

欧几里得算法

证明 $\gcd(a,b) = \gcd(b,a \bmod b)$:

设 $a = g \times k_1$,$b = g \times k_2$,其中 $k_1,k_2$ 互质。

要证明 $\gcd(a,b) = \gcd(b,a\bmod b)$,即证 $g = \gcd(g \times k_2, (g \times k_1 )\bmod (g \times k_2))$.

将 $g$ 消去,即证 $\gcd(k_2,k_1 \bmod k_2) = 1$.

假设有公因子 $x$,则:

$k_2 = x \times n_1$,$k_1 \bmod k_2 = x \times n_2$.

可以写出:$k_1 = k \times k_2 + k_1 \bmod k_2$,其中 $k$ 为整数。

化简得出:$k_1 = k \times (x \times n_1) + x \times n_2$,则 $k_1,k_2$ 有公因子 $x$,而 $k_1,k_2$ 应当互质,矛盾。

所以 $\gcd(k_2,k_1 \bmod k_2) = 1$.

则 $\gcd(a,b) = \gcd(b,a\bmod b)$.

拓展欧几里得算法

证明 $ax + by = \gcd(a,b)$ 一定有解。$(a,b \neq 0)$.

由上文结论,该式等价于:

$bx + (a \bmod b)y = \gcd(b,a \bmod b)$.

辗转相除,当 $b = 0$ 时,方程为:$ax = \gcd(a,b) = a$,即用欧几里得算法求最大公因数的最后一步。

有 $x = 1,y = 0$.

向上还原,即可获得一组解。


拓展欧几里得算法用于求出 $ax + by = \gcd(a,b)$ 的一组特解:

$ax_1 + by_1 = \gcd(a,b)$ 可转化为:

$bx_2 + (a \bmod b)y_2 = \gcd(b,a \bmod b)$.

其中:$\gcd(a,b) = \gcd(b,a \bmod b)$.

可以得到:$ax_1 + by_1 = bx_2 + (a \bmod b)y_2$.

即:$ax_1 + by_1 = bx_2 + (a – \lfloor \dfrac{a}{b} \rfloor \times b)y_2$.

展开、合并同类项得到:$ax_1 + by_1 = ay_2 + b(x_2 – y_2 \times \lfloor \dfrac{a}{b} \rfloor)$

即:

$x_1 = y_2$,$y_1 = x_2 – y_2 \times \lfloor \dfrac{a}{b} \rfloor = x_2 – x_1 \times \lfloor \dfrac{a}{b} \rfloor$.

可以据此写出拓展欧几里得算法。可以使用引用传值的技巧简化代码。

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

此时回归到解 $ax + by = c$ 的问题上来。若 $\gcd(a,b) \mid c$,则原方程有整数解。

$ax’ + by’ = \gcd(a,b)$,则:

$$x = x’ \times \dfrac{c}{\gcd(a,b)}$$

$$y = y’ \times \dfrac{c}{\gcd(a,b)}$$

若要得到任意组解呢?

记求出的特解为 $x_0,y_0$.

$ax + by = ax_0 + by_0$.

$a(x – x_0) = b(y_0-y)$.

$\dfrac{a}{b} \times (x – x_0) = y_0 – y$.

$x – x_0$ 需要为整数,则 $\dfrac{b}{\gcd(a,b)} \mid (x’ – x_0′)$,设为 $x’ – x_0′ = \dfrac{b}{\gcd(a,b)} \times t$.

可以写出:

$x’ = x_0′ + \dfrac{b}{\gcd(a,b)} \times t$

$y’ = y_0′ – \dfrac{a}{\gcd(a,b)} \times t$

由此即可求出任意解,再带回原式即可。

例如求 $x’$ 的最小正整数解,可以发现 $x’$ 任意加减 $b$ 都是合法解,因此可以直接 $((x \bmod b) + b) \bmod b$.


拓展欧几里得算法还可以用于求解线性同余方程。

形如 $ax \equiv c \pmod{b}$.

可以写成 $ax + by = c$,随后按二元一次不定方程的方法求解即可。


例题:青蛙的约会

两只青蛙初始坐标为 $x,y$,每次跳 $m,n$ 米,定义相遇为模 $L$ 意义下坐标相等,求一起跳跃多少次后相遇。

形式地,$x + m \times t = y + n \times t \pmod L$.

$$x + m \times t + L \times k = y + n \times t$$

$$(m-n) \times t + L \times k = y – x$$

求 $t$ 的最小正整数解。

套上板子即可。

#include <algorithm>
#include <cstdio>
#define ll long long
ll sx, sy, m, n, L;
inline ll gcd(ll a, ll b) { return b == 0 ? a : gcd(b, a % b); }
inline void exgcd(ll a, ll b, ll& x, ll& y)
{
    if (b == 0) return (void)(x = 1, y = 0);
    exgcd(b, a % b, y, x), y -= (a / b) * x;
}
int main()
{
    scanf("%lld %lld %lld %lld %lld", &sx, &sy, &m, &n, &L);
    if(m - n < 0) std::swap(m,n),std::swap(sx,sy);
    ll g = gcd(m - n, L);
    ll a = m - n, b = L, x, y, c = sy - sx;
    if (c % g) { puts("Impossible"); return 0; }
    exgcd(a, b, x, y);
    x *= c / g, y *= c / g;
    ll mod = b / g;
    x = ((x % mod) + mod) % mod;
    printf("%lld\n", x);
    return 0;
}

中国剩余定理

$$\begin{cases} x \equiv a_1 \pmod {p_1} \\ x \equiv a_2 \pmod {p_2} \\ \ldots \\ x \equiv a_i \pmod {p_i} \end{cases}$$

中国剩余定理用于求解模数两两互质的线性同余方程组。

对于这种 $x \equiv a_i \pmod{p_i}$ 的线性同余方程组,在 $\bmod \prod \limits _{i=1}^n p_i$ 的意义下有唯一解。


设 $M = \prod \limits _{i=1}^n p_i$,而 $M_i = \dfrac{M}{p_i}$.

则解 $x = \sum \limits _{i=1}^{n} M_i \times K_i$,其中 $K_i$ 满足:

$$M_i \times K_i \equiv a_i \pmod{p_i}$$

对于任意的 $M_i \times K_i$,由于 $M_i$ 是所有剩余模数的积,则 $M_i \times K_i \equiv 0 \pmod{p_j} (j \neq i)$,因此在 $x$ 中加上 $M_i \times K_i$ 对满足 $x \equiv a_j \pmod{p_j} (j \neq i)$ 无影响。

而每项 $M_i \times K_i$ 都能保证第 $i$ 个线性同余方程得到满足。

只需要解出 $K_i$ 即可。

$$M_i \times K_i + p_i \times t = a_i$$

利用扩展欧几里得定理求出 $K_i$.

由于模数两两互素,$\gcd(M_i,p_i) = 1$,即解:

$M_i \times k_i + p_i \times t = 1$ 的 $k_i$,然后 $K_i = a_i \times k_i$.

可以发现,$k_i$ 就是 $M_i$ 在 $\bmod p_i$ 意义下的逆元。

那么就可以得到一个解:$x = \sum \limits _{i=1}^n M_i \times k_i \times a_i$


而通解为:$x = k \times M + \sum \limits _{i=1}^n M_i \times k_i \times a_i$

证明:设两个解为 $x_1,x_2$,有:

$$\forall i \in \{1,2,\ldots,n\},x_1 – x_2 \equiv 0 \pmod{p_i}$$

而 $p_i$ 两两互素,则 $M \mid x_1 – x_2$,那么 $x_1 – x_2 = k \times M$,则在 $\bmod M$ 意义下有唯一解。


inline ll CRT()
{
    ll prod = 1,ans = 0;
    for (int i = 1; i <= n; ++i) prod *= a[i];
    for (int i = 1; i <= n; ++i)
    {
        ll x,y;
        exgcd(prod / a[i], a[i], x, y);
        ans = (ans + prod / a[i] * b[i] * x) % prod;
    }
    return (ans + prod) % prod;
}

例题:[TJOI2009]猜数字

给出 $a_i,b_i$,要求最小的 $n$ 使得 $b_i \mid (n – a_i)$.

改写为 $n – a_i \equiv 0 \pmod {b_i}$.

即 $n \equiv a_i \pmod{b_i}$

求 $n$ 的最小整数解。

这题会爆 long long,需要用到龟速乘的技巧。

#include <cstdio>
#define ll long long
const int maxk = 20;
ll k, a[maxk], b[maxk];
void exgcd(ll a, ll b, ll& x, ll& y)
{
    if (b == 0) return (void)(x = 1, y = 0);
    exgcd(b, a % b, y, x), y -= (a / b) * x;
}
ll mul(ll a,ll b,ll mod)
{
    ll res = 0;
    while(a)
    {
        if(a & 1) res = (res + b) % mod;
        b = (b + b) % mod;
        a >>= 1;
    }
    return res % mod;
}
ll crt()
{
    ll prod = 1, ans = 0;
    for (int i = 1; i <= k; ++i) prod *= b[i];
    for (int i = 1; i <= k; ++i)
    {
        ll x, y, mi = prod / b[i];
        exgcd(mi, b[i], x, y);
        x = ((x % b[i]) + b[i]) % b[i];
        ans = (ans + mul(mul(mi, x, prod), a[i], prod)) % prod;
    }
    return (ans + prod) % prod;
}
int main()
{
    scanf("%lld",&k);
    for (int i = 1; i <= k; ++i) scanf("%lld", a + i);
    for (int i = 1; i <= k; ++i) scanf("%lld", b + i);
    printf("%lld\n", crt());
    return 0;
}

拓展中国剩余定理

然而如果不保证模数互质,就需要使用拓展中国剩余定理。

模数不互质的话,之前的巧妙构造行不通了。尝试更平凡的列式子考虑:

$$\begin{cases}x \equiv r_1 \pmod{p_1} \\ x \equiv r_2 \pmod{p_2}\end{cases}$$

可以写成:

$$\begin{cases} x = r_1 + p_1 \times k_1 \\ x = r_2 + p_2 \times k_2 \end{cases}$$

即:

$r_1 + p_1 \times k_1 = r_2 + p_2 \times k_2$

移项得到:

$p_1 \times k_1 – p_2 \times k_2 = r_2 – r_1$

等价于解 $a = p_1$,$b = -p_2$,$c = r_2 – r_1$ 的 $ax + by = c$ 二元一次不定方程。

可以用拓展欧几里得定理解出可行的 $k_1,k_2$.

解出一组特解 $k_1′,k_2’$ 后,利用前文知识,求出通解:

记 $d=\gcd(p_1,-p_2)$

$$\begin{cases} k_1 = k_1′ + t \times \dfrac{-p_2}{d}\\ k_2 = k_2′ – t \times \dfrac{p_1}{d}\end{cases}$$

$x = r_1 + p_1 \times k_1$ 可以写为:

$x = r_1 + p_1 \times k_1′ + p_1 \times t \times \dfrac{-p_2}{d}$

其中只有 $t$ 是自变量,该式子等价于:

$x \equiv r_1 + p_1 \times k_1′ \pmod{|p_1 \times \dfrac{-p_2}{d}|}$

其中 $|p_1 \times \dfrac{-p_2}{d}| = \operatorname{lcm}(p_1,p_2)$

通过一系列的变换,我们成功将两个线性同余方程合并为一个线性同余方程。

依次类推,即可求出最后的一个线性同余方程,随后求出原线性同余方程组的解。

在计算过程中注意取模,将可取模的都取模以避免溢出。


#include <cstdio>
#define ll long long
const int maxn = 1e6 + 100;
int n;
ll p[maxn], r[maxn];
ll exgcd(ll a, ll b, ll& x, ll& y)
{
    if (b == 0)
    {
        x = 1, y = 0;
        return a;
    }
    ll g = exgcd(b, a % b, y, x);
    y -= (a / b) * x;
    return g;
}
ll mul(ll a, ll b, ll mod)
{
    ll res = 0;
    while (b > 0)
    {
        if (b & 1) res = (res + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return ((res % mod) + mod) % mod;
}
inline ll exCRT()
{
    ll lcm = p[1], lastr = r[1], k1, k2;
    for (int i = 2; i <= n; ++i)
    {
        ll a = lcm, b = p[i], c = (((r[i] - lastr) % p[i]) + p[i]) % p[i];
        ll g = exgcd(a, b, k1, k2);
        ll mod = p[i] / g;
        k1 = ((k1 % mod) + mod) % mod;
        k1 = mul(k1, c / g, mod);

        lastr += k1 * lcm;
        lcm = lcm / g * p[i];
        lastr = ((lastr % lcm) + lcm) % lcm;
    }
    return lastr;
}
int main()
{
    scanf("%d", &n);
    for (int i = 1; i <= n; ++i) scanf("%lld %lld", p + i, r + i);
    printf("%lld\n", exCRT());
    return 0;
}

BSGS

全名为 BabyStepGiantStep 的算法,用于求解离散对数问题。

$$A^x \equiv B \pmod{C}$$

其中 $A,C$ 互质。

考虑暴力枚举,$x = 0$ 时 $A^x \equiv 1$,而由费马小定理,$A^{C-1} \equiv 1$,往后就出现重复。那么只需要枚举 $[0,C-1)$ 即可。

而利用分块思想可以简化计算。

设 $x = i \times \sqrt{C} – j,i \in [1,\sqrt{C}],j \in [0,\sqrt{C}]$

$A^{i \times \sqrt{C} – j} \equiv B \pmod{C}$

$\dfrac{A^{i \times \sqrt{C}}}{A^j} \equiv B \pmod {C}$

$\dfrac{A^{i \times \sqrt{C}}}{A^j \times B} \equiv 1 \pmod{C}$

即 $\bmod C$ 意义下,$A^{i \times \sqrt{C}} = A^j \times B$

预处理出所有 $A^j \times B$ 的值对应的 $j$,枚举 $i$ 即可快速查询到对应的 $j$,从而计算出解。

可以利用哈希表实现。


#include <cstdio>
#include <list>
#include <cmath>
#include <cstring>
#define ll long long
ll C,A,B;
struct BetterListHashTable
{
    static const int mod = 1000103;
    static const int maxn = 1e6 + 100;
    struct node
    {
        int next, key, val;
    } A[maxn + 10];
    int head[mod + 100],tot;
    int hash(int x)
    {
        int res = x % mod;
        if (res < 0) res += mod;
        return res;
    }
    int & operator [](int key)
    {
        int pos = hash(key);
        for (int p = head[pos]; p; p = A[p].next)
            if (A[p].key == key) return A[p].val;
        A[++tot].next = head[pos], A[tot].key = key, A[tot].val = 0, head[pos] = tot;
        return A[tot].val;
    }
    int* find(int key)
    {
        int pos = hash(key);
        for(int p = head[pos];p;p=A[p].next)
            if(A[p].key == key) return &A[p].val;
        return NULL;
    }
}HT;
inline ll fastpow(ll x,ll k)
{
    ll res = 1;
    while(k)
    {
        if(k & 1) res = (res * x) % C;
        x = (x * x) % C;
        k >>= 1;
    }
    return res % C;
}
int main()
{
    scanf("%lld %lld %lld",&C,&A,&B);
    ll sqrtc = std::ceil(std::sqrt(C));
    ll t = B % C;
    for(int j = 0;j<=sqrtc;++j)
    {
        HT[t] = j;
        t = t * A % C;
    }
    ll sq = fastpow(A,sqrtc),ans;
    t = sq;
    for(int i = 1;i<=sqrtc;++i)
    {
        int *ans = HT.find(t);
        if(ans != NULL)
        {
            printf("%lld\n",i * sqrtc - *ans);
            return 0;
        }
        t = (t * sq) % C;
    }
    puts("no solution");
    return 0;
}

例题:多少个 $1$?

给整数 $K$ 与质数 $m$,求最小的 $N$ 使得 $11\cdots111 (N个1) \equiv K \pmod{m}$

两边同时 $\times 9 + 1$,可以化为:

$10^N \equiv 9K + 1 \pmod{m}$

然后套上板子即可。

exBSGS

$$A^x \equiv B \pmod{P}$$

然而当 $A,P$ 不互质时,就需要使用 EXBSGS.

设 $d_1 = \gcd(A,P)$,若 $d_1 \nmid B$,原方程无解。

证明:

证明

将 $A^{x-1},k$ 看做自变量,则当 $\gcd(A,P) \mid B$ 时该二元一次不定方程有解。


推式子

若 $\gcd(A,\dfrac{P}{d_1}) \neq 1$ ,则继续除。

令 $D = \prod \limits_{i=1}^m d_i$

推式子

此时 $\gcd(A,\dfrac{P}{D}) = 1$,可以使用 BSGS 求解。

随后加上 $m$ 即可。


#include <cstdio>
#include <cmath>
#include <cstring>
#define ll long long
ll fastpow(ll x,ll k,ll mod)
{
    ll res = 1;
    for (; k; k >>= 1)
    {
        if (k & 1) res = (res * x) % mod;
        x = (x * x) % mod;
    }
    return ((res % mod) + mod) % mod;
}
struct BetterListHashTable
{
    static const int mod = 100103;
    static const int maxn = 1e6 + 100;
    struct node
    {
        int next, key, val;
    } A[maxn + 10];
    int head[mod + 100], tot;
    void clear() { memset(head, 0, sizeof(head)), tot = 0; }
    int hash(int x)
    {
        int res = x % mod;
        if (res < 0) res += mod;
        return res;
    }
    int& operator[](int key)
    {
        int pos = hash(key);
        for (int p = head[pos]; p; p = A[p].next)
            if (A[p].key == key) return A[p].val;
        A[++tot].next = head[pos], A[tot].key = key, A[tot].val = 0, head[pos] = tot;
        return A[tot].val;
    }
    int* find(int key)
    {
        int pos = hash(key);
        for (int p = head[pos]; p; p = A[p].next)
            if (A[p].key == key) return &A[p].val;
        return NULL;
    }
} HT;
inline ll BSGS(ll A, ll B, ll P)
{
    //A^x = B (mod P) (A,P) = 1
    HT.clear();
    ll t = B;
    ll sqrtp = std::ceil(std::sqrt(P));
    for (int i = 0; i <= sqrtp; ++i) HT[t] = i, t = (t * A) % P;
    ll sq = fastpow(A, sqrtp, P);
    t = sq;
    for (int i = 1; i <= sqrtp; ++i)
    {
        int* ans = HT.find(t);
        if (ans != NULL) return i * sqrtp - *ans;
        t = (t * sq) % P;
    }
    return -1;
}
ll gcd(ll a, ll b) { return b == 0 ? a : gcd(b, a % b); }
void exgcd(ll a, ll b, ll& x, ll& y)
{
    if (b == 0) return (void)(x = 1, y = 0);
    exgcd(b, a % b, y, x), y -= (a / b) * x;
}
ll inv(ll x, ll P)
{
    //xinv + Pk = 1
    ll k1, k2;
    exgcd(x, P, k1, k2);
    k1 = ((k1 % P) + P) % P;
    return k1;
}
inline ll exBSGS(ll A,ll B,ll P)
{
    ll m = 0,D = 1,g;
    if(B == 1) return 0;
    while((g = gcd(A,P)) != 1)
    {
        if (B % g) return -1;
        P /= g, ++m;
    }
    B = (B * inv(fastpow(A, m, P), P)) % P;  //inv(a^m)
    ll res = BSGS(A, B, P);
    if (res == -1) return -1;
    return res + m;
}

例题:【模板】拓展BSGS

本文标题:[笔记]同余方程
本文作者:Clouder
本文链接: https://www.codein.icu/congruence/
转载请标明。
暂无评论

发送评论 编辑评论

|´・ω・)ノ
ヾ(≧∇≦*)ゝ
(☆ω☆)
(╯‵□′)╯︵┴─┴
 ̄﹃ ̄
(/ω\)
∠( ᐛ 」∠)_
(๑•̀ㅁ•́ฅ)
→_→
୧(๑•̀⌄•́๑)૭
٩(ˊᗜˋ*)و
(ノ°ο°)ノ
(´இ皿இ`)
⌇●﹏●⌇
(ฅ´ω`ฅ)
(╯°A°)╯︵○○○
φ( ̄∇ ̄o)
ヾ(´・ ・`。)ノ"
( ง ᵒ̌皿ᵒ̌)ง⁼³₌₃
(ó﹏ò。)
Σ(っ °Д °;)っ
( ,,´・ω・)ノ"(´っω・`。)
╮(╯▽╰)╭
o(*////▽////*)q
>﹏<
( ๑´•ω•) "(ㆆᴗㆆ)
😂
😀
😅
😊
🙂
🙃
😌
😍
😘
😜
😝
😏
😒
🙄
😳
😡
😔
😫
😱
😭
💩
👻
🙌
🖕
👍
👫
👬
👭
🌚
🌝
🙈
💊
😶
🙏
🍦
🍉
😣
Source: github.com/k4yt3x/flowerhd
颜文字
Emoji
小恐龙
花!
上一篇
下一篇