[Acwing] 算法基础课总结 四 数学知识

目录

质数1.试除法判断质数2.分解质因数3.筛质数 约数1. 试除法求约数2.约数个数3. 约数之和4.最大公约数 欧拉函数1.欧拉函数2.筛法求欧拉函数 快速幂1.快速幂2.快速幂求逆元 扩展欧几里得1.扩展欧几里得算法2.线性同余方程 中国剩余定理高斯消元1. 高斯消元解线性方程组2.高斯解异或线性方程组 组合数1.求组合数 I 2.求组合数 II3.求组合数 III4. 求 组合数IV5.满足条件的01序列 容斥原理1.能被整除的数 前言题意思路code

质数

1.试除法判断质数

这个是最朴素的判断质数的方法, O ( n ) O(\sqrt n) O(n ​)的因为大于 n \sqrt n n ​的不可能为约数

bool is_prime(int x){if (x < 2) return false;for (int i = 2; i <= x / i; i ++ )if (x % i == 0)return false;return true;}

2.分解质因数

根据算术基本定理,我们可以将一个数拆分成多个质数的次幂的乘积,因此我们只需要遇到约数,让后将其除尽即可

void divide(int x){for (int i = 2; i <= x / i; i ++ )if (x % i == 0){int s = 0;///计算质因数的指数while (x % i == 0)x /= i, s ++ ;cout << i << ' ' << s << endl;}///2的情况if (x > 1)cout << x << ' ' << 1 << endl;cout << endl;}

3.筛质数

普通筛法 :

O ( n l o g n ) O(nlogn) O(nlogn)

先塞除质数,然后利用质数的倍数将所有倍数打上标记,因此没有标记到的就是质数了。

void get_primes2(){for(int i=2;i<=n;i++){if(!st[i]) primes[cnt++]=i;//把素数存起来for(int j=i;j<=n;j+=i){//不管是合数还是质数，都用来筛掉后面它的倍数st[j]=true;}}}

埃氏筛法 :

O ( n l o g l o g n ) O(nloglogn) O(nloglogn)

我们没必要在没有判断出质数的时候就进行倍数去重了,我们可以在遇到质数的时候在这样干

void get_primes1(){for(int i=2;i<=n;i++){if(!st[i]){primes[cnt++]=i;for(int j=i;j<=n;j+=i) st[j]=true;//可以用质数就把所有的合数都筛掉；}}}

线性筛法 :

O ( n ) O(n) O(n)

建议直接背过

void get_primes(){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;}}}

约数

1. 试除法求约数

首先约数总是成对出现的, 因此我们在判断 x % i x\%i x%i的时候,可以一起判断 x / i x/i x/i

因此可以在 O ( n + M l o g M ) O(\sqrt n+MlogM) O(n ​+MlogM)

vector<int> get_divisors(int x){vector<int> res;for (int i = 1; i <= x / i; i ++ )if (x % i == 0){res.push_back(i);if (i != x / i) res.push_back(x / i);}sort(res.begin(), res.end());return res;}

2.约数个数

unordered_map<int,int> primes;//映射函数while(n--){int x;scanf("%d",&x);for(int i=2;i<=x/i;i++)while(x%i==0){primes[i]++;x/=i;//方便求得约数的数量}if(x>1) primes[x]++;//x的最大公约数可能大于sqrt(x);}long long res=1;for(auto p:primes) res=res*(p.second+1)%mod;//将统计出来的数按照由图中公式所得出来的结论得出答案printf("%lld\n",res);

3. 约数之和

const int N = 110, mod = 1e9 + 7;int main(){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;primes[i] ++ ;}if (x > 1) primes[x] ++ ;}LL res = 1;for (auto p : primes){LL a = p.first, b = p.second;LL t = 1;while (b -- ) t = (t * a + 1) % mod;res = res * t % mod;}cout << res << endl;return 0;}

4.最大公约数

int gcd(int a, int b){return b ? gcd(b,a%b):a; }

欧拉函数

欧拉函数,就是 1 − N 1-N 1−N中与 N N N互质的个数就是欧拉函数

因此我们是通过筛出N N N的质因子之后,进行操作欧拉函数的

1.欧拉函数

int phi(int x){int res = x;for(int i=2; i<=x/i; i++)if(x%i == 0 ){/// N * (1-1/p)的变形res = res/i * (i-1);while(x%i == 0 ) x/=i;}if(x>1) res = res /x *(x-1);return res;}

2.筛法求欧拉函数

因为欧拉函数与质因子息息相关,因此我们筛选质数的同时可以处理欧拉函数

const int N = 1000010;int primes[N], cnt;int euler[N];bool st[N];void get_eulers(int n){euler[1] = 1;for (int i = 2; i <= n; i ++ ){if (!st[i]){primes[cnt ++ ] = i;euler[i] = i - 1;}for (int j = 0; primes[j] <= n / i; j ++ ){int t = primes[j] * i;st[t] = true;if (i % primes[j] == 0){euler[t] = euler[i] * primes[j];break;}euler[t] = euler[i] * (primes[j] - 1);}}}

快速幂

1.快速幂

我们可以将一个数的次幂进行二进制化, 然后我们对二进制进行操作即可

Code :

int qmi(int x){int res = 1;while(x){if(x&1) res = res*a;a = a*a;x >>=1;}return res;}

2.快速幂求逆元

当 p p p是质数的时候,由费马小定理可以得知

x [ b ] = b p − 2 % p x[b] = b^{p-2}\%p x[b]=bp−2%p

qmi(a, p - 2, p)

扩展欧几里得

1.扩展欧几里得算法

最重要的一个就是 g c d ( a , b ) = g c d ( b , a % b ) gcd(a,b)= gcd(b,a\%b) gcd(a,b)=gcd(b,a%b)

因此我们可以通过递归,求出下一层的 x ′ y ′ x'y' x′y′然后回代

int exgcd(int a, int b, int &x, int &y){//返回gcd(a,b) 并求出解(引用带回)if(b==0){x = 1, y = 0;return a;}int x1,y1,gcd;gcd = exgcd(b, a%b, x1, y1);x = y1, y = x1 - a/b*y1;return gcd; }

2.线性同余方程

a ∗ x + y ∗ m = b \begin{aligned} &a*x+y*m=b\\ \end{aligned} ​a∗x+y∗m=b​

显然当 g c d ( a , m ) ∣ b gcd(a,m)|b gcd(a,m)∣b的时候有解

因此我们不妨先用扩展欧几里得求解出

a ∗ x 0 + y 0 ∗ m = g c d ( a , m ) \begin{aligned} &a*x_0+y_0*m=gcd(a,m)\\ \end{aligned} ​a∗x0​+y0​∗m=gcd(a,m)​

然后我们将式子左右两边同时乘上 b g c d ( a , m ) \frac{b}{gcd(a,m)} gcd(a,m)b​

( a ∗ x 0 + y 0 ∗ m ) ∗ b g c d ( a , m ) = g c d ( a , m ) ∗ b g c d ( a , m ) a ∗ x 0 ∗ b g c d ( a , m ) + y 0 ∗ b g c d ( a , m ) = b \begin{aligned} &(a*x_0+y_0*m)*\frac{b}{gcd(a,m)}=gcd(a,m)*\frac{b}{gcd(a,m)}\\ &a*x_0*\frac{b}{gcd(a,m)}+y_0*\frac{b}{gcd(a,m)}=b \end{aligned} ​(a∗x0​+y0​∗m)∗gcd(a,m)b​=gcd(a,m)∗gcd(a,m)b​a∗x0​∗gcd(a,m)b​+y0​∗gcd(a,m)b​=b​

因此最后对于上面由扩展欧几里得求出的 x 0 x_0 x0​可以用来推导出所求式的 x x x

x = x 0 ∗ g c d ( a , m ) / b % m x = x_0*gcd(a,m)/b\%m x=x0​∗gcd(a,m)/b%m

#include <iostream>#include <algorithm>using namespace std;typedef long long LL;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;scanf("%d", &n);while (n -- ){int a, b, m;scanf("%d%d%d", &a, &b, &m);int x, y;int d = exgcd(a, m, x, y);if (b % d) puts("impossible");else printf("%d\n", (LL)b / d * x % m);}return 0;}

中国剩余定理

证明同扩展欧几里得,求的是一个最小的x满足

x ≡ m i m o d ( a i ) x\equiv m_i\ mod(a_i) x≡mi​mod(ai​)

code :

#include<iostream>using namespace std;typedef long long LL;//数据范围比较大,所以用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;LL a1,m1;cin>>n>>a1>>m1;LL x=0;for(int i=1;i<n;i++){LL a2,m2;cin>>a2>>m2;LL k1,k2;LL d=exgcd(a1,a2,k1,k2);if((m2-m1)%d){x=-1;break;}k1*=(m2-m1)/d;//因为此时k1是k1*a1+k2*a2=d的解,所以要乘上(m2-m1)/d的倍数大小LL t=abs(a2/d);k1=(k1%t+t)%t;//数据比较极端,所以只求k的最小正整数解m1=k1*a1+m1;//m1在被赋值之后的值为当前"x"的值，此时赋值是为了方便下一轮的继续使用a1=abs(a1*a2/d);//循环结束时a1的值为当前所有的a1,a2,……an中的最小公倍数}if(x!=-1)x=(m1%a1+a1)%a1;//当循环结束时，此时的值应该与最小公倍数取模,以求得最小正整数解printf("%lld\n",x);return 0;}作者：灰之魔女链接：/solution/content/23099/来源：AcWing著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

高斯消元

目前还不清到底有什么当作板子吧

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

typedef priority_queue<int,vector<int>,greater<int>> Pri_m;typedef pair<int,int> pii;typedef vector<int> VI;map<int,int> mp;const int N = 110;const double eps = 1e-6;double a[N][N];int n;int gauss(){int c,r;//c代表列,r代表行for(c = 0 , r= 0 ; 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;//如果当前列的最大数是0 那么没必要再算了 因为他的约束方程可能是上面几个for(int i = c;i<n+1;i++){swap(a[t][i],a[r][i]);//把当前这一行,换到最上面 //换到第r行去}for(int i=n;i>=c;i -- )a[r][i]/=a[r][c];//把这一行 第一个数的系数变为1for(int i = r+1;i<n;i++) //把当前列下所有数都变为1if(fabs(a[i][c]) > eps)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)return 2;//说明 无解return 1;//否则,0 = 0 那么也就是无穷解因为有多个C}//唯一解//从下往上回代,得到方程的解for(int i = n-1;i>=0;i -- )for(int j = i+1;j<n;j++){a[i][n] -=a[j][n] * a[i][j];}return 0;}void solve(){cin>>n;for(int i = 0;i<n;i++){for(int j = 0;j<n+1;j++){cin>>a[i][j];}}int t = gauss();if(t == 0){for(int i = 0 ;i<n;i++){if(fabs(a[i][n]) < eps){cout<<"0.00"<<endl;continue;}printf("%.2lf\n",a[i][n]);}}else if(t == 1){cout<<"Infinite group solutions"<<endl;}else cout<<"No solution"<<endl;}

2.高斯解异或线性方程组

#include <iostream>#include <algorithm>using namespace std;const int N = 110;int n;int a[N][N];int gauss(){int c, r;for (c = 0, r = 0; c < n; c ++ ){int t = r;for (int i = r; i < n; i ++ )if (a[i][c])t = i;if (!a[t][c]) continue;for (int i = c; i <= n; i ++ ) swap(a[r][i], a[t][i]);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 ++ )if (a[i][n])return 2;return 1;}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(){cin >> n;for (int i = 0; i < n; i ++ )for (int j = 0; j < n + 1; 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;}

组合数

1.求组合数 I

这个是一个可以开二维数组,然后使用 d p dp dp预处理的组合数

因为 C a b = C a − 1 b − 1 + C a − 1 b C_a^b=C_{a-1}^{b-1}+C_{a-1}^b Cab​=Ca−1b−1​+Ca−1b​

状态表示 :

从 a a a中选 b b b个的方案

状态计算 :

选第 b b b个,那么只需要从 a − 1 a-1 a−1选 b − 1 b-1 b−1个即可不选第 b b b个,那么需要在 a − 1 a-1 a−1中选 b b b个

for(int i=0;i<=2000;i++){for(int j=0;j<=i;j++){if(!j) f[i][j]=1;else f[i][j]=(f[i-1][j-1]+f[i-1][j])%mod;}}

2.求组合数 II

这里就不能开二维数组了,但是可以将二维数组转换成一维数组

我们可以通过预处理阶乘的逆元将其转换为一维

#include <iostream>#include <algorithm>using namespace std;typedef long long LL;const int N = 100010, mod = 1e9 + 7;int fact[N], infact[N];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 main(){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] * qmi(i, mod - 2, mod) % mod;}int n;scanf("%d", &n);while (n -- ){int a, b;scanf("%d%d", &a, &b);printf("%d\n", (LL)fact[a] * infact[b] % mod * infact[a - b] % mod);}return 0;}

3.求组合数 III

这个很难了,这个是使用卢卡斯定理,优化的 p p p是质数,数据范围很大的情况下用到的,目前还不熟悉记作板子

#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){if (b > a) return 0;int res = 1;for (int i = 1, j = a; i <= b; i ++, j -- ){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);return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p;}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;}

4. 求 组合数IV

这是一个需要高精度的板子,因为不对 p p p取余了

#include <iostream>#include <algorithm>#include <vector>using namespace std;const int N = 5010;int primes[N], cnt;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] <= n / i; j ++ ){st[primes[j] * i] = true;if (i % primes[j] == 0) break;}}}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;}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);}vector<int> res;res.push_back(1);for (int i = 0; i < cnt; i ++ )for (int j = 0; j < sum[i]; j ++ )res = mul(res, primes[i]);for (int i = res.size() - 1; i >= 0; i -- ) printf("%d", res[i]);puts("");return 0;}

5.满足条件的01序列

容斥原理

1.能被整除的数

前言

传送门 :

题意

给定 n , m n,m n,m,同时给定 p 1 , p 2 . . . . p m p_1,p_2....p_m p1​,p2​....pm​

询问 1 − n 1-n 1−n中能被 p 1 , p 2 . . . p m p_1,p_2...p_m p1​,p2​...pm​中至少一个整除的数有多少个

思路

我们设 S i S_i Si​为 1 − n 1-n 1−n中能被 p i p_i pi​整除的集合

⋃ i = 1 m S i = S 1 + S 2 + . . . . S m − ( S 1 ⋂ S 2 + S 1 ⋂ S 3 . . . . S m − 1 ⋂ S m ) + ( S 1 ⋂ S 2 ⋂ S 3 + . . . . S m − 2 ⋂ S m − 1 ⋂ S m ) + . . . + ( − 1 ) m − 1 ( ⋂ i = 1 m S ) \begin{aligned} \bigcup\limits_{i=1}^mS_i=S_1+S_2+....S_m-(S_1\bigcap S_2+S_1\bigcap S_3....S_{m-1}\bigcap S_m)+(S_1\bigcap S_2 \bigcap S_3+....S_{m-2}\bigcap S_{m-1}\bigcap S_m)+...+(-1)^{m-1}(\bigcap_{i=1}^mS) \end{aligned} i=1⋃m​Si​=S1​+S2​+....Sm​−(S1​⋂S2​+S1​⋂S3​....Sm−1​⋂Sm​)+(S1​⋂S2​⋂S3​+....Sm−2​⋂Sm−1​⋂Sm​)+...+(−1)m−1(i=1⋂m​S)​

(为什么会有上式,下面例举出来了一个较为简易的图解)

但是通过上式子,我们需要知道的结论是

如果选中的集合个数是奇数,那么前面系数就是 + + +否则就是 − - −

当然现在还有两个疑问就是

如何求集合的大小如何求集合的 ⋂ \bigcap ⋂

1.如何求集合的大小

对于一个 p i p_i pi​,在 1 − n 1-n 1−n中有多少个能被整除的数呢 ？

这个显然我们可以通过枚举倍数 1 ∗ p i , 2 ∗ p i , 3 ∗ p i . . . . k ∗ p i > n 1*p_i,2*p_i,3*p_i....k*p_i>n 1∗pi​,2∗pi​,3∗pi​....k∗pi​>n

因此显然我们可以得知 k = n / p i k=n/p_i k=n/pi​ (先说一下这种推导方式很常用)

扩展 :

例如,如何求 1 − n 1-n 1−n中能被开 x x x次方的数有多少个

我们设 a x > n a^x>n ax>n则 a < = n ∗ p o w ( 1 , 1.0 / x ) a<=n*pow(1,1.0/x) a<=n∗pow(1,1.0/x)

如何求集合的 ⋂ \bigcap ⋂

因为 p i p_i pi​都是质数,因此他们的交际一定是他们的最小公倍数的所有倍数

即 S 1 ⋂ S 2 = n / ( p 1 ∗ p 2 ) S_1\bigcap S_2=n/(p_1*p_2) S1​⋂S2​=n/(p1​∗p2​)

好了现在所有计算都已经会弄了,我们考虑如何枚举|遍历

因为 m m m很小,所以我们考虑使用二进制枚举进行遍历,然后我们判断每一个位置即可

code

#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，集合数量+1t *= 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;}