文章目录
前言
在考试中,经常会出现一类问题,它们不涉及很深的算法,却和数学息息相关。这样的问题通常难度不大,也不需要特别的数学知识,只要掌握简单的数理逻辑即可。——《算法笔记》胡凡
质数(素数)
质数又称素数,是指除了1和本身外,不能被其它的任何数整除。
质数的判定
如何判断给定的正整数n是否是质数?
时间复杂度: O ( s q r t ( n ) ) O(sqrt(n)) O(sqrt(n))
bool isPrime(int n) {
if (n < 2) return false;
for (int i = 2; i <= n / i; i++) // i <= n / i
if (n % i == 0) return false;
return true;
}
分解质因数
n
中有且仅有一个大于sqrt(n)
的质因子。(如果有两个,两者相乘将大于n
,明显不符合)
据此,我们可以先枚举小于sqrt(n)
的所有质因子
如果最后n
仍大于1,那么这个n
就是那个大于sqrt(n)
的质因子
最坏时间复杂度: O ( s q r t ( n ) ) O(sqrt(n)) O(sqrt(n))
最好时间复杂度:
O
(
l
o
g
(
n
)
)
O(log(n))
O(log(n)),(若n
为
2
k
2^k
2k,那么在i=2
时,所有因子就被分解出来了,所以时间应该是
l
o
g
(
2
k
)
=
k
log(2^k)=k
log(2k)=k)
void divide(int n) {
for (int i = 2; i <= n / i; i++)
if (n % i == 0) {
int s = 0;
while (n % i == 0) {
s++;
n /= i;
}
printf("%d %d\n", i, s);
}
if (n > 1) printf("%d %d", n, 1);
}
筛质数
筛选出从1-n
之间质数的个数
数据量小于等于 1 0 6 10^6 106时,以下两个算法效率不相上下
埃氏筛法思想
一个数是p
,如果2~p-1
中的质数都不是它的质因数的话,那么p
是质数(以此对算法进行优化)
朴素筛法
O ( n l o g ( n ) ) O(nlog(n)) O(nlog(n))
时间复杂度计算,当n=2时,执行n/2次;当n=3时,执行n/3次
n 2 + n 3 + n 4 + . . . + n n \frac{n}{2}+\frac{n}{3}+\frac{n}{4}+...+\frac{n}{n} 2n+3n+4n+...+nn
= n ∗ ( 1 2 + 1 3 + 1 4 + . . . + 1 n ) =n*(\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+...+\frac{1}{n}) =n∗(21+31+41+...+n1)
= n ∗ lim x → + ∞ 1 2 + 1 3 + 1 4 + . . . + 1 n =n*\lim_{x \to +\infty}{\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+...+\frac{1}{n}} =n∗limx→+∞21+31+41+...+n1
= n ∗ l n ( n ) < n ∗ l o g ( n ) = n l o g ( n ) =n*ln(n)<n*log(n)=nlog(n) =n∗ln(n)<n∗log(n)=nlog(n)
int primes[N];
bool st[N];
void get_primes(int n) {
for (int i = 2; i <= n; i++) {
if (!st[i]) primes[cnt++] = i;
for (int j = i + i; j <= n; j += i) st[j] = true;
}
}
cout << cnt << endl; // 输出质数个数
优化
O ( n l o g ( l o g ( n ) ) ) O(nlog(log(n))) O(nlog(log(n)))。相当于线性的时间复杂度
1~n
中有
n
l
n
(
n
)
\frac{n}{ln(n)}
ln(n)n个质数
int primes[N];
bool st[N];
void get_primes(int n) {
for (int i = 2; i <= n; i++) {
if (!st[i]) {
primes[cnt++] = i;
// 是质数的倍数的话,判定为合数
for (int j = i + i; j <= n; j += i) st[j] = true;
}
}
}
cout << cnt << endl; // 输出质数个数
线性筛法(*常用)
n
只会被最小质因子筛掉
-
i % primes[j] == 0
p r i m e s [ j ] primes[j] primes[j]一定是 i i i的最小质因子, p r i m e s [ j ] primes[j] primes[j]一定是 p r i m e s [ j ] ∗ i primes[j]*i primes[j]∗i的最小质因子
-
i % primes[j] != 0
p r i m e s [ j ] primes[j] primes[j]一定小于 i i i的所有质因子, p r i m e s [ j ] primes[j] primes[j]也一定是 p r i m e s [ j ] ∗ i primes[j]*i primes[j]∗i的最小质因子
int primes[N]; // 用于存储质数
bool st[N]; // false时为质数
void get_primes(int n ) {
for (int i = 2; i <= n; i++) {
if (!st[i]) primes[cnt++] = i;
for (int j = 0; primes[j] <= n / i; j++) {
st[primes[j] * i] = true;
if (i % primes[j] == 0) break; // primes[j]一定是i的最小质因子
}
}
}
约数
12的约数为1,12,3,4,2,6
试除法求约数
O ( s q r t ( n ) ) O(sqrt(n)) O(sqrt(n))
vector<int> get_divisors(int n) {
vector<int> res;
for (int i = 1; i <= n / i; i++) // 只枚举约数里最小的
if (n % i == 0) {
res.push_back(i);
if (i != n / i) res.push_back(n / i); // 避免出现5*5的情况
}
sort(res.begin(), res.end());
return res;
}
约数个数
若要求出N
的约数个数,首先将N
质因数分解,p
代表质数,
α
\alpha
α为p
的次数
N = p 1 α 1 + p 2 α 2 + p 3 α 3 + . . . + p n α n N=p_1^{\alpha_1}+p_2^{\alpha_2}+p_3^{\alpha_3}+...+p_n^{\alpha_n} N=p1α1+p2α2+p3α3+...+pnαn
然后求出每个质因数的约数个数,如
p 1 α 1 : p 1 0 , p 1 1 , . . . , p 1 α 1 − 1 , p 1 α 1 p_1^{\alpha_1}:p_1^0,p_1^1,...,p_1^{\alpha_1-1},p_1^{\alpha_1} p1α1:p10,p11,...,p1α1−1,p1α1
所以 p 1 α 1 p_1^{\alpha_1} p1α1共 α 1 + 1 \alpha_1+1 α1+1个约数,同理 p 2 α 2 p_2^{\alpha_2} p2α2共 α 2 + 1 \alpha_2+1 α2+1个约数……
根据乘法原则,N的约数个数为 ( α 1 + 1 ) ( α 2 + 1 ) ∗ . . . ∗ ( α n + 1 ) (\alpha_1+1)(\alpha_2+1)*...*(\alpha_n+1) (α1+1)(α2+1)∗...∗(αn+1)
约定 n n n个整数 a i a_i ai,请你输出这些数的乘积的约数个数,答案对 1 0 9 + 7 10^9+7 109+7取模
数据范围
1 ≤ n ≤ 100 1 \le n \le 100 1≤n≤100
1 ≤ a i ≤ 2 ∗ 1 0 9 1 \le a_i \le 2*10^9 1≤ai≤2∗109
输入格式
3
2
6
8
输出格式
12
#include <iostream>
#include <algorithm>
#include <unordered_map>
using namespace std;
typedef long long LL;
const int mod = 1e9 + 7;
int n;
int main() {
cin >> n;
unordered_map<int, int> primes;
while (n -- ) {
int x;
cin >> x;
// 分解质因数
for (int i = 2; i <= x/i; i++)
while (x % i == 0) {
x /= i;
primes[i]++;
}
// 记得收尾
if (x > 1) primes[x]++;
}
LL res = 1;
// 代约数个数公式
for (auto prime: primes) res = res * (prime.second + 1) % mod;
cout << res << endl;
return 0;
}
约数之和
若要求出N
的约数之和,首先将N
质因数分解,p
代表质数,
α
\alpha
α为p
的次数
N = p 1 α 1 + p 2 α 2 + p 3 α 3 + . . . + p n α n N=p_1^{\alpha_1}+p_2^{\alpha_2}+p_3^{\alpha_3}+...+p_n^{\alpha_n} N=p1α1+p2α2+p3α3+...+pnαn
约数之和公式为
( p 1 0 + p 1 1 + . . . + p 1 α 1 ) ( p 2 0 + p 2 1 + . . . + p 2 α 2 ) . . . ( p n 0 + p n 1 + . . . + p n α n ) (p_1^0+p_1^1+...+p_1^{\alpha_1})(p_2^0+p_2^1+...+p_2^{\alpha_2})...(p_n^0+p_n^1+...+p_n^{\alpha_n}) (p10+p11+...+p1α1)(p20+p21+...+p2α2)...(pn0+pn1+...+pnαn)
至于为什么是这样计算呢?可以利用乘法分配律进行验证,两两相乘后可以得到N
的所有约数(因为这些质因数本就是由约数拆解而来,这步相当于还原了用质数还原了约数),进行加和就是答案
约定 n n n个整数 a i a_i ai,请你输出这些数的乘积的约数之和,答案对 1 0 9 + 7 10^9+7 109+7取模
#include <iostream>
#include <algorithm>
#include <unordered_map>
using namespace std;
typedef long long LL;
const int mod = 1e9 + 7;
int n;
int main() {
cin >> n;
unordered_map<int,int> primes;
// 质因数分解
while (n -- ) {
int x;
cin >> x;
for (int i = 2; i <= x/i; i++)
while (x % i == 0) {
x /= i;
primes[i]++;
}
if (x > 1) primes[x]++;
}
LL res = 1;
// 约数之和公式
for(auto prime: primes) {
int p = prime.first, a = prime.second;
LL t = 1;
while(a--) t = (t * p + 1) % mod;
res = (res * t) % mod;
}
cout << res << endl;
return 0;
}
最大公约数
欧几里得算法,即辗转相除法,是将两个数辗转相除直到余数为0
O ( l o g ( n ) ) O(log(n)) O(log(n))
int gcd(int a, int b) {
return b? gcd(b, a % b):a;
}
欧拉函数
在数论,对正整数n,欧拉函数是小于等于n的正整数中与n互质的数的数目
互质:若N个整数的最大公因数是1,则称这N个整数互质。
若要求出N
的欧拉函数,首先将N
质因数分解,p
代表质数,
α
\alpha
α为p
的次数
N = p 1 α 1 + p 2 α 2 + p 3 α 3 + . . . + p n α n N=p_1^{\alpha_1}+p_2^{\alpha_2}+p_3^{\alpha_3}+...+p_n^{\alpha_n} N=p1α1+p2α2+p3α3+...+pnαn
那么欧拉函数公式如下:
ϕ ( n ) = N ∗ ( 1 − 1 p 1 ) ∗ ( 1 − 1 p 2 ) ∗ . . . ∗ ( 1 − 1 p k ) \phi(n)=N*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k}) ϕ(n)=N∗(1−p11)∗(1−p21)∗...∗(1−pk1)
证明如下(容斥原理):
N同时代表个数,首先我们减去N中所有质因数倍数的个数
N − N p 1 − N p 2 − . . . − N p k N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k} N−p1N−p2N−...−pkN
但这样减会出现一个问题,有的数同时是p1、p2的倍数,就会被减去两次,所以我们需要加回来
N − N p 1 − N p 2 − . . . − N p k + N p 1 ∗ p 2 + N p 2 ∗ p 3 + . . . N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}+\frac{N}{p_1*p_2}+\frac{N}{p_2*p_3}+... N−p1N−p2N−...−pkN+p1∗p2N+p2∗p3N+...
一些数同时是p1、p2、p3的倍数,那么减去3次的同时,又加上了3次,相当于没有计算,因此
N − N p 1 − N p 2 − . . . − N p k + N p 1 ∗ p 2 + N p 2 ∗ p 3 + . . . − N p 1 ∗ p 2 ∗ p 3 − N p 2 ∗ p 3 ∗ p 4 N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}+\frac{N}{p_1*p_2}+\frac{N}{p_2*p_3}+...-\frac{N}{p_1*p_2*p_3}-\frac{N}{p_2*p_3*p_4} N−p1N−p2N−...−pkN+p1∗p2N+p2∗p3N+...−p1∗p2∗p3N−p2∗p3∗p4N
根据以上式子,可以看出一定的规律,将它们进行整理即可得到欧拉函数
求一个数的欧拉函数
int phi(int x) // 欧拉函数
{
int res = x;
for (int i = 2; i <= x / i; i ++ )
if (x % i == 0)
{
res = res / i * (i - 1);
while (x % i == 0) x /= i;
}
if (x > 1) res = res / x * (x - 1);
return res;
}
筛选法求欧拉
给定一个数n
,求出1-n
中每个数欧拉函数的和。这里使用到了前面线性筛质数的方法
-
若
i
为质数,那么i
的欧拉函数 ϕ ( i ) \phi(i) ϕ(i)为i-1
-
若
i
不为质数,求 ϕ ( p r i m e s [ j ] ∗ i ) \phi(primes[j]*i) ϕ(primes[j]∗i)分两种情况-
若
i % primes[j] == 0
,意味着primes[j]
是i
的一个质因子此时 ϕ ( i ) = i ∗ ( 1 − 1 p 1 ) ∗ ( 1 − 1 p 2 ) ∗ . . . ∗ ( 1 − 1 p k ) \phi(i)=i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k}) ϕ(i)=i∗(1−p11)∗(1−p21)∗...∗(1−pk1)
则 ϕ ( p r i m e s [ j ] ∗ i ) = p r i m e s [ j ] ∗ i ∗ ( 1 − 1 p 1 ) ∗ ( 1 − 1 p 2 ) ∗ . . . ∗ ( 1 − 1 p k ) \phi(primes[j]*i)=primes[j]*i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k}) ϕ(primes[j]∗i)=primes[j]∗i∗(1−p11)∗(1−p21)∗...∗(1−pk1)
对比两个式子发现,后者只比前者多乘了一个 p r i m e s [ j ] primes[j] primes[j],因此我们得到公式为
ϕ ( p r i m e s [ j ] ∗ i ) = p r i m e s [ j ] ∗ ϕ ( i ) \phi(primes[j]*i)=primes[j]*\phi(i) ϕ(primes[j]∗i)=primes[j]∗ϕ(i)
-
若
i % primes[j] != 0
,意味着意味着primes[j]
不是i
的一个质因子此时 ϕ ( i ) = i ∗ ( 1 − 1 p 1 ) ∗ ( 1 − 1 p 2 ) ∗ . . . ∗ ( 1 − 1 p k ) \phi(i)=i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k}) ϕ(i)=i∗(1−p11)∗(1−p21)∗...∗(1−pk1)
则 ϕ ( p r i m e s [ j ] ∗ i ) = p r i m e s [ j ] ∗ i ∗ ( 1 − 1 p 1 ) ∗ ( 1 − 1 p 2 ) ∗ . . . ∗ ( 1 − 1 p k ) ∗ ( 1 − 1 p j ) \phi(primes[j]*i)=primes[j]*i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})*(1-\frac{1}{p_j}) ϕ(primes[j]∗i)=primes[j]∗i∗(1−p11)∗(1−p21)∗...∗(1−pk1)∗(1−pj1)
对比两个式子,因为此时的
primes[j]
不是i
的一个质因子,所以在求primes[j]*i
时,需要比原先式子多乘了primes[j]
和(1-1/pj)
,化简得最后结果如下ϕ ( p r i m e s [ j ] ∗ i ) = p [ j ] ∗ ϕ ( i ) ∗ ( 1 − 1 p j ) = ϕ ( i ) ∗ ( p j − 1 ) \phi(primes[j]*i)=p[j]*\phi(i)*(1-\frac{1}{p_j})=\phi(i)*(p_j-1) ϕ(primes[j]∗i)=p[j]∗ϕ(i)∗(1−pj1)=ϕ(i)∗(pj−1)
-
int get_eulers(int n) {
phi[1] = 1;
for (int i = 2; i <= n; i++) {
if (!st[i]) {
primes[cnt++] = i;
phi[i] = i - 1;
}
for (int j = 0; primes[j] <= n / i; j++) {
st[primes[j]*i] = true;
if (i % primes[j] == 0) {
phi[i*primes[j]] = primes[j] * phi[i];
break;
}
phi[i*primes[j]] = (primes[j]-1) * phi[i];
}
}
LL res = 0;
for (int i = 1; i <= n; i++) res += phi[i];
return res;
}
欧拉数论定理
内容
若
a
a
a与
n
n
n互质(gcd(a,n)=1
),则$a^{\phi(n)} \equiv 1 ; mod ; n $
其中 ϕ ( n ) \phi(n) ϕ(n)称为对模 n n n缩系的元素个数
- 当n为质数时,$a^{n-1} \equiv 1 ; mod ; n $(费马小定理)
证明
1~n
中,与n
互质的数为
a
1
,
a
2
,
a
3
,
.
.
.
,
a
ϕ
(
n
)
a_1,a_2,a_3,...,a_{\phi(n)}
a1,a2,a3,...,aϕ(n)(注意这里的a与前面的a不同含义)
将上述的数同时乘以a
,得
a a 1 , a a 2 , a a 3 , . . . , a a ϕ ( n ) aa_1,aa_2,aa_3,...,aa_{\phi(n)} aa1,aa2,aa3,...,aaϕ(n)
因a1,a2,a3
与n
互质并且a
与n
也互质,所以两个数的乘积也和n
互质。将两组分别内部同时相乘,得
a 1 ∗ a 2 ∗ a 3 ∗ . . . ∗ a ϕ ( n ) ≡ a a 1 ∗ a a 2 ∗ a a 3 ∗ . . . ∗ a a ϕ ( n ) a_1*a_2*a_3*...*a_{\phi(n)} \equiv aa_1*aa_2*aa_3*...*aa_{\phi(n)} a1∗a2∗a3∗...∗aϕ(n)≡aa1∗aa2∗aa3∗...∗aaϕ(n)
a 1 ∗ a 2 ∗ a 3 ∗ . . . ∗ a ϕ ( n ) ≡ a ϕ ( n ) ( a 1 ∗ a 2 ∗ a 3 ∗ . . . ∗ a ϕ ( n ) ) a_1*a_2*a_3*...*a_{\phi(n)} \equiv a^{\phi(n)}(a_1*a_2*a_3*...*a_{\phi(n)}) a1∗a2∗a3∗...∗aϕ(n)≡aϕ(n)(a1∗a2∗a3∗...∗aϕ(n))
约去相同部分,得
a ϕ ( n ) ≡ 1 ( m o d n ) a^{\phi(n)} \equiv 1 (mod \;n) aϕ(n)≡1(modn)
因此我们得出结论
$a^{\phi(n)} \equiv 1 ; mod ; n $
应用
这个定理可以用来简化幂的模运算。比如计算
7
222
7^{222}
7222的个位数,实际是求
7
222
7^{222}
7222被10
除的余数。7
和10
互质,且φ(10)=4。由欧拉定理知
7
4
≡
1
(
m
o
d
10
)
7^4\equiv 1\ (mod\ 10)
74≡1 (mod 10)。所以
7
222
=
(
7
4
)
55
∗
7
2
≡
49
≡
9
(
m
o
d
10
)
7^{222}=(7^{4})^{55}*7^2\equiv 49\equiv 9\ (mod 10)
7222=(74)55∗72≡49≡9 (mod10)
快速幂
1 ≤ a , k , p ≤ 1 0 9 1 \le a,k,p \le 10^9 1≤a,k,p≤109
以 l o g ( k ) log(k) log(k)的时间复杂度,快速求出 a k m o d p a^k \ mod\ p ak mod p的结果
- 要快速求出上述的结果,我们需要先预处理出 a 2 0 , a 2 1 . . . ( m o d p ) a^{2^0},a^{2^1}...(mod\ p) a20,a21...(mod p)的值
- 然后将 a k a^k ak进行二进制分解, a k = a 2 0 ∗ a 2 1 ∗ . . . = a 2 0 + 2 1 + . . . ( m o d p ) a^k=a^{2^0}*a^{2^1}*...=a^{2^0+2^1+...}(mod\ p) ak=a20∗a21∗...=a20+21+...(mod p)
- 最后用原先预处理的结果,相乘取模得出答案
此算法主要的时间花在了预处理
O
(
l
o
g
(
k
)
)
O(log(k))
O(log(k))和k
的填二进制拆解
O
(
l
o
g
(
k
)
)
O(log(k))
O(log(k))上
虽然上述过程有一些复杂,但在最后的算法实现非常简洁
// a^k%p
int qmi(int a, int k, int p) {
int res = 1;
// 同时处理了两个过程
while(k) {
if (k & 1) res = (LL)res * a % p;
k >>= 1;
a = (LL)a * a % p;
}
return res;
}
扩展欧几里得算法
求出一组 x i , y i x_i,y_i xi,yi,使其满足 a i ∗ x i + b i ∗ y i = g c d ( a i , b i ) a_i*x_i+b_i*y_i=gcd(a_i,b_i) ai∗xi+bi∗yi=gcd(ai,bi)
推导如下:
b y + ( a m o d b ) x = d by+(a\ mod\ b)x=d by+(a mod b)x=d
⇒ b y + ( a − [ a b ] ∗ b ) x = d \Rightarrow by+(a-[\frac{a}{b}]*b)x=d ⇒by+(a−[ba]∗b)x=d
⇒ a x + b ( y − [ a b ] ∗ x ) = d \Rightarrow ax+b(y-[\frac{a}{b}]*x)=d ⇒ax+b(y−[ba]∗x)=d
int exgcd(int a, int b, int& x, int& y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a%b, y, x);
y -= a / b * x;
return d;
}
线性同余方程
题目描述
给定 n n n组数据 a i , b i , m i a_i,b_i,m_i ai,bi,mi,对于每组数求出一个 x i x_i xi,使其满足 a i ∗ x i ≡ b i ( m o d m i ) a_i*x_i \equiv b_i(mod\ m_i) ai∗xi≡bi(mod mi)
推导
a i ∗ x i ≡ b i ( m o d m i ) a_i*x_i \equiv b_i(mod\ m_i) ai∗xi≡bi(mod mi)
⇒ a i ∗ x i / m i = y i . . . b i \Rightarrow a_i*x_i / m_i = y_i...b_i ⇒ai∗xi/mi=yi...bi
⇒ a i ∗ x i − m i ∗ y i = b i \Rightarrow a_i*x_i - m_i*y_i=b_i ⇒ai∗xi−mi∗yi=bi
令 y i ′ = y i y_i\prime = y_i yi′=yi
⇒ a ∗ x + m ∗ y ′ = b \Rightarrow a*x + m*y\prime=b ⇒a∗x+m∗y′=b
如果b
不是最大公约数d
的倍数,那么一定无解
根据扩展欧几里得算法,对于任意的整数a
和m
,必定存在一组x
和y'
,使得ax+my'
等于gcd(a,m)
所以可以求得 a ∗ x + m ∗ y ′ = d a*x + m* y\prime =d a∗x+m∗y′=d
⇒ a ∗ x ∗ b d + m ∗ y ′ ∗ b d = b \Rightarrow a*x*\frac{b}{d} + m*y\prime*\frac{b}{d}=b ⇒a∗x∗db+m∗y′∗db=b
所以要求解的xi
=x*b/d
中国剩余定理
中国剩余定理(Chinese remainder theorem, CRT)
小学数学时大家都学过这样的问题:有一群人,3个人一组,剩1人,5个人一组,剩3人,10人一组,剩7人。
那么有多少人呢?CRT的出现就是用来解决上述问题
以上的问题可以看成一个线性同余方程组
CRT公式: x = ( ∑ i = 1 n a i ∗ M i − 1 ∗ M i ) m o d M x=(\sum_{i=1}^n a_i*M_i^{-1}*M_i)\ mod \ M x=(∑i=1nai∗Mi−1∗Mi) mod M
以上问题可以转换成以下的线性同余方程组表达
{ x ≡ a 1 ( m o d m 1 ) x ≡ a 2 ( m o d m 2 ) ⋯ ⋯ x ≡ a k ( m o d m k ) \begin{cases} x\equiv a_1(mod \ m_1)\\ x\equiv a_2(mod \ m_2) \\ \cdots \cdots \\ x\equiv a_k(mod \ m_k)\end{cases} ⎩ ⎨ ⎧x≡a1(mod m1)x≡a2(mod m2)⋯⋯x≡ak(mod mk)
M = m 1 ∗ m 2 ∗ ⋯ ∗ m k M=m_1*m_2*\cdots*m_k M=m1∗m2∗⋯∗mk
M i = M / m i = ⋯ ∗ m i − 1 ∗ m i + 1 ∗ ⋯ M_i=M/m_i=\cdots*m_{i-1}*m_{i+1}*\cdots Mi=M/mi=⋯∗mi−1∗mi+1∗⋯
此外,还需要求出 M i M_i Mi的逆元 M i − 1 M_i^{-1} Mi−1,使其满足条件 M i ∗ M i − 1 ≡ 1 ( m o d m i ) M_i*M_i^{-1} \equiv 1(mod \ m_i) Mi∗Mi−1≡1(mod mi)
而对于本组中的其它i
,满足条件
M
j
∗
M
j
−
1
≡
0
(
m
o
d
m
i
)
M_j*M_j^{-1} \equiv 0(mod \ m_i)
Mj∗Mj−1≡0(mod mi)
根据 x = a i ∗ M i − 1 ∗ M i x=a_i*M_i^{-1}*M_i x=ai∗Mi−1∗Mi,求出其中的一个值
对每一线性同余方程执行上述过程,最后将结果对M
取模,得出最后答案
举个例子,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
即,一个整数除以3余2,除以5余3,除以7余2,求这个整数。
M = 3 ∗ 5 ∗ 7 = 105 M=3*5*7=105 M=3∗5∗7=105
首先求解 x ≡ 2 ( m o d 3 ) x \equiv 2 (mod \ 3) x≡2(mod 3)
M 1 = M / m 1 = 105 / 3 = 35 M_1=M/m_1=105/3=35 M1=M/m1=105/3=35
M 1 − 1 ∗ 35 ≡ 1 ( m o d 3 ) M_1^{-1}*35 \equiv 1 (mod \ 3) M1−1∗35≡1(mod 3),此步利用扩展欧几里得算法,求出 M 1 − 1 = 2 M_1^{-1} = 2 M1−1=2
x = a i ∗ M i − 1 ∗ M i = 2 ∗ 2 ∗ 35 = 2 ∗ 70 = 140 x=a_i*M_i^{-1}*M_i=2*2*35=2*70=140 x=ai∗Mi−1∗Mi=2∗2∗35=2∗70=140
同理,可以求出 x ≡ 3 ( m o d 5 ) x \equiv 3 (mod \ 5) x≡3(mod 5)和 x ≡ 2 ( m o d 7 ) x \equiv 2 (mod \ 7) x≡2(mod 7)
x = 15 ∗ 2 + 21 ∗ 3 + 70 ∗ 2 = 233 m o d M = 233 m o d 105 = 23 x=15*2+21*3+70*2=233 \ mod \ M = 233 \ mod \ 105 = 23 x=15∗2+21∗3+70∗2=233 mod M=233 mod 105=23
这个余数23就是合乎条件的最小数.
扩展CRT
根据题目所述,这里的 a 1 , a 2 , ⋯ , a n a_1,a_2,\cdots,a_n a1,a2,⋯,an不满足两两互质的条件,所以需要重新推导CRT
{ x ≡ m 1 ( m o d a 1 ) x ≡ m 2 ( m o d a 2 ) \begin{cases} x\equiv m_1(mod \ a_1)\\ x\equiv m_2(mod \ a_2) \\ \end{cases} {x≡m1(mod a1)x≡m2(mod a2)
化简得,
{ x = k 1 a 1 + m 1 x = k 2 a 2 + m 2 \begin{cases} x= k_1a_1+m_1\\ x= k_2a_2+m_2 \\ \end{cases} {x=k1a1+m1x=k2a2+m2
所以 k 1 a 1 + m 1 = k 2 a 2 + m 2 k_1a_1+m_1=k_2a_2+m_2 k1a1+m1=k2a2+m2
⇒ k 1 a 1 − k 2 a 2 = m 2 − m 1 \Rightarrow k_1a_1-k_2a_2=m_2-m_1 ⇒k1a1−k2a2=m2−m1,这个式子有解则满足 g c d ( a 1 , a 2 ) ∣ m 2 − m 1 gcd(a_1,a_2) | m_2-m_1 gcd(a1,a2)∣m2−m1
运用扩展欧几里得算法,能够求出k1,k2的通解,非齐次=齐次通解+非齐次特解
{ k 1 = k 1 + k a 2 d k 2 = k 2 + k a 1 d \begin{cases} k_1= k_1+k\frac{a_2}{d}\\ k_2= k_2+k\frac{a_1}{d} \\ \end{cases} {k1=k1+kda2k2=k2+kda1
x = k 1 a 1 + m 1 x=k_1a_1+m_1 x=k1a1+m1
⇒ x = ( k 1 + k a 2 d ) a 1 + m 1 \Rightarrow x=(k_1+k\frac{a_2}{d})a_1+m_1 ⇒x=(k1+kda2)a1+m1
⇒ x = a 1 k 1 + m 1 + k a 1 a 2 d \Rightarrow x=a_1k_1+m_1+k\frac{a_1a_2}{d} ⇒x=a1k1+m1+kda1a2
其中a1a2/d
就是a1a2
的最小公倍数,因此
⇒ x = a 1 k 1 + m 1 + k [ a 1 a 2 ] \Rightarrow x=a_1k_1+m_1+k[a_1a_2] ⇒x=a1k1+m1+k[a1a2]
通过对比最前面化简后的式子,我们可以把
a
1
k
1
+
m
1
a_1k_1+m_1
a1k1+m1看成m
,
k
[
a
1
a
2
]
k[a_1a_2]
k[a1a2]看成ka
⇒ x = x 0 + k a \Rightarrow x=x_0+ka ⇒x=x0+ka
以这样的形式重复更新合并后来的线性同余方程,可以得到最终答案
#include <iostream>
#include <cmath>
using namespace std;
typedef long long LL;
int n;
LL exgcd(LL a,LL b, LL& x, LL& y) {
if (!b) {
x = 1, y = 0;
return a;
}
LL d = exgcd(b, a%b, y, x);
y -= a / b * x;
return d;
}
int main() {
cin >> n;
LL a1, m1;
cin >> a1 >> m1;
bool has_answer = true;
for (int i = 0; i < n-1; i ++ ) {
LL a2, m2;
cin >> a2 >> m2;
LL k1, k2;
LL d = exgcd(a1, a2, k1, k2);
if ((m2-m1) % d) { // 当m2-m1不是d的倍数时,无解
has_answer = false;
break;
}
k1 *= (m2 - m1)/d; // 欧几里得求出的解是以最小公约数d为结果,转换为m2-m1
LL t = a2 / d;
k1 = (k1 % t + t) % t; // 最小正整数解
m1 = a1*k1+m1; // 为下一轮循环计算m1,a1
a1 = abs(a1 / d * a2); // 不能a1*a2就近相乘,否则会溢出
}
if (has_answer) cout << (m1 % a1 + a1) % a1 << endl;
else puts("-1");
return 0;
}
高斯消元
- 枚举每一列
- 找到该列中绝对值最大的那个数所在的行
- 将该行移动到最上面(次上、次次上)
- 将该行的第一个数变为1
- 用这个数把下面行该列的数消成0