Nickel_Angel's nest

Nickel_Angel's nest

The key to magic is inside of us.

[数学]求解一类方程的方法

posted on 2019-07-05 20:31:47 | under 学习笔记 |

在做了几天图论和一些奇怪的题后,突然发现寒假培训期间的拓展欧几里得算法等算法还未掌握,加上最近学习的BSGS算法,算是一个关于在数论中解一类方程问题的总结吧QwQ

(注:本文中的字母若无特殊说明均可视为整数)

1.线性同余方程

1.方程简介

定义:线性同余方程是形如 $ax \equiv c \pmod{b}$ 的方程(其中 $x$ 为未知数)。

引理:对于方程 $ax \equiv c \pmod{b}$,其等价于方程 $ax+by = c$,其中 $y = \lfloor \frac c b \rfloor - \lfloor \frac {ax} {b} \rfloor$

证明:

将原方程的同余符号转化为等号得:

$$ax\,mod\,c = b\,mod\,c$$

由模的定义($a\,mod\,b = a - b \times \lfloor \frac a b \rfloor$)得:

$$ax - b \times \lfloor \frac {ax} {b} \rfloor = c - b \times \lfloor \frac c b \rfloor$$

移项并合并同类项得:

$$ax + (\lfloor \frac c b \rfloor - \lfloor \frac {ax} {b} \rfloor)b = c$$

令 $y = \lfloor \frac c b \rfloor - \lfloor \frac {ax} {b} \rfloor$,即可得 $ax + by = c$,从而得证。

2.判断方程有整数解的条件

引理:方程 $ax + by = c$ 有整数解的充要条件是 $gcd(a,b) \mid c$

在证明这个命题之前,先证明如下定理:

引理(裴蜀定理):对于 $\forall a, b \in Z,\exists x,y \in Z$ 使得 $ax + by = gcd(a,b)$

证:

首先,显然有 $gcd(a,b) \mid a$,$gcd(a,b) \mid b$ 。由整除的性质可知 $\forall x,y \in Z, gcd(a,b) \mid (ax + by)$ 。设 $s$ 为 $ax + by$ 的最小正值,故有 $gcd(a,b) \mid s$。

设 $k = \lfloor \frac {a} {s} \rfloor, r = a\,mod\,s$,则有 $r = a - k(ax + by) = a(1 - kx) + b(-ky)$,发现 $r$ 也为 $(a,b)$ 的线性组合,而 $r$ 又在模 $s$ 剩余系中,故 $r \in [0, s - 1]$,又由于 $s$ 为最小正值,故 $r = 0$。

所以 $s \mid a$,同理 $s \mid b$,根据两个数的公约数必为这两个数的最大公约数的约数,有 $s \mid gcd(a,b)$,而上文已证明$gcd(a,b) \mid s$ ,故 $s = gcd(a,b)$ 。

进而得证,故 $ax + by = gcd(a,b)$ 这个方程一定有解……

考虑将 $ax + by = gcd(a,b)$ 变为 $ax+by = c$:

为了区分这两个方程,我们将前者改为 $a'x + b'y = gcd(a',b')$

令 $k = gcd(a',b')$,方程两边同除 $k$,得:

$$\frac {a'} {k} x + \frac {b'} {k} y = 1$$

方程两边同乘$c$得:

$$\frac {a'c} {k} x + \frac {b'c} {k} y = c$$

只要通过换元,即令 $a = \frac {a'c} {k},b = \frac {b'c} {k}$,即可将方程转化为 $ax + by = c$ 的形式

考虑到 $gcd(\frac {a'c} {k},\frac {b'c} {k}) = c = gcd(a, b)$,且对于一般的方程:$atx + bty = c'(c' = ct, t \in Z)$

故只有 $c$ 为 $gcd(a,b)$ 倍数时原方程有解,即 $gcd(a,b) \mid c$,从而得证。

3.解方程的方法(拓展欧几里得算法)

我们可在“判断方程有解的条件”的证明中得到启发:解方程 $ax + by = c$,不妨先解方程 $a'x + b'y = gcd(a',b')$

根据辗转相除法求两个数的最大公约数可得:$gcd(a,b) = gcd(b, a\,mod\,b)$

故有:

$$\begin{aligned}a'x + b'y & = gcd(a',b') \\[2ex] & = gcd(b', a'\,mod\,b') \\[2ex] & = b'x + (a'\,mod\,b')y \\[2ex] & = b'x + (a' - \frac {a'} {b'} \times b')y \\[2ex] & = a'y + b'(x - \frac {a'} {b'} \times y) \\[2ex] \end{aligned}$$

令 $x' = y,y' = x - \frac {a'} {b'} \times y$,我们得到了一个 $a'x' + b'y' = gcd(a',b')$ 的方程……

不难发现这个方程可以像上文那样迭代下去,直到 $b' = 0$ 时,得出 $x' = 1, y' = 0$ 然后递归回去得出原方程的解。

假设我们得到了 $a'x + b'y = gcd(a',b')$ 的一组解 $x_0,y_0$,从上文得知 $a = \frac {a'c} {gcd(a',b')},b = \frac {b'c} { gcd(a',b')}$,故原方程的解为 $x = \frac {x_0c} {gcd(a',b')}$

Code:

int exgcd(int a, int b, int &x, int &y)
//使用方法:初始调用本函数时调用原方程的a和b 
{
    if (b == 0)//递归终点
    {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);//注意这里y与x是反过来写的
    y -= a / b * x;
    return d;
}

4.例题

P1082 同余方程

P1516 青蛙的约会

2.线性同余方程组

1.方程简介

定义:线性同余方程组为形如:

$$\begin{cases} x \equiv b_1 \pmod{a_1} \\ x \equiv b_2 \pmod{a_2} \\ \vdots \\ x \equiv b_n \pmod{a_n} \end{cases}$$

的方程组。

2.解方程的方法

由于我们一般求的是最小正整数解,故这里只讨论最小整数数解的求法。

我们先来讨论 $a_i$ 彼此互质的情况(我们称这种算法为 CRT (Chinese Remainder Theorem),即中国剩余定理(算法导论中将其翻译为“中国余数定理”))。

我们先来证明该方程组在 $[0, \prod_{i = 1}^n a_i)$ 区间内有唯一解。

引理:(中国剩余定理)若两个互质的正整数 $n, m$,若两个整数 $a,b$ 满足 $0 \leq a < n, 0 \leq b < m$,存在且仅存在一个 $x$ 满足 $0 \leq x < nm$ 且 $x = a \text{ } mod \text{ } n,x = b \text{ } mod \text{ } m$ 。

证:

考虑反证法,设存在两个不相等的整数 $0 < x_1, x_2 \leq nm$,使得 $x_1, x_2$ 关于 $n, m$ 同余。

不妨设 $x_1 < x_2,k = x_2 - x_1$,显然有 $nm \mid k$,故当 $k$ 为最小值时 $k_{min} = lcm(n,m)$,又因为 $n,m$ 互质,故 $lcm(n, m) = nm$,进而有 $k_{min} = nm$ 。

故 $x_2 = x_1 + k \geq x_1 + nm \geq nm$,与假设矛盾,故原命题成立。

由这个定理数学归纳可知上文的方程组在 $[0, \prod_{i = 1}^n a_i)$ 区间内有唯一解。

考虑构造一个解满足此方程组,我们通过一些玄学的方法构造出了一个解为:

$$x = \sum_{i = 1}^n \limits k_i \frac {\prod_{j = 1}^n \limits a_j} {a_i}$$

其中 $k_i$ 为满足 $k_i \frac {\prod_{j = 1}^n a_j} {a_i} \equiv b_i \pmod{a_i}$ 的最小非负整数。

我们可以将 $x$ 代入原方程,根据取模运算对加法的分配律可知,$x$ 为原方程的一个解。

考虑如何计算 $k_i$,为了方便,我们设 $m_i = \frac {\prod_{j = 1}^n a_j} {a_i}$ 。注意到 $a_i$ 彼此互质,则 $m_i$ 一定与 $a_i$ 互质。故在模 $a_i$ 意义下一定存在 $m_i$ 的乘法逆元,我们设其为 $t_i$,根据乘法逆元的定义,有 $t_im_i \equiv 1 \pmod{a_i}$,再根据同余的同乘性,有 $b_it_im_i \equiv b_i \pmod{a_i}$,故 $b_it_i \equiv k_i \pmod{a_i}$,由于 $k_i$ 的最小性,$k_i = b_it_i \text{ } mod \text{ } a_i$ 。

故原方程的一个通解为:

$$x = \sum_{i = 1}^n \limits (b_it_i \text{ } mod \text{ } a_i) \times m_i$$

显然我们可以通过拓展欧几里得算法,对于每一个 $m_i$ 解线性同余方程得到所有的 $t_i$,进而轻易得到方程的解。

Code:

typedef long long ll;

inline ll CRT(int *a, int *b)
{
    ll A = 1, x, y, ans = 0;
    for (int i = 1; i <= n; ++i)
    {
        A *= a[i];
        b[i] = (b[i] % a[i] + a[i]) % a[i];//防止 b[i] 出现负数
    }
    for (int i = 1; i <= n; ++i)
    {
        ll res = A / a[i];//求 m[i]
        exgcd(res, a[i], x, y);//求 m[i] 逆元
        x = (x % a[i] + a[i]) % a[i];
        ans = (ans + mul(mul(x, res, A), b[i], A)) % A;
        //mul 为龟速乘防爆 long long
    }
    return (ans + A) % A;
}

然而在不保证 $a_i$ 两两互质时,又是另一种情况。

由于 $a_i$ 之间可能不互质,故我们可能找不到某个 $m_i$ 在模 $a_i$ 意义下的逆元,故我们得另辟蹊径。

不妨从小规模的问题开始考虑,考虑对于仅由两个线性同余方程组成的方程组:

$$\begin{cases} x \equiv b_1 \pmod{a_1} \\ x \equiv b_2 \pmod{a_2} \end{cases}$$

我们考虑其本质是什么,我们可以将其化成等式的形式:

$$\begin{cases} x = b_1 + k_1a_1 \\ x = b_2 + k_2a_2 \end{cases}$$

不难发现 $b_1 + k_1a_1 = b_2 + k_2a_2$ 。

故 $k_1a_1 - k_2a_2 = b_1 - b_2$ 。

由于 $a_1,a_2,b_1,b_2$ 已知,故我们可以将其看做一个线性同余方程,用拓展欧几里得算法即可求出 $k_1, k_2$ 。

而由该线性同余方程有解条件可知 $gcd(a_1, a_2) \mid b_1 - b_2$ 。

设 $x_0 = b_1 + k_1a_1$,而由于 $lcm(a_1, a_2) \text{ } mod \text{ } a_1 = 0, lcm(a_1, a_2) \text{ } mod \text{ } a_2 = 0$ ,故原方程组有一个通解为 $x = x_0 + k \times lcm(a_1, a_2)$,我们可以将两个方程组合并得 $x \equiv x_0 + lcm(a_1,a_2) \pmod{lcm(a_1,a_2)}$ 。

(这里要加上一个 $lcm(a_1,a_2)$ 是为了保证答案非负)

故对于 $n$ 个方程,我们顺次合并两个方程即可。

typedef long long ll;

inline ll exCRT(ll *a, ll *b, int n)
{
    ll ans = b[1];//注意初始化
    ll x, y, m = a[1];//这里的 m 是当前已经合并了的所有方程的 a[i] 的 lcm
    for (int i = 2; i <= n; ++i)
    {
        ll B = ((b[i] - ans) % a[i] + a[i]) % a[i];//求两个要合并的方程 b 之差
        ll d = exgcd(m, a[i], x, y);
        x = mul(x, B / d, a[i]);//求出线性同余方程的解
        ans += m * x;//求出 x0
        m *= a[i] / d;//利用两个数的乘积等于这两个数 gcd 与 lcm 的乘积的原理更新 lcm
        ans = (ans + m) % m;//合并方程的解
    }
    return ans;
}

3.高次不定方程

1.方程简介

定义:高次不定方程是形如 $A^x \equiv B \pmod{p}$ 的方程(其中 $x$ 为未知数)

除此之外就似乎没什么了

2.判断方程有整数解的条件

分为两种情况,$A, p$ 互质和 $A, p$ 不互质的情况……

先说当 $A, p$ 互质的情况

(1)当 $p \mid A$ 且 $B \neq 0$ 时无解

证明:

由于 $p \mid A$,故$A^x \equiv 0 \pmod{p}$

而 $B \neq 0$,故无解……

(2)当通过枚举等方法发现其在 $[0,p - 1]$ 内无整数解,则此方程无整数解(关于当 $p$ 为质数时,方程一定在 $[0,p - 1]$有一个整数解的证明见下文)我能怎么办,我也很绝望啊 QAQ

$A, p$ 不互质时,有解的条件将在下文提出

3.解方程的方法

分为两种情况,$A, p$ 互质和 $A, p$ 不互质的情况…… 人类的本质是……

对于 $A, p$ 互质的情况,可以用 BSGS 算法来求解,而 $A, p$ 不互质情况则需用到 ExBSGS……


这里是对 $A, p$ 互质时的解法QwQ

(1)朴素的算法

在介绍新算法前,不妨考虑一下朴素的算法。

最容易想到的是枚举正整数 $x$,将其带入方程,直到找到一个 $x$ 使方程成立,我们就找到了一个解……

没错就是这么简单粗暴 QwQ

我们现在来估计一下其时间复杂度,不妨考虑一下解的范围问题。

引理:(费马小定理)若 $p$ 为质数,且整数 $a$ 满足 $p \nmid a$,则有 $a^{p-1} \equiv 1 \pmod{p}$

在证明这条定理之前,我们先再看两个引理:

引理1:对于整数 $a,b,c$,$d$ 为正整数,其满足 $gcd(c,d) = 1,ac \equiv bc \pmod{d}$,则有 $a \equiv b \pmod{d}$

证明:

$\because ac \equiv bc \pmod{d}$

$\therefore ac - bc \equiv 0 \pmod{d}$

$\therefore (a - b)c \equiv 0 \pmod{d}$

其等价于:

$$(a - b)c = kd(k \in Z)$$

由于 $gcd(c,d) = 1$,即$c,d$互质,故$c,d$无公共因子,进而得知要使等式成立,必有 $c \mid k$,即 $\frac k c$ 为整数

故等式两边同除 $c$ 得:

$$a - b = \frac {k} {c} d$$

考虑到 $\frac k c$ 为整数,故方程可重新写为:

$$a - b \equiv 0 \pmod{d}$$

即:

$$a \equiv b \pmod{d}$$

进而得证。

引理2:设 $m$ 为大于 1 的整数,$b$ 是一个整数且满足 $gcd(b,m) = 1$;若 $a_1,a_2,\cdots,a_m$ 是模 $m$ 的一个完全剩余系,则有正整数 $b$ 使得 $a_1b,a_2b,\cdots,a_mb$ 也构成模 $m$ 的一个完全剩余系

证明:

考虑反证法,若存在两个整数 $i,j \in [1,m]$ 且 $i \neq j$ 使得 $a_ib$ 与 $a_jb$ 在模 $m$ 意义下同余,即 $a_ib \equiv a_jb \pmod{m}$,根据引理 1 则有 $a_i \equiv a_j \pmod{m}$,而这与 $a_1,a_2,\cdots,a_m$ 是模 $m$ 的一个完全剩余系矛盾,故 $a_1b,a_2b,\cdots,a_mb$ 也构成模 $m$ 的一个完全剩余系,进而得证。

在证明了两个引理后,我们接下来来证明费马小定理:

考虑构造模 $p$ 的完全剩余系: $P = \{0,1,2,\cdots,p - 1\}$

由 $p$ 为质数且 $p \nmid a$ 得 $a,p$ 互质,即 $gcd(a,p) = 1$,由引理 2 得 $A = \{0, a, 2a,\cdots,(p - 1)a\}$ 也是模 $p$ 的一个完全剩余系;

根据完全剩余系的性质,可得:

$$1 \times 2 \times \cdots \times (p - 1) \equiv a \times 2a \times \cdots \times (p - 1)a \pmod{p}$$

即:

$$(p - 1)! \equiv (p - 1)! \times a^{p - 1} \pmod{p}$$

由于$p$为质数,由质数的定义可知,$p$ 与 $1,2,\cdots,p - 1$ 中所有数均互质,故 $p$ 与 $(p - 1)!$ 互质;

故由引理 1 得:

$$a^{p - 1} \equiv 1 \pmod{p}$$

进而得证。

在证明了费马小定理后,不难发现,对于一个高次不定方程 $A^x \equiv B \pmod{p}$ 来说,当 $p$ 为质数时,有:

$$A^{p - 1} \equiv 1 \pmod{p}$$

而 $A^0 \equiv 1 \pmod{p}$

故 $A^0 \equiv A^{p - 1} \pmod{p}$

故在 $0,1,2,\cdots,p - 1$ 必定存在一个长度不超过$p$的循环节,方程一定在 $[0,p - 1]$ 有一个整数解

故朴素枚举的时间复杂度为 $O(p)$。

(2)BSGS(Baby-Step-Giant-Step) 算法

朴素枚举的 $O(p)$ 的时间复杂度在 $p$ 很大的时候效率较低,故考虑优化算法

不妨令 $x = am - b$,故原方程变为 $A^{am - b} \equiv B \pmod{p}$ 。

方程两边同乘 $A^b$ 得:

$$A^{am} \equiv A^bB \pmod{p}$$

然而还是考虑枚举……将 $m$ 固定(这样也同时固定了枚举的范围),先枚举 $b$,将枚举过程中将对应的 $A^bB$ 的值存进哈希表;后枚举 $a$,并同时通过查表找到一个对应的 $b$ 使得方程成立,不难发现,整数 $b$ 的枚举范围为 $[0,m]$,整数 $a$ 的枚举范围为 $[0,\frac {p} {m}]$,故当 $m = \sqrt{p}$ 时,时间复杂度最优,为$O(\sqrt{p})$。

故 BSGS 算法的时间复杂度为 $O(\sqrt{p})$。

Code:(由于自己写的哈希表常数太大,这里用 unordered_map 代替了哈希表)

unordered_map<int, int> V;

inline int BSGS(int A, int B, int p)
{
    if (A % p == 0 && B) return -1;//判断无解的情况,这里无解返回-1
    A %= p, B %= p;//小优化
    if (B == 1) return 0;//当B为1时,A的0次方恰为1,与B在模p意义下同余
    int m = floor(sqrt(1.0 * p)) + 1//注意sqrt()计算出的结果需上取整
    V.clear();//清空哈希表
    for (int i = 0, t = B; i < m; ++i, t = 1ll * t * A % p)
        V[t] = i;//枚举a,并将对应的值插入哈希表
    int j = power(A, m, p);//这里power为快速幂
    for (int i = 1, t = j; i <= m + 1; ++i, t = 1ll * t * j % p)//注意乘1ll防爆int
    {
        if (V.count(t)) return i * m - V[t];//枚举b,并查表
    }
    return -1;//枚举完还未找到解,可知方程无解
}

下面是 ExBSGS QwQ

当 $A, p$ 不互质时,考虑将同余式 $A^x \equiv B \pmod{p}$ 转化成等式 $A^x + kp = B$ 。

设 $d = gcd(A, p)$,将等式两边同除 $d$,式子将转化为 $\frac{A} {d} A^{x - 1} + \frac{p} {d} k = \frac {B} {d}$,令 $p' = \frac {p} {d}$,若此时 $A, p'$ 仍不互质,则等式两边继续同除 $d' = gcd(A, p')$,这样一直除到 $A, p_i$ 互质为止(为了方便,我们将等式两边第 $i$ 次同除 $p$ 的值设为 $p_i$)。

假设我们两边同除了 $y$ 次,设 $d_i$ 为第 $i$ 次两边同除,$gcd(A, p_{i - 1})$ 的值,则最终得到的式子为:

$$\frac {A^y} {\prod_{i = 1}^y \limits d_i} A^{x - y} \equiv \frac {B} {\prod_{i = 1}^y \limits d_i} \pmod{\frac{p} {\prod_{i = 1}^y \limits d_i}}$$

若同除过程中出现 $\frac {B} {\prod d_i}$ 不为整数的情况,则原方程无解。

我们不妨设 $C = \frac {A^y} {\prod_{i = 1}^y \limits d_i}, B' = \frac {B} {\prod_{i = 1}^y \limits d_i}, P' = \frac{p} {\prod_{i = 1}^y \limits d_i}, x' = x - y$ 。

故原方程被我们转化成了:

$$CA^{x'} \equiv B' \pmod{P'}$$

其中 $A, P'$ 互质,而 $C$ 不影响我们解方程,故直接套用上文的 BSGS 即可。

Code:

unordered_map<int, int> V;

inline int exBSGS(int A, int B, int p)
{
    A %= p, B %= p;
    if (B == 1) return 0;
    if (!B && !A) return 1;
    if (!A) return -1;
    if (!B)//特判 B == 0 的情况
    {
        int res = 0, d = gcd(A, p);
        while (d != 1)
        {
            ++res, p /= d;
            d = gcd(A, p);
            if (p == 1) return res;
        }
        return -1;
    }
    int res = 0, d = gcd(A, p), C = 1;
    while (d != 1)//计算 y(res)
    {
        if (B % d != 0) return -1;
        p /= d, B /= d, C = 1ll * A / d * C % p;
        ++res, d = gcd(A, p);
        if (C == B) return res;
    }
    V.clear();
    int m = floor(sqrt(p)) + 1;//最后用 BSGS 解方程
    for (int i = 0, t = B % p; i < m; ++i, t = 1ll * t * A % p)
        V[t] = i;
    int j = power(A, m, p);
    for (int i = 1, t = 1ll * j * C % p; i <= m; ++i, t = 1ll * t * j % p)
    //注意这里枚举时要多乘上一个 C
    {
        if (V.count(t)) return i * m - V[t] + res;//记得加上 y(res)
    }
    return -1;
}

4.例题

P3846 [TJOI2007]可爱的质数

P2485 [SDOI2011]计算器

P4195 【模板】exBSGS/Spoj3105 Mod

4.后记

本文部分参考了百度百科上的证明方法。