数论——线性同余方程、乘法逆元
众所周知:
说明
除非特殊说明,以下提到的 exgcd
函数均定义为:
// ax + by = gcd(a, b)
ll exgcd(ll a, ll b, ll &x, ll &y, ll d = 0)
{
if (b == 0) x = 1, y = 0, d = a;
else d = exgcd(b, a % b, y, x), y -= a / b * x;
return d;
}
以下提到的 quick_pow
函数均定义为:
// pow
int quick_pow(int a, int b, int p, int res = 1)
{
for (; b; b >>= 1)
b & 1 ? res = res * a % p, a = a *a % p : a = a * a % p;
return res;
}
如果爆 int 请自行开 long long 或边读边模,或高精度处理。
同余
定义
若 \(a \bmod m = b \bmod m\),则称 \(a\) 与 \(b\) 关于模 \(m\) 同余,记为 \(a \equiv b \pmod m\).
同余的性质
- 反身性:\(a \equiv a \pmod m\);
- 对称性:若 \(a \equiv b \pmod m\),则 \(b \equiv a \pmod m\);
- 传递性:若 \(a \equiv b \pmod m\)、\(b \equiv c \pmod m\),则 \(a \equiv c \pmod m\);
- 同余式相加:若 \(a \equiv b \pmod m\)、\(c \equiv d \pmod m\),则 \(a \pm c \equiv b \pm d \pmod m\);
- 同余式相乘:若 \(a \equiv b \pmod m\)、\(c \equiv d \pmod m\),则 \(a \times c \equiv b \times d \pmod m\);
- 乘方:若 \(a \equiv b \pmod m\),则 \(a^k \equiv b^k \pmod m\);
- 除法1:若 \(ka \equiv kb \pmod{km}\),则 \(a \equiv b \pmod m\);
- 除法2:若 \(ka \equiv kb \pmod m\),则 \(a \equiv b \pmod{m / \gcd(k, m)}\);
线性同余方程
形式
关于 \(x\) 的方程,形如 \(ax \equiv n \pmod b\),则称之为线性同余方程(Linear Congruence Equation)。
一般要求求出特解,或 \(x \in [0, b - 1]\) 的通解。
求解方法
方程 \(ax \equiv n \pmod b\) 可以理解为 \(ax + by = n\),其中 \(y\) 为一个整数。
-
证明如下:
因为 \(ax + by = n\),
所以 \((ax + by) \bmod b = n \bmod b\),即 \(ax \bmod b + by \bmod b = n \bmod b\),
因为 \(by \bmod b = 0\),所以 \(ax \bmod b = n \bmod b\),
转换为同余方程的形式就是 \(ax \equiv n \pmod b\).
因此原方程转化为 \(ax + by = n\),接下来就是扩展欧几里得算法的事情了;
详见:https://www.cnblogs.com/RainPPR/p/gcd-exgcd.html
解的判断
扩展欧几里得算法只能求解 \(ax + by = \gcd(a, b)\) 的情况,
所以只有当 \(n = k \times \gcd(a, b)\) ,其中 \(k\) 为任意正整数,才可以用扩展欧几里得算法求解。
-
证明如下:
可以求出一组 \(x_0\)、\(y_0\),使得 \(ax_0 + by_0 = \gcd(a, b)\),
等式两边同时乘以 \(k\),便得到:\(akx_0 + bky_0 = k \times \gcd(a, b)\),因此可得到 \(\left\{\begin{matrix} x = kx_0 \\ y = ky_0\end{matrix}\right.\),
此时便有 \(ax + by = n\).
特解到通解
下面假设有解:
我们已经将 \(ax \equiv n \pmod b\) 转化为 \(ax + by = n\),并通过扩展欧几里得算法解出来一个通解 \(x_0\).
然后请看:https://www.cnblogs.com/RainPPR/p/gcd-exgcd.html
设 \(t = \dfrac{\operatorname{lcm}(a, b)}{a}\),则有通解 \(x = x_0 + kt\),其中 \(k \in \mathbb{Z}\)。
特殊化的线性同余方程
考虑方程 \(ax \equiv 1 \pmod b\),也就是上面的 \(n = 1\)。
此时存在解的条件为 \(k \times \gcd(a, b) = n = 1\),也就是 \(k = 1\) 时的情况:
\(\gcd(a, b) = 1\),即 \(a\) 与 \(b\) 互质,这样就可以用扩展欧几里得算法求解 \(ax + by = 1\) 了。
所以此时的通解公式也可以化简为 \(x = x_0 + kb\):
证明:\(t = \dfrac{\operatorname{lcm}(a, b)}{a}\),而 \(\operatorname{lcm}(a, b) = ab\),所以就有 \(t = b\)。
代码实现
// ax = 1 (mod b)
int lieu1(int a, int b)
{
int x, y;
int d = exgcd(a, b, x, y);
if (d != 1)
return -1;
return (x % b + b) % b;
}
// ax = n (mod b)
int lieu(int a, int b, int n)
{
int x, y;
int d = exgcd(a, b, x, y);
if (n % d != 0)
return -1;
int t = b / d;
return (x % t + t) % t;
}
乘法逆元
有理数取模
加减法:\((a \pm b) \bmod p = (a \bmod p \pm b \bmod p) \bmod p\).
乘法:\((a \times b) \bmod p = (a \bmod p \times b \bmod p) \bmod p\).
那除法呢?举例可知 \(\dfrac{a}{b} \bmod p\) 不一定等于 \(\dfrac{a \bmod p}{b \bmod p}\).
如 \(\dfrac{10}{2} \bmod 3 = 5 \bmod 3 = 2\),而 \(\dfrac{10 \bmod 3}{2 \bmod 3} = \dfrac{1}{2}\).
乘法逆元的定义
若 \(\dfrac{a}{b} \bmod p = (a \times x) \bmod p\),
则称 \(x\) 为 \(b\) 的模 \(p\) 意义下的乘法逆元(或 \(x\) 为 \(b \bmod p\) 的逆元),记作 \(x = b^{-1}\)
思路
根据 \(\dfrac{a}{b} \bmod p = (a \times x) \bmod p\) 可以写出同余方程:\(\dfrac{a}{b} \equiv a \times x \pmod{p}\)
两边同时乘以 \(\dfrac{b}{a}\) 可以得到:\(bx \equiv 1 \pmod{p}\);或者可以理解为 \(x\) 在模 \(p\) 意义下等价于 \(\dfrac{1}{b}\)。
转化一下就是 \(xb + kp = 1\),而 \(xb + kp = \gcd(b, p)\),
因此逆元并不是普遍存在的,条件是 \(\gcd(b, p) = 1\),也就是 \(b\) 与 \(p\) 互质。
扩展欧几里得算法求逆元
上面已经得到了 \(bx \equiv 1 \pmod{p}\) 及 \(xb + kp = 1\),而这就是上面讲到的特殊化的线性同余方程,可以使用扩展欧几里得算法求逆元。
详见上面:线性同余方程。
快速幂求逆元
前置知识:快速幂、费马小定理
若 \(p\) 为素数,\(\gcd(a, p) = 1\),则 \(a^{p - 1} \equiv 1 \pmod{p}\)。
证明见:https://oi-wiki.org/math/number-theory/fermat/
仅当 \(p\) 是质数时,即 \(\gcd(b, p) = 1\) 时,也可以用快速幂求逆元:
上面已得 \(bx \equiv 1 \pmod{p}\),
根据费马小定理,\(b^{p - 1} \equiv 1 \pmod p\),
可以转化为 \(b \times b^{p - 2} \equiv 1 \pmod p\),
而我们要求的是 \(bx \equiv 1 \pmod p\)。
因此可得 \(x = b^{p - 2}\)。
代码实现
// s1: exgcd
int inv1(int a, const int p) {
int x, y;
exgcd(a, p, x, y);
return (x % p + p) % p;
}
// s2: pow
int inv2(int a, const int p) {
return quick_pow(a, p - 2, p);
}
线性求逆元
线性求任意 n 个数的逆元
给定长度为 \(n\) 的序列 \(a\)(\(1 \le a_i < p\)),求序列每个数的逆元。
- \(a_i\),表示原序列,即给定的序列;
- \(\displaystyle s_i = \prod_{i = 1}^n a_i\),表示原序列的前缀积。
- \(inv_i = {a_i}^{-1}\),表示原序列的乘法逆元,即待求的序列;
- \(\displaystyle sv_i = {s_i}^{-1} = \prod_{i = 1}^n sv_i\),表示原序列前缀积的乘法逆元,根据逆元性质也等于原序列乘法逆元的前缀积。
- 计算给定序列 \(a_i\) 的前缀积,记为 \(s_i\);
- 使用快速幂或扩展欧几里得法计算 \(s_n\) 的逆元,记为 \(sv_n\);
- 因为 \(sv_n\) 是 \(n\) 个数的积的逆元,所以当我们把它乘上 \(a_n\) 时,就会和 \(a_n\) 的逆元抵消;这样就得到了 \(a_1\) 到 \(a_{n - 1}\) 的积逆元,记为 \(sv_{n - 1}\);
- 同理我们可以依次计算出所有的 \(sv_i\),于是 \({a_i}^{-1}\) 就可以用 \(s_{i - 1} \times sv_i\) 求得。
所以我们就在 \(O(n + \log p)\) 的时间内计算出了 \(n\) 个数的逆元。
// 计算前缀积
s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % p;
// 计算全部乘法逆元的前缀积
sv[n] = quick_pow(s[n], p - 2, p);
// 递推前缀积、求序列的乘法逆元
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % p
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;
来自 OI-Wiki 的代码
s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * a[i] % p;
sv[n] = qpow(s[n], p - 2);
// 当然这里也可以用 exgcd 来求逆元,视个人喜好而定.
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * a[i] % p;
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;
特化:线性求 1 ~ n 的逆元
即原序列 \(a_i = i\),此时有更加快速的方法,但是这里不讲(见 OI-Wiki 内)。
我们在此就简化原程序。
s[0] = 1;
for (int i = 1; i <= n; ++i) s[i] = s[i - 1] * i % p;
sv[n] = quick_pow(s[n], p - 2, p);
for (int i = n; i >= 1; --i) sv[i - 1] = sv[i] * i % p;
for (int i = 1; i <= n; ++i) inv[i] = sv[i] * s[i - 1] % p;
例题
1、线性同余方程:P1082 同余方程
点击查看代码
#include <bits/stdc++.h>
#define rr read()
using namespace std;
typedef long long ll;
inline ll read()
{
ll num = 0, flag = 1;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-')
flag = -1;
for (; isdigit(ch); ch = getchar())
num = (num << 3) + (num << 1) + ch - '0';
return num * flag;
}
ll exgcd(ll a, ll b, ll &x, ll &y, ll d = 0)
{
if (b == 0) x = 1, y = 0, d = a;
else d = exgcd(b, a % b, y, x), y -= a / b * x;
return d;
}
int main()
{
ll a = rr, b = rr;
ll x = 0, y = 0;
exgcd(a, b, x, y);
printf("%lld", (x % b + b) % b);
return 0;
}
2、快速幂求逆元:P2613 有理数取余
点击查看代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD = 19260817;
inline ll read()
{
ll num = 0;
char ch = getchar();
while (!isdigit(ch))
ch = getchar();
for (; isdigit(ch); ch = getchar())
num = (num << 3) + (num << 1) + ch - '0', num %= MOD;
return num;
}
ll qpow(ll a, ll b, const ll p, ll res = 1)
{
for (; b; b >>= 1)
b & 1 ? res = res * a % p, a = a *a % p : a = a * a % p;
return res;
}
int main()
{
ll a = read(), b = read();
if (b == 0)
printf("Angry!\n"), exit(0);
ll res = a * qpow(b, MOD - 2, MOD) % MOD;
printf("%lld\n", res);
return 0;
}
3、线性求 1 ~ n 的逆元:P3811 模意义下的乘法逆元
点击查看代码
#include <bits/stdc++.h>
#define rr read()
using namespace std;
typedef long long ll;
const int N = 3e6 + 10;
inline ll read()
{
ll num = 0, flag = 1;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (num == -1)
flag = -1;
for (; isdigit(ch); ch = getchar())
num = (num << 3) + (num << 1) + ch - '0';
return num * flag;
}
ll s[N], sv[N];
ll qpow(ll a, ll b, const ll p, ll r = 1)
{
for (; b; b >>= 1)
b & 1 ? r = r * a % p, a = a * a % p : a = a * a % p;
return r;
}
int main()
{
const int n = rr;
const ll p = rr;
s[0] = 1;
for (int i = 1; i <= n; ++i)
s[i] = s[i - 1] * i % p;
sv[n] = qpow(s[n], p - 2, p);
for (int i = n; i; --i)
sv[i - 1] = sv[i] * i % p;
for (int i = 1; i <= n; ++i)
printf("%lld\n", sv[i] * s[i - 1] % p);
return 0;
}
4、线性求 n 数的逆元:P5431 模意义下的乘法逆元 2
求:\(\sum\limits_{i = 1}^n \frac{k^i}{a_i}\).
点击查看代码
#include <bits/stdc++.h>
#define rr read()
using namespace std;
typedef long long ll;
const int N = 5e6 + 10;
inline ll read()
{
ll num = 0, flag = 1;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (num == -1)
flag = -1;
for (; isdigit(ch); ch = getchar())
num = (num << 3) + (num << 1) + ch - '0';
return num * flag;
}
ll a[N];
ll s[N], sv[N];
ll qpow(ll a, ll b, const ll p, ll r = 1)
{
for (; b; b >>= 1)
b & 1 ? r = r * a % p, a = a *a % p : a = a * a % p;
return r;
}
int main()
{
const int n = rr;
const ll p = rr, k = rr;
s[0] = 1;
for (int i = 1; i <= n; ++i)
a[i] = rr, s[i] = s[i - 1] * a[i] % p;
sv[n] = qpow(s[n], p - 2, p);
for (int i = n; i; --i)
sv[i - 1] = sv[i] * a[i] % p;
ll res = 0, kt = k;
for (int i = 1; i <= n; ++i)
res = (res + kt * (sv[i] * s[i - 1] % p) % p) % p, kt = kt * k % p;
printf("%lld\n", res);
return 0;
}
Reference
[1] https://oi-wiki.org/math/number-theory/fermat/
[2] https://oi-wiki.org/math/number-theory/inverse/
[3] https://oi-wiki.org/math/number-theory/linear-equation/