复杂数论

学习
字数统计: 8.2k   |   阅读时长≈ 39 分钟

复杂数论

最没有自信的一章哈哈😅

1. 质数判断

时间复杂度 $O(n)→O(\sqrt n)$

1
2
3
4
5
6
7
8
9
// 1:朴素试除
bool isPrime(int n) {
	if(n<2)
		return false;
	for(int i=2;i<=n;i++) {
		if(n%i==0) 
			return false;
	return true;
}
1
2
3
4
5
6
7
8
9
10
// 2:枚举次数优化→sqrt:这个函数本身速度较慢,不推荐
// 包括 i*i<=n 这种写法也不推荐,因为i*i可能小于等于n,但是(i+1)*(i+1)可能会溢出int型最大值,如果溢出了就会变成一个负数,从而影响最终答案的判断
bool isPrime(int n) {
	if(n<2)
		return false;
	for(int i=2;i<=sqrt(n);i++) {
		if(n%i==0) 
			return false;
	return true;
}
1
2
3
4
5
6
7
8
9
// 3:最佳写法应该是i<=n/i,速度快并且不会溢出
bool isPrime(int n) {
	if(n<2)
		return false;
	for(int i=2;i<=n/i;i++) {	
		if(n%i==0) 
			return false;
	return true;
}

2. 分解质因数

质因数的底数和指数：底数指质数的基数，比如分解12=2×2×3时，2和3就是底数，而指数指的是各个基数出现的次数，2的指数是2，3的指数是1

2.1. 试除法

时间复杂度 $O(n)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void divede(int n) {
	for(int i=2;i<=n;i++) {
		// 先找到这个质数
		if(n%i==0) {
			// 再求质数出现的次数
			int s=0;
			// 当n%i不等于0了,说明不再被整除,指数求完了
			while(n%i==0) {
				n/=i;
				s++;
			}
			cout<<i<<' '<<s<<endl;
		}
	}
}

2.2. 枚举次数优化

  • 对上述算法进行优化，我们可以证明n中最多包含一个大于sqrt(n)的质因子（假设有两个的话，两个大于sqrt(n)的质因子相乘的结果就大于n了，所以不成立），所以可以和优化求质数时一样的做法来减少枚举次数。

时间复杂度 $[O(log(n)),\ O(sqrt(n))]$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void divede(int n) {
	for(int i=2;i<=n/i;i++) {
		if(n%i==0) {
			int s=0;
			while(n%i==0) {
				n/=i;
				s++;
			}
			cout<<i<<' '<<s<<endl;
		}
	}
	// 因为我们证明了大于sqrt(n)的数最多只有一个
	// 那么n如果还>1,说明没除完,那剩下的那个必定就是大于sqrt(n)的质因子
	if(n>1)
		cout<<n<<' '<<1<<endl;
}

3. 质数筛

3.1. 朴素法

时间复杂度约等于 $O(nln^n)<O(nlog^n)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
const int N=1e6+10; // 筛选1~1e6+10中的质数
// primes:记录质数小于等于N的所有质数,cnt:记录质数个数
// primes中下标为0~cnt-1存储的就是所有小于等于N的质数
int primes[N],cnt; 
bool st[N]; // 记录每个数的访问状态,为true的说明不是质数

void get_primes(int n) {
	// 从2开始遍历数字,并把其对应倍数的数字全部设置成已访问
	for(int i=2;i<=n;i++) {
		// 如果没有访问过
		if(!st[i]) {
			primes[cnt++]=i;
		}
		// 从i*2开始,只要是i的倍数都设置成已访问
		for(int j=i+i;j<=n;j+=i)
			st[j]=true;
	}
}

3.2. 埃及筛

枚举优化，时间复杂度 $O(nlg^{lg^n})$，约等于 $O(n)$

  • 质数定理：1~n当中有 $n/ln(n)$ 个质数，我们在遍历的时候只把质数选出来
1
2
3
4
5
6
7
8
9
10
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;
		}
	}
}

3.3. 欧拉筛

时间复杂度 $O(n)$，真正意义上的线性筛

  • 线性筛法就在欧拉筛的基础上改进了一个东西，$n$ 只会被最小质因子筛掉，欧拉筛、埃及筛都是建立在2是公认的最小的质数上做的，建议疯狂看注释理解嗷
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
const int N=1e6+10;
bool st[N];
int primes[N],cnt;

void get_primes(int n) {
	for(int i=2;i<=n;i++) {
		if(!st[i])
			primes[cnt++]=i;
		// primes[j](质数)*i<=n,遍历primes[j]即为:primes数组中当前的所有质数
		for(int j=0;primes[j]<=n/i;j++) {
			// primes[j]一定是i的最小质因子,primes[j]一定也是primes[j]*i的最小质因子
			st[primes[j]*i]=true;
			// 对于一个合数x,假设primes[j]是x的最小质因子,当i枚举到x/primes[j]时就break
			// 啥子意思,就是说我i是个合数的话就只筛一遍,这样才能达到近似O(n)的效果
			// 比如i=6时,质数集中有2,3,5这三个数,那么我只筛2*6=12,3*6和5*6不筛
			// 为什么不筛?因为我在算3*6=18时,i=9,9*2也为18,其实也筛掉了
			// 所以对于合数,我们只用其最小质因子来筛,提前break就可以避免重复筛
			if(i%primes[j]==0)
				break;
		}
	}
}

4. 约数

4.1. 试除法求约数个数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// 求取约数的函数
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);
			// 如果i和n/i是同一个数字,那么只用放一遍
			if(i!=n/i)
				res.push_back(n/i);
		}
	}
	sort(res.begin(),res.end()); // 从小到大排序
	return res;
}

4.2. 公式求约数个数

  • 基于算术基本定理，$N$ 的约数的个数就和 $[β_1,β_k]$ 的取法是一模一样的；所以 $[β_1,β_k]$ 就是 $α_1$ 的选法，共 $(α_1+1)$ 种，也表示 $α2$ 的选法，共 $(α_2+1)$ 种，根据乘法定理，所以：对于一个数 $N$，约数个数一定是 $(α_1+1)×(α_2+1)×(α_3+1)×···(α_k+1)$
  • 只要我们能求解出一个数分解质因数的结果，就能求出其约数个数

image-20231201152205955

  • PS：int范围内，约数个数最多的数的约数个数大概是1500左右。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

// 把因式分解每个基数的底数和指数存到一个map中
const int mod=1e9+7;

int main() {
	// 求n个数乘积的约数的个数,并对mod取余
	int n;
	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;
				// 对这n个数的约数都用primes来存,最后统一算出约数个数
				primes[i]++; 
			}
		}
		// 如果x还大于1,把剩下大的那个质因数给加上
		if(x>1)
			primes[x]++;
	}
	
	// primes中存储了所有的底数和其对应的指数
	ll res=1;
	for(auto prime:primes)
		res=res*(prime.second+1)%mod;
	cout<<res<<endl;
	return 0;
}

4.3. 约数之和

  • 只要我们能求解出一个数分解质因数的结果，就能求得其约数之和，把n的所有约数展开，公式就是下图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

// 把因式分解每个基数的底数和指数存到一个map中
const int mod=1e9+7;

int main() {
	// 求n个数乘积的约数的个数,并对mod取余
	int n;
	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;
				// 对这n个数的约数都用primes来存,最后统一算出约数个数
				primes[i]++; 
			}
		}
		// 如果x还大于1,把剩下大的那个质因数给加上
		if(x>1)
			primes[x]++;
	}
	
	// primes中存储了所有的底数和其对应的指数
	ll res=1;
	for(auto prime:primes) {
		int p=prime.first; // 底数
		int a=prime.second; // 指数
		ll t=1;
        // 每一次提一个p出来
        // 1+p1(p1^0+p1^1+···+p1^(α1-1))
		while(a--)
			t=(t*p+1)%mod;
		res=res*t%mod; // 结果同样对mod取余
	}
	cout<<res;
	return 0;
}

4. 最大公约数

  • 欧几里得算法，也叫辗转相除法

  • 如果a能除d，b能除d，那么(a+b)也能整除d，比如6能被2整除，14能被2整除，那么20能被2整除，甚者还有 a·x+b·y 能被d整除

  • 根据辗转相除法，a和b的最大公约数就等于b和a mod b的最大公约数。怎么推的？

  • $a%b=a-floor(a/b)×b$ ，写成 $a-c×b$，那么 $(a,b)=(b, a - c × b)$

  • 对于左侧，$(a,b)$，假设 $d$ 能整除 $a$ ，$d$ 也能整除 $b$，对于右侧，$d$ 能整除 $b$，现在只需要证明 $d$ 也能整除 $a$ ，就能说明左右两侧等价了，$d$ 能整除 $a - c × b$，只需要在右侧加上一个 $c × b$ (相当于 $a×x+b×y$ )，就能说明 $d$ 能整除 $a$ 了，所以左右两侧等价。证毕，$(a,b)=(b,a%b)$

1
2
3
4
int gcd(int a,int b) {
    // 当b=0时,(a,0)的最大公约数一定是a,所以返回值的第二项是a
	return b?gcd(b,a%b):a;
}

5. 欧拉函数

时间复杂度$O(sqrt(n))$，来源于分解质因数的$O(sqrt(n))$

5.1. 定义法

  • $φ(n)$表示1~n中与n互质的个数，比如$φ(6)$，$[1,6]$ 中有1,5与6互质（互质：公约数只有1的两个整数，叫做互质整数），所以$φ(6)$=2

image-20231201162530871

  • 根据容斥原理可以推导，容斥原理：

  1. 从 $[1,N]$ 中去掉 $p_1, p_2, … , p_k$ 的所有倍数，而 $[1,N]$ 当中 $p_k$ 的倍数等于$N-N/p_1-N/p_2-···-N/p_k$，此时有些数可能既是$p_1$的倍数，也是$p_2$的倍数，就多去除了一次，所以此时进行第二步；

  2. 加上所有$p_i×p_j$的倍数，i和j是1~K中任意两个数，得到下式：

    $N-N/p_1-N/p_2-···-N/p_k+N/p_1×p_2+N/p_1×p_3+···N/p_{k-1})×p_k$

  3. 这时候又会遇到一个问题，如果一个数既是$p_1$、$p_2$，也是$p_3$的倍呢？它就会被p1减一次，p2减一次，p3再减一次，再在第二步中加上三次，所以第三步；这个时候相当于加回来了，而我们是想把它删去，所以还要减去所有pi×pj×pk的倍数，得到下式：$N-N/p1-N/p2-···-N/pk+N/p1×p2+N/p1×p3+···N/p(k-1)×pk-N/p1×p2×p3-N/p1×p2×p4-···-N/p(k-2)×p(k-1)×pk$，以此类推，最后得到的公式其实就是下面这个函数的展开，所以欧拉函数其实就是容斥原理的展开：

image-20231201163755176

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void get_ola(int n)
{
	int res=n;
   	// 找n的所有质因子
	for(int i=2;i<=n/i;i ++) {
		if(n%i==0) {
			// res=res*(1-1/a):这是欧拉公式
			// 像下面这个公式一样转换是为了避免出现小数,因为这是一个整数运算,保证只有整数
			res=res/i*(i-1);
			while(n%i==0)
				n/=i;
		} 
	}
	if(n>1)
		res=res/n*(n-1);
	cout<<res<<endl;
}

5.2. 筛法

  • 在某些情况下需要求出 $[1,N]$ 中每一个数的欧拉函数，这种情况下用公式就非常的慢，可以借助之前线性筛的思想，用 $O(n)$ 的时间复杂度求出每个数的欧拉函数。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

const int N=1e6+10;
int primes[N];
bool st[N];
int cnt;
int phi[N]; // phi[i]:1~i中有多个数和i互斥
ll res;

// 在线性筛的过程中求欧拉函数
ll get_eulers(int n) {
	phi[1]=1; // 1~1中和1互斥的只有自己
	for(int i=2;i<=n;i++) {
		if(!st[i]) {
			primes[cnt++]=i;
			// 如果i是质数,那么前i-1个数都和i互质
			phi[i]=i-1;
		}
		for(int j=0;primes[j]<=n/i;j++) {
			st[primes[j]*i]=true;
			if(i%primes[j]==0) {
				// i mod pj==0时说明pj是i的质因子
				// phi(pj*i)就等于pj*phi(i)
				phi[primes[j]*i]=primes[j]*phi[i];
				break;
			}
			// i mod pj!=0, pj 一定是 i*pj 的质因子
			// phi(pj*i)就等于p(j-1)*phi(i)
			phi[primes[j]*i]=phi[i]*(primes[j]-1);
		}
	}
	ll res=0;
	for(int i=1;i<=n;i++)
		res+=phi[i];
	return res;
}

int main() {
	// 求1~n中每个数的欧拉函数之和
	int n;
	cin>>n;
	cout<<get_eulers(n)<<endl;
	return 0;
}

5.3. 欧拉定理

  • 如果 $a$ 与 $n$ 互质，那么 $a^phi(n)\ mod\ n$ 同余 $1$，例如：$a=5$, $n=6$ 时有：$5^phi(6) = 5^2 =25，25 mod 6 =1

​ 欧拉定理推论：若 a^phi(p) mod p =1，假设p是质数，则有 a^(p-1) mod p=1，这也是费马(小)定理。如果n=5，那我们可以把1和6视为同样的值，因为1%5 == 6%5 == 1。

6. 快速幂

能在 $O(log^k)$ 的时间下快速求出 $a^k\ mod\ p$ 的结果，其中 $a,\ p,\ k<=10^9$

  • 对于 $a^k\ mod\ p$，快速幂就是预计算出 $a^{(2^0)}\ mod\ p$ 的结果、 $a^{(2^1)}\ mod\ p$ 的结果、… 、 $a^{2^{(log^k)}}\ mod\ p$ 的结果，一共 $log^k$ 个数，如果要把 $a^k$ 组合出来，其实就是一种二进制优化
  • 比如我们要计算 $4^5\ mod\ 10$ ，就预算出 $4^{2^0}\ mod\ 10$、 $4^{2^1}\ mod\ 10$、 $4^{2^2}\ mod\ 10$ 的结果，$5$ 的二进制又等于 $(101)_2$，所以就是 $4^ {2^0} mod 10 × 4^{2^2}\ mod\ 10$

6.1. 快速幂

1
2
3
4
5
6
7
8
9
10
11
12
int qmi(int a,int k,int p) {
	int res=1;
	while(k) {
		// 取2进制最后一位,如果为1,则乘以这一位的位权
		if(k&1)
			res=(ll)res*a%p;
		// 去掉最后一位
		k>>=1;
		a=(ll)a*a%p; // 位权提升
	}
	return res;
}

6.2. 快速幂求逆元

$a/b = a×x\ (mod\ m)$，找到这个 $x$ ，就可以把除法变为乘法，等式两端约去 $a$，得到 $b×x=1(mod\ m)$，即只要找到数字 $x$ ，使得 $b×x\ mod\ m$ 的值等于 $1$，$x$ 就是逆元。

image-20231205080313203

  • 如果 $p$ 是质数，那么这个式子会变成一个费马定理，所以 $b^{p-1}=1(mod p)$，把 $b^(p-1)$ 写成 $b×b^(p-2)$ ， 那么 $b^(p-2)$ 就是 $b\ mod\ p$ 的逆元，又因为最小的质数是 $2$ ，所以 $p-2$ 一定是 $≥0$ 的，所以成立。所以题目变成一个快速幂，求取 $b^{p-2}\ mod\ p$
  • 注意：如果 $p$ 和 $b$ 之间存在倍数关系（不互质）的话，那么 $b×x\ mod\ p$ 一定等于 $0$ ，此时逆元无解；当没有倍数关系时，一定是有解的。只有当保证 $p$ 是质数时，才能使用快速幂求逆元！
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

int qmi(int a,int k,int p) {
	int res=1;
	while(k) {
		// 取2进制最后一位,如果为1,则乘以这一位的位权
		if(k&1)
			res=(ll)res*a%p;
		// 去掉最后一位
		k>>=1;
		a=(ll)a*a%p; // 位权提升
	}
	return res;
}

int main() {
	int n;
	cin>>n;
	while(n--) {
		int a,p;
		cin>>a>>p;
		int res=qmi(a,p-2,p);
		// 如果a是p的倍数,则没有逆元,因为 a*x mod p一定等于0
		if(a%p)
			cout<<res<<endl;
		else
			cout<<"impossible"<<endl;
	}
	return 0;
}

7. 扩展欧几里得算法

7.1. 裴蜀定理

  • 对于任意一对正整数 $a,b$ ，那么一定存在非零整数 $x,y$ ，使得 $a×x+b×y=gcd(a,b)$ ，并且 $a$ 和 $b$ 的最大公约数是 $a$ 和 $b$ 凑出来的最小的正整数
  • 如何证明？$gcd(a,b)$ 是 $a$ 和 $b$ 的最大公约数，那么 $a$ 和 $b$ 的线性组合一定是最大公约数的倍数，那么最小的倍数就是 $1$ 倍呗，所以一定能凑出：$a×x+b×y=gcd(a,b)$

7.2. 扩展欧几里得

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

// 扩展欧几里得算法
int exgcd(int a,int b,int &x,int &y) {
	// 如果b=0,那么gcd(a,b)==a,线性组合即是(1,0)
	if(!b) {
		x=1,y=0;
		return a;
	}

	int d=exgcd(b,a%b,y,x);
	// ①b*y+(a mod b)*x=gcd(a,b)
	// ②又 a mod b=a-floor(a/b)*b→证明:假设r是a mod b的余数,现在我们要求这个r,a/b=floor(a/b)+r,∴r=a-floor(a/b)*b	
	// 二式相结合,有: a*x+b*(y-floor(a/b)*x)=gcd(a,b),即x的系数不变,但是y要减去那一坨
	y-=a/b*x;
	return d;
}

int main() {
	int n;
	cin>>n;
	while(n--) {
		int a,b,x,y;
		cin>>a>>b;
		exgcd(a,b,x,y);
		cout<<x<<' '<<y<<endl; // 注意:x和y不唯一
	}
	return 0;
}

7.3. 【应用】线性同余方程

给定 $a,b,m$，求整数 $x$ 使得 $a×x\ mod\ m=b$（如果 $b$ 为 $1$ 就是求逆元了是吧，如果 $m$ 再为质数，就可以用快速幂了），注意 $x$ 的解不唯一。

$$
ax\equiv b(mod\ m) \ 等价于存在整数y∈Z,\ st. \ ax=m*y+b \ ,即m的若干倍再加上b
\
化简得\ ax-my=b, \ 把-y令成y’，即有\ ax+my’=b, \ 就变成了扩展欧几里得的方程
\
所以根据扩展欧几里得算法可以求出一组系数x,y’使得ax+my’=b, \ d是最大公约数
$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

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;
}

int main() {
	int n;
	cin>>n;
	while(n--) {
		// a*x mod m =b
		int a,b,m;
		cin>>a>>b>>m;
		int x,y;
		int d=exgcd(a,m,x,y);
		// 如果b不是d的倍数方程无解,因为a*x无法通过%m同余等于b
		if(b%d)
			cout<<"impossible"<<endl;
		else
			// 注意这个解不唯一
			// 当b是d的倍数时,方程两边同时初一d,得到 a/d*x=b/d(mod m/d),由于
			// a/d和m/d是互质的,可以找到a/d的逆元x',使得(a/b)*x'=1(mod m/d)
			// 所以让 b/d乘以这个米远x'得到x*(b/d),这就是x的一个特解
			cout<<(ll)x*(b/d)%m<<endl;
	}
	return 0;
}

7.4. 【应用】中国剩余定理

image-20231205095348192

  • 注意：$m_1,\ ···\ ,m_n$ 之间是两两互质的
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

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()
{
	int n;
	cin>>n;
	LL x=0,m1,a1;
	cin>>a1>>m1;
	for(int i=0;i<n-1;i++){
		LL m2,a2;
		cin>>a2>>m2;
		LL k1,k2;
		LL d=exgcd(a1,a2,k1,k2);
		if((m2-m1)%d){
			x=-1;
			break;
		}
		
		//更新状态
		k1*=(m2-m1)/d;
		LL t=a2/d;
		// 将解变成一个最小的正整数解
		k1=(k1%t+t)%t;
		
		x=k1*a1+m1;
		// 更新a和m,k只是个变量,不用管,取余的时候会自动消失
		m1=k1*a1+m1;
		a1=abs(a1/d*a2);
	}
	
	if(x!=-1) x=(m1%a1+a1)%a1;
	
	cout<<x<<endl;
	return 0;
}

8. 高斯消元

8.1. 高斯消元解线性方程组

image-20231205234751191

  • 这种题一般有 $SPJ$ 优化，即 $0.00$ 和 $-0.00$ 都认为是 $0$，即都认为是相同的内容，因为在数学上 $+0$ 和 $-0$ 代表的概念和意义都是一样的。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath> // fabs() = float abs = 浮点数的绝对值，相对应的，整数的绝对值即为abs()
using namespace std;

// 高斯消元的四部操作:
// 1:找到绝对值最大的一行(为了代码的稳定性)
// 2:将该行移到最上面
// 3:将改行第一个数变为1
// 4:将最上面一行下面的所有行的第c列变为0
// 以上4步均可用初等行变换完成

const int N = 110;
const double eps = 1e-8; // 防止出现精度问题

double a[N][N];
int n;

int gauss()
{
	int c = 0, r = 0;
	for(; c < n; c++)
	{
		int t = r;
		for (int i = r; i < n; i ++ ) // 第一步，找绝对值最大的行
		{
			if(fabs(a[i][c]) > fabs(a[t][c]))
				t = i; 
		}
		
		if(fabs(a[t][c]) < eps) continue; // 代表这一列已经被处理过了
		
		for(int i = c; i <= n; i++) // 第二步 
			swap(a[t][i], a[r][i]); // 将这一行调到最上面
		
		for(int i = n; i >= c; i--) // 第三步
			a[r][i] /= a[r][c]; // 要倒着算，否则会影响后面的数
		
		for(int i = r + 1; i < n; i++) // 第四步
		{
			if(fabs(a[i][c]) > eps)  // 如果是0就不用操作了
			{
				for(int j = n; j >= c; j--)
				{
					a[i][j] -= a[r][j] * a[i][c];
				}
			}
		}
		r++;
	}
	
	if(r < n) // 方程数 < n
	{
		for(int i = r; i < n; i++)
		{
			if(fabs(a[i][n]) > eps) // 0 != 0
			{
				return 2; // 无解
			}
		}
		
		return 1; // 无数解, 0=0
	}
	
	for(int i = n - 1; i >= 0; i --) // 有解从下往上回代
	{
		for(int j = i + 1; j < n; j++)
		{
			a[i][n] -= a[i][j] * a[j][n];
		}
	}
	return 0;
	
}


int main()
{
	scanf("%d", &n);
	
	for (int i = 0; i < n; i ++ )
		for (int j = 0; j < n + 1; j ++ )
			scanf("%lf", &a[i][j]);
	
	int t = gauss();
	
	if(t == 0)  // 有唯一解
	{
		// 最后输出的就是每个未知解的值:>
		for (int i = 0; i < n; i ++ )
		{
			if(fabs(a[i][n]) < eps) a[i][n] = 0.00; // 避免输出-0.00
			printf("%.2lf\n", a[i][n]);
			
		}
	}
	
	else if(t == 1) puts("Infinite group solutions"); // 无数解
	else puts("No solution"); // 无解
	
}

8.2. 高斯消元解异或线性方程组

image-20231205235928312
  • 系数和常数的值都为 $0$ 或 $1$ ，每个未知数的取值也为 $0$ 或 $1$ ，方程中间是异或运算：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

// 核心思想:异或→不进位的加法
// 等式与等式之间的异或要一起进行才能保证等式左右两边依然是相等关系
// a^b^c=x,d^f=y,则有a^b^d^c^f=x^y
// 1. 左下角消0
// 1.1 枚举列
// 1.2 找第一个非零行
// 1.3 交换
// 1.4 把同列下面行消零(异或)
// 2. 判断3种情况
// 2.1 唯一解
// 2.2 秩<n
// 2.2.1 有矛盾 无解
// 2.2.2 无矛盾 无穷多解

const int N = 110; // 最多110个系数

int n;
int a[N][N];
int gauss()
{
	int c,r;
	for(c=0,r=0;c<n;c++)
	{
		// 找主元
		int t=-1;
		for(int i=r;i<n;i++)
			if(a[i][c])
			{
				t=i;
				break;
			}
		if(t==-1) continue;
		// 交换主元行
		for(int j=c;j<=n;j++) swap(a[r][j],a[t][j]);
		// 左下角消
		for(int i=r+1;i<n;i++)
			if(a[i][c])//漏了
				for(int j=n;j>=c;j--)//漏了
					a[i][j] ^= a[r][j];
		r++;
	}
	// 判断
	if(r<n)
	{
		for(int i=r;i<n;i++) // i=r
			if(a[i][n])
				return 2;
		return 1;
	}
	// 消右上角
	for(int i=n-1;i>=0;i--)
		for(int j=i+1;j<n;j++)
			//如果是0 就不用下面的a[j][j] 来^a[i][j]了
			//如果不是0 才需要用第j行第j列a[j][j]来^第i行第j列a[i][j] 
			//进而进行整行row[i]^row[j] 间接导致 a[i][n]^a[j][n]
			if(a[i][j])
				a[i][n]^=a[j][n];
	return 0;
}

int main()
{
	cin >> n;
	for(int i=0;i<n;i++)
		for(int j=0;j<=n;j++)
			cin >> a[i][j];
	int t = gauss();
	if(t==0)
	{
		for(int i=0;i<n;i++) cout << a[i][n] << endl;
	}
	else if(t==1) puts("Multiple sets of solutions");
	else puts("No solution");
	return 0;
}

9. 组合数学

9.1. 组合数

9.1.1. 朴素法

时间复杂度 $O(n^2)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
const int mod=1e9+10;
// 初始化C数组
void init() {
	// 从0开始赋值,是为了让C11=C[0][0]+c[0][1]==1,正确的得到1
	// 其实就是把C[N][1~N]全部算出来了而已
	for(int i=0;i<N;i++) {
		for(int j=0;j<=i;j++) {
			if(!j)
				c[i][j]=1;
			else
				c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
		}
	}
}

9.1.2. 快速幂

时间复杂度 $O(a×log^{mod})$

image-20231206090536947
  • 三个阶乘级别的数都可以用快速幂迅速算出，在对组合数进行运算时，我们要对组合数进行模运算，直接用除法是不行的，因为除法得到的结果不一定总为整数（可能是浮点数），此时取模就没有意义，需要找到分母的乘法逆元，然后将除法转换为乘法，快速幂求逆元不就是费马小定理吗，注意模数必须是质数：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include<bits/stdc++.h>
#define x first
#define y second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> PII;

const int mod=1e9+7; // 这是一个质数,所以与2~1e9+6都互质
const int N=1e5+10; // 数字范围
ll fac[N]; // 存储阶乘,基于数组实现记忆化
ll infac[N]; // 存储i!的逆元,infac[i]=(i!)^(-1) mod p

// 快速幂模板
ll qmi(ll a,ll k,ll p) {
	ll res=1;
	while(k) {
		if(k&1==1)
			res=res*a%p;
		k>>=1;
		a=a*a%p;
	}
	return res;
}

int main() {
	int n;
	fac[0]=infac[0]=1; // 0的阶乘以及逆元都是1
	for(int i=1;i<=N;i++) {
		fac[i]=fac[i-1]*i%mod;
		infac[i]=(ll)infac[i-1]*qmi(i,mod-2,mod)%mod; // 根据费马小定理求到i!逆元
	}
	cin>>n;
	while(n--) {
		int a,b;
		cin>>a>>b;
		// cab=a!/(b!)*(a-b)!
		cout<<(ll)fac[a]*infac[b]%mod*infac[a-b]%mod<<endl; // 通过乘法逆元转换为乘法运算
	}
	return 0;
	
}

9.1.3. 卢卡斯定理

时间复杂度 $O(plog_pN)$

  • 对任意非负整数 $a,b$ 和**质数$p$**，有：
  • 将 $m$ 和 $n$ 用 $p$ 进制表示，复杂度为 $O(log_pn)$，如果再用快速幂求逆元计算组合数，复杂度为$O(p)$，故总时间复杂度为$O(plog_pN)$，这里的 $N$ 是 $a,b$ 的时间复杂度约为$10^{18}$，如果预处理阶乘和阶乘逆元，时间复杂度能下降至$O(p+log_pN)$
image-20231206093200137
  • 其证明需要用到两个重要等式，这里不再证明：
image-20231206093244681 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include<iostream>
#include<algorithm>

using namespace std;
typedef long long LL;

int qmi(int a,int k,int p)
{
	int res = 1;
	while(k)
	{
		if(k&1)res = (LL)res*a%p;
		a = (LL)a*a%p;
		k>>=1;
	}
	return res;
}

int C(int a,int b,int p)//自变量类型int
{
	if(b>a)return 0;//漏了边界条件
	int res = 1;
	// a!/(b!(a-b)!) = (a-b+1)*...*a / b! 分子有b项
	for(int i=1,j=a;i<=b;i++,j--)//i<=b而不是<
	{
		res = (LL)res*j%p;
		res = (LL)res*qmi(i,p-2,p)%p;
	}
	return res;
}
//对公式敲
int lucas(LL a,LL b,int p)
{
	if(a<p && b<p)return C(a,b,p);//lucas递归终点是C_{bk}^{ak}
	return (LL)C(a%p,b%p,p)*lucas(a/p,b/p,p)%p;//a%p后肯定是<p的,所以可以用C(),但a/p后不一定<p 所以用lucas继续递归
}

int main()
{
	int n;
	cin >> n;
	while(n--)
	{
		LL a,b;
		int p;
		cin >> a >> b >> p;
		cout << lucas(a,b,p) << endl;
	}
	return 0;
}

9.1.4. 线性筛+高精乘

时间复杂度近似 $O(n)$

  • 题目变得更加简单粗暴，实现方法：由于 $a$ 和 $b$ 的最大取值是 $5000$，所以用线性筛在 $[1,5000]$ 的范围内筛选出质数，求每个质数的次数，再用高精度乘法把所有质因子乘上：
image-20231206144822955 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include<iostream>
#include<algorithm>
#include<vector>

using namespace std;

const int N=5010;

int primes[N],cnt; // primes[N]:依次存放1~N中出现的质数
int sum[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=0;primes[j]*i<=n;j++)
		{
			st[primes[j]*i]=true;
			if(i%primes[j]==0)
				break; // ==0每次漏
		}
	}
}

// 对p的各个<=a的次数算整除下取整倍数
int get(int n,int p)
{
	int res =0;
	while(n)
	{
		res+=n/p;
		n/=p;
	}
	return res;
}

// 高精度乘
vector<int> mul(vector<int> a, int b)
{
	vector<int> c;
	int t = 0;
	for (int i = 0; i < a.size(); i ++ )
	{
		t += a[i] * b;
		c.push_back(t % 10);
		t /= 10;
	}
	while (t)
	{
		c.push_back(t % 10);
		t /= 10;
	}
	// while(C.size()>1 && C.back()==0) C.pop_back();//考虑b==0时才有pop多余的0 b!=0不需要这行
	return c;
}

int main()
{
	int a,b;
	cin >> a >> b;
	get_primes(a);
	
	for(int i=0;i<cnt;i++)
	{
		int p = primes[i];
		sum[i] = get(a,p)-get(a-b,p)-get(b,p); // 是a-b不是b-a
	}
	
	vector<int> res;
	res.push_back(1);
	
	for (int i = 0; i < cnt; i ++ )
		for (int j = 0; j < sum[i]; j ++ ) // primes[i]的次数
			res = mul(res, primes[i]);
	
	for (int i = res.size() - 1; i >= 0; i -- ) 
		printf("%d", res[i]);
	puts("");
	
	return 0;
}

9.2. 卡特兰数

image-20231206150423021

  • 若把 $0$ 和 $1$ 置于坐标轴中，起点定于原点，若 $0$ 表示向右走，$1$ 表示向上走，那么任何前缀中 $0$ 的个数不少于 $1$ 的个数就转换为：对于路径上的任意一点，横坐标大于等于纵坐标，题目所求即为这样的合法路径数量：
Catalan.png
  • 由图可知，任何一条不合法的路径（黑色路径），都对应一条从 $(0,0)$ 走到 $(n-1,n+1)$ 的路径（灰色路径），而任何一条从 $(0,0)$ 走到 $(n-1,n+1)$ 的路径，也对应一条从$(0,0)$ 走到 $(n,n)$ 的不合法路径：$ ans=C_{2n}n-C_{2n}{n-1}=C_{2n}n/n+1$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <bits/stdc++.h>
using namespace std;

typedef long long LL;
const int N = 200010, mod = 1e9 + 7;

int n;
int fact[N], infact[N]; // fact[i]:i!的阶乘,infact[i]:i!阶乘逆元

int ksm(int a, int k) {
	int res = 1;
	while (k) {
		if (k & 1) res = (LL)res * a % mod;
		a = (LL)a * a % mod;
		k >>= 1;
	}
	return res;
}

void init() {
	fact[0] = infact[0] = 1;
	for (int i = 1; i < N; i++) {
		fact[i] = (LL)fact[i - 1] * i % mod;
		infact[i] = (LL)infact[i - 1] * ksm(i, mod - 2) % mod;
	}
}

int main() {
	init();
	cin >> n;
	int res = (LL)fact[2 * n] * infact[n] % mod * infact[n] % mod * ksm(n + 1, mod - 2) % mod;
	cout << res << endl;
	return 0; 
}

10. 容斥原理

  • 注意欧拉函数的推导运用到了容斥原理，欧拉函数求取的是 $[1,n]$ 当中与 $n$ 互质的数的个数：

image-20231206152700721

  • 比如 $n=10,\ p1=2,\ p2=3$，求$[1,10]$中满足能整除 $p_1$ 或 $p_2$ 的个数即：$2,\ 3,\ 4,\ 6,\ 8,\ 9,\ 10$ 共 $7$ 个
image-20231206153027100 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include<iostream>
using namespace std;
typedef long long LL;

const int N = 20;
int p[N], n, m;

int main() {
	cin >> n >> m;
	for(int i = 0; i < m; i++) cin >> p[i];
	int res = 0;
	// 枚举从1 到 1111...(m个1)的每一个集合状态,(至少选中一个集合)
	for(int i = 1; i < 1 << m; i++) {
		int t = 1; // 选中集合对应质数的乘积
		int s = 0; // 选中的集合数量
		// 枚举当前状态的每一位
		for(int j = 0; j < m; j++){
			// 选中一个集合
			if(i >> j & 1) {
				// 乘积大于n,则n/t = 0,跳出这轮循环
				if((LL)t * p[j] > n) {    
					t = -1;
					break;
				}
				s++; // 有一个1,集合数量+1
				t *= p[j];
			}
		}
		if(t == -1) 
			continue;  	
		if(s & 1) 
			res += n / t; // 选中奇数个集合, 则系数应该是1, n/t为当前这种状态的集合数量
		else 
			res -= n / t; // 反之则为 -1
	}
	cout << res << endl;
	return 0;
}
打赏
  • 版权声明： 本博客所有文章除特别声明外，著作权归作者所有。转载请注明出处！
  • Copyrights © 2015-2024 John Doe
  • 访问人数: | 浏览次数:

请我喝杯咖啡吧~

支付宝
微信