组合数学专题

专题简介

本专题包含了一些组合数学中常见的套路和方法,如拉格朗日插值,动态规划,容斥原理,狄利克雷卷积,线性筛,杜教筛 等等.

目录

1. 2018 四川省赛GRISAIA (数论分块)
2. HDU 6428 Calculate (狄利克雷卷积,线性筛)
3. BZOJ4559 成绩比较 (动态规划,拉格朗日插值)
4. BZOJ 2633 已经没有什么好害怕的了 (容斥森林,动态规划)
5. BZOJ 3930 选数 (容斥原理,动态规划)
6. 徐州网络赛 D Easy Math (积性函数,线性筛)
7. NWERC2015 Debugging (动态规划,数论分块)

---

四川省赛 GRISAIA

问题描述

求${\sum }_{i=1}^n{\sum }_{j=1}^{i}n\%(i\ast j)$,其中$1\le n\le 10^{11}$

引入的一些性质

$\lceil \frac{n}{xy}\rceil =\lceil \frac{\lceil \frac{n}{x}\rceil }{y}\rceil$
$\lfloor \frac{n}{xy}\rfloor =\lfloor \frac{\lfloor \frac{n}{x}\rfloor }{y}\rfloor$

解法

先对原问题进行变形:

${\sum }_{i=1}^n{\sum }_{j=1}^{i}n\%(i\ast j)$

$={\sum }_{i=1}^n{\sum }_{j=1}^i(n-\lfloor (\frac{n}{ij})\rfloor \ast ij)$

只需要求:

${\sum }_{i=1}^n{\sum }_{j=1}^i\lfloor (\frac{n}{ij})\rfloor \ast ij$

$=\frac{1}{2}{\sum }_{i=1}^n{\sum }_{j=1}^n(\lfloor (\frac{n}{ij})\rfloor \ast ij)+\frac{1}{2}{\sum }_{i=1}^n\lfloor \frac{n}{i^2}\rfloor \ast i^2$

只需要求:

${\sum }_{i=1}^n{\sum }_{j=1}^n(\lfloor (\frac{n}{ij})\rfloor \ast ij)$

这里我们就要使用引入的性质了:

$={\sum }_{i=1}^n{\sum }_{j=1}^n(i\ast \lfloor \frac{\lfloor \frac{n}{i}\rfloor }{j}\rfloor \ast j)$

$={\sum }_{i=1}^{n}i{\sum }_{j=1}^n(\lfloor \frac{\lfloor \frac{n}{i}\rfloor }{j}\rfloor \ast j)$

记录

$f(n)={\sum }_{i=1}^n\lfloor \frac{n}{i}\rfloor \ast i$

显然

${\sum }_{i=1}^n{\sum }_{j=1}^n(i\ast \lfloor \frac{\lfloor \frac{n}{i}\rfloor }{j}\rfloor \ast j)={\sum }_{i=1}^{n}i\ast f(\lfloor \frac{n}{i}\rfloor )$

显然$f(n)$的计算可以数论分块,然后上面部分的计算也可以分块.

---

HDU6428 [Calculate]

知识清单

1. 欧拉函数
2. 莫比乌斯函数
3. 莫比乌斯反演
4. 狄利克雷卷积
5. 线性筛

题解

要求计算${\sum }_{i=1}^A{\sum }_{j=1}^b{\sum }_{k=1}^C(\varphi (gcd(i,j^2,k^3)))$
我们先对式子进行变换,比较套路一步是改为枚举$gcd$,然后统计$gcd$出现的次数,即:

$={\sum }_{d=1}\varphi (d)({\sum }_{i=1}^A{\sum }_{j=1}^b{\sum }_{k=1}^C[gcd(i,j^2,k^3)=d])$

对上述式子中的${\sum }_{i=1}^A{\sum }_{j=1}^b{\sum }_{k=1}^C[gcd(i,j^2,k^3)=d]$再进行莫比乌斯反演:

设$g(d)={\sum }_{i=1}^A{\sum }_{j=1}^b{\sum }_{k=1}^C[gcd(i,j^2,k^3)=d]$

并设$G[n]={\sum }_{n|d}g(d)$

那么显然有$g(n)={\sum }_{n|d}\mu (\frac{d}{n})G(d)$

那么要求的式子可以写成

$={\sum }_{d=1}\varphi (d)g(d)={\sum }_{d=1}\varphi (d){\sum }_{d|t}\mu (\frac{t}{d})G(t)$

由于$\varphi (x)$ 与$\mu (x)$都是积性函数,故将$G[t]$和$\varphi (d)$的位置做交换,使之形成狄利克雷卷积.

$={\sum }_{t=1}G(t){\sum }_{d|t}\mu (\frac{t}{d})\varphi (d)$

在这里,记$F(x)={\sum }_{d|t}\mu (\frac{t}{d})\varphi (d)$

$F(x)$是积性函数的狄利克雷卷积,因此$F(x)$也是积性函数,可以用线性筛法求得.

线性筛所能用到的信息:

$F(p)=p-2$

$F(p^x)=F(p^{x-1})\ast p$

$F(p^2)=p^2-2\ast p+1$

有关$G(n)$的计算:

若$n=\prod p_i^{a_i}$

则定义

$g_2(n)=\prod p_i^{\lceil \frac{a_i}{2}\rceil }$

$g_3(n)=\prod p_i^{\lceil \frac{a_i}{3}\rceil }$

那么$G(n)=[\frac{A}{n}][\frac{B}{g_2(n)}][\frac{C}{g_3(n)}]$

这样的话式子$={\sum }_{t=1}G(t)F(t)$中的$G(t)$和$F(t)$均可在$O(n)$线性时间复杂度内预处理出来.

另外:这道题卡常数,真的恶心.

代码

#include <bits/stdc++.h>using namespace std;typedef long long ll;
const ll MOD = 1<<30;
const int maxn = 10000000;
int zhi[maxn+10],prime[maxn+10],pcnt,F[maxn+10];
int low[maxn+10],lowm[maxn+10],g2[maxn+10],g3[maxn+10];
void sieve(){zhi[1] = 1;low[1] = F[1] = g3[1] = g2[1] = 1;lowm[1] = 0;for(int i = 2;i <= maxn;++i){if(!zhi[i]) {F[i] = i-2,prime[pcnt++] = i;low[i] = i;lowm[i] = 1;g3[i] = g2[i] = i;}for(int j = 0;j < pcnt && prime[j]*i <= maxn;++j){zhi[i*prime[j]] = 1;if(i % prime[j] == 0){low[i*prime[j]] = low[i] * prime[j];lowm[i*prime[j]] = lowm[i]+1;if(i == low[i]){if(i == prime[j]) F[i*prime[j]] = prime[j]*prime[j]-2*prime[j]+1;else F[i*prime[j]] = F[i] * prime[j];}else{F[i*prime[j]] = F[i/low[i]]*F[low[i]*prime[j]];}g2[i*prime[j]] = g2[i];g3[i*prime[j]] = g3[i];if(lowm[i*prime[j]] % 2 == 1) g2[i*prime[j]] *= prime[j];if(lowm[i*prime[j]] % 3 == 1) g3[i*prime[j]] *= prime[j];break;}else{F[i*prime[j]] = F[i] * F[prime[j]];low[i*prime[j]] = prime[j];lowm[i*prime[j]] = 1;g2[i*prime[j]] = g2[i] * prime[j];g3[i*prime[j]] = g3[i] * prime[j];}}}
}int T ,A,B,C;
int main(){sieve();cin>>T;while(T--){scanf("%d%d%d",&A,&B,&C);ll ans = 0;for(int n = 1;n <= max(A,max(B,C));++n){ans = (ans + (A/n)*(B/g2[n])%MOD*(C/g3[n])%MOD*F[n]%MOD) % MOD;}cout<<ans<<endl;}return 0;
}

---

BZOJ4559 [成绩比较]

知识清单

1. 动态规划
2. 拉格朗日插值

题解

分两部分来解决:

1.先不考虑具体分数,只考虑满足相对排名计算得到方案数.

定义$dp[i][j]$表示仅考虑前$i$节课,B神恰好碾压$j$个人.
转移方法:
显然$dp[i][j]$由$dp[i][j+k]$转移过来.
具体转移方法:
考虑第$i$门课,从$j+k$中选出$k$个,他们的分数大于B神,剩下$j$个放在B神的下面,表示继续被碾压.方案总数$C_{k+j}^k$.
剩下的还有$N-1-(j+k)$名同学,这些同学中取出$Rank[i]-k$名同学放在B神上面,填满$Rank[i]$.方案总数$C_{N-1-j-k}^{Rank[i]-k}$.
因此转移方程:
$dp[i][j]={\sum }_{k=0}^{k+j\le N-1}dp[i][j+k]C_{k+j}^k{C}_{N-1-j-k}^{Rank[i]-k}$
dp初始化:
$dp[0][N-1]=1$

该部分具体实现代码:

// dp[i][j]表示仅考虑前i门课程,其中有j名同学被B神碾压,的相对排名的方案数.  
dp[0][N-1] = 1;
for(int i = 1;i <= M;++i){for(int j = N-1;j >= 0;--j){for(int k = 0;k+j<=N-1;++k){ll res = C[j+k][k] * C[N-1-(j+k)][R[i]-1-k] % MOD;res = res * dp[i-1][j+k] % MOD;dp[i][j] = (dp[i][j] + res) % MOD;}}
}

2.求出满足题目给出的排名条件的分数具体分配方案数.

对于第$i$门课,设B神的排名为$Rank[i]$,则该门课产生的方案数就是枚举该门课B神的分数得到的方案数之和:
$L_{Rank[i]}[U_x]={\sum }_{t=1}^{Ux}{t}^{N-Rank[i]}(U_x-t{)}^{Rank[i]-1}$
这个式子是关于$U_x$的$N$次多项式,而$Ux$最大则可能达到$1e9$,显然不可能直接来求,我们考虑拉格朗日插值.

---

$Lagrange$插值简介

对于一个$n$次的多项式,只需要$n+1$个点即可唯一确定,而用$n+1$个点恢复这个$n$次多项式的方法就是$Lagrange$插值法.

设$(x_i,y_i)$为我们已经获得的$N+1$个点.

定义$f_i[x]={\prod }_{i!=j}\frac{x-x_j}{x_i-x_j}$

那么$L[x]={\sum }_{i=1}^{N+1}{y}_i{f}_i[x]$

---

如果这道题用$Lagrange$插值求得了$L[x]$,那么第$i$门课的答案即是$L_{Rank[i]}[U_x]$.

//拉格朗日实现具体代码
ll Lagrange(int id){//拉格朗日插值for(int xi = 1;xi <= N+1;++xi){x[xi] = xi;y[xi] = 0;for(int i = 1;i <= xi;++i){y[xi] = (y[xi] + (mod_pow(i,N-R[id])*mod_pow(xi-i,R[id]-1)%MOD)) % MOD;}}ll ans = 0;for(int xi = 1;xi <= N+1;++xi){ll res = y[xi];for(int xj = 1;xj <= N+1;++xj){if(xj == xi) continue;res = res * (U[id] - xj + MOD) % MOD;res = res * div(xi - xj + MOD) % MOD;}ans = (ans + res) % MOD;}return ans;
}

综合起来,这道题的答案就是$ans=dp[M][K]\ast {\prod }_{i=1}^M{L}_{Rank[i]}[U_x]$

###实现代码

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <stack>using namespace std;#define pr(x) cout<<#x<<":"<<x<<endl
typedef long long ll;
int N,M,K;
const ll MOD = 1e9+7;
ll mod_pow(ll x,ll n){if(n == 0) return 1;ll res = 1;while(n){if(n&1) res = res * x %MOD;x = x * x % MOD;n >>= 1;}return res;
}
const int maxn = 107;
ll R[maxn],U[maxn],dp[maxn][maxn];
ll C[maxn][maxn];
ll div(ll x){return mod_pow(x,MOD-2);
}
void initC(){C[0][0] = 1;for(int i = 1;i < maxn;++i){C[i][0] = 1;for(int j = 1;j <= i;++j){C[i][j] = (C[i-1][j] + C[i-1][j-1])%MOD;}}
}
ll x[maxn],y[maxn];
ll Lagrange(int id){//拉格朗日插值for(int xi = 1;xi <= N+1;++xi){x[xi] = xi;y[xi] = 0;for(int i = 1;i <= xi;++i){y[xi] = (y[xi] + (mod_pow(i,N-R[id])*mod_pow(xi-i,R[id]-1)%MOD)) % MOD;}}ll ans = 0;for(int xi = 1;xi <= N+1;++xi){ll res = y[xi];for(int xj = 1;xj <= N+1;++xj){if(xj == xi) continue;res = res * (U[id] - xj + MOD) % MOD;res = res * div(xi - xj + MOD) % MOD;}ans = (ans + res) % MOD;}return ans;
}
int main(){initC();ios::sync_with_stdio(false);cin>>N>>M>>K;for(int i = 1;i <= M;++i) cin>>U[i];for(int i = 1;i <= M;++i) cin>>R[i];// dp[i][j]表示仅考虑前i门课程,其中有j名同学被B神碾压,的相对排名的方案数.dp[0][N-1] = 1;for(int i = 1;i <= M;++i){for(int j = N-1;j >= 0;--j){for(int k = 0;k+j<=N-1;++k){ll res = C[j+k][k] * C[N-1-(j+k)][R[i]-1-k] % MOD;res = res * dp[i-1][j+k] % MOD;dp[i][j] = (dp[i][j] + res) % MOD;}}}ll ans = dp[M][K];for(int i = 1;i <= M;++i){ans = ans * Lagrange(i) % MOD;}cout<<ans<<endl;return 0;
}

---

BZOJ2633 [已经没有什么好害怕的了]

知识清单

1. 容斥原理
2. 动态规划

题意

给出$n$个数 $a_i$和$n$个数$b_i$,将$a_i$和$b_i$进行匹配,使得匹配中$a_i>b_i$的对数比$a_i<b_i$的匹配数大$k$组.

题解

首先将$a,b$数组排序,然后对于每个$i$求出满足$b_j<a_i(1\le j\le n)$的$j$的个数,并记为$cnt[i]$.

定义$f[i][j]$表示仅考虑前$i$个$a$中,满足有$j$对匹配使得$a_u>b_v$成立的方案数.(显然$1\le j\le i$) (注意这个定义中不要求剩下的$a$和$b$进行匹配.)

那么我们可以得到转移方程:
$f[i][j]=f[i-1][j]+f[i-1][j-1]\ast (cnt[i]-(j-1))$

看起来似乎$f[n][j]\ast (n-j)!$就是我们要找的答案?

但显然并不是,因为$(n-j)!$表示的剩下未匹配的元素进行随机搭配的时候,有时候也会搭配出满足$a_u>b_v$的匹配,像这种就是不合法的,因此我们还要进行去重操作.

设$g[i]$表示在$n$个匹配中,恰好有$i$个匹配满足$a_u>b_v$的方案数.显然这回$g[(n+k)/2]$就是我们要的答案.

容易想到转移方程:
$g[i]=f[n][i]\ast (n-i)!-g[i+1]\ast C_{i+1}^i-g[i+2]\ast C_{I+2}^i-...-g[n]\ast C_n^i$

解释一下为什么要乘以$C_j^i$这个系数:因为$f[n][i]\ast (n-i)!$产生的满足有$j$个$a_u>b_v$的匹配数的方案数可以由固定$i$个匹配数来得到(剩下的$j-i$个匹配由$(j-i)!$部分负责).而固定$i$个匹配数的方案数刚好就是$C_j^i$

得到这个式子倒着转移就可以了.

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <stack>using namespace std;
const ll MOD = 1e9+9;
int n,k;
const int maxn = 2007;
ll f[maxn][maxn],g[maxn];
ll a[maxn],b[maxn];
ll fac[maxn];
ll C[maxn][maxn];
ll cnt[maxn];
int main(){C[0][0] = fac[0] = 1;for(int i = 1;i <= 2000;++i) {fac[i] = fac[i-1] * i % MOD;C[i][0] = 1;}for(int i = 1;i <= 2000;++i) {for(int j = 1;j <= i;++j){C[i][j] = (C[i-1][j] + C[i-1][j-1]) % MOD;}}cin>>n>>k;if((n+k)%2 != 0) return puts("0");for(int i = 1;i <= n;++i) scanf("%lld",&a[i]);for(int i = 1;i <= n;++i) scanf("%lld",&b[i]);sort(a+1,a+1+n);sort(b+1,b+1+n);for(int i = 1;i <= n;++i){cnt[i] = lower_bound(b+1,b+1+n,a[i])-b-1;}f[0][0] = 1;for(int i = 1;i <= n;++i){f[i][0] = 1;for(int j = 1;j <= i;++j){f[i][j] = f[i-1][j] + (f[i-1][j-1]*max(cnt[i]-j+1,0LL)%MOD);f[i][j] %= MOD;}}for(int i = n;i >= (n+k)/2;--i){g[i] = f[n][i]*fac[n-i]%MOD;for(int j = i+1;j <= n;++j){g[i] = (g[i] + MOD - (C[j][i] * g[j] % MOD)) % MOD;}}cout<<g[(n+k)/2]<<endl;return 0;
}

---

BZOJ3930 [选数]

知识清单

1. 容斥原理 (或莫比乌斯反演)

题解

题目给的一个很关键的信息就是$H-L\le 1e5$.这意味着如果在区间$[L,H]$中选出$N$个不完全相同的整数,那么他们的$GCD\le H-L\le 1e5$
证明:$GCD(x,y)=GCD(x-y,y)<=x-y$

如果将$[L,H]$中的各个元素除以$K$,就相当于在$(\frac{L-1}{K},\frac{H}{K}]$中找$N$个不完全相同的数使得其$GCD=1$

设$F[x]$表示满足在$(\frac{L-1}{K},\frac{H}{K}]$中找$N$个不完全相同的元素,使得其$GCD=1$的方案数.

显然$F[x]=([\frac{H}{Kx}]-[\frac{L-1}{Kx}]{)}^N-([\frac{H}{Kx}]-[\frac{L-1}{Kx}])-F[2\ast x]-F[3\ast x]-....$

遇见这种递推式子直接倒着进行转移就可以了.
注意,如果 $K$ 存在于区间$[L,H]$之间,那么$ans++$,因为这种属于$N$个数字全部相等的情况,我们的算法是没有考虑进去的.

另外,
这题实际上可以使用莫比乌斯反演,而且莫比乌斯反演形式特别直接,但是复杂度不对.

重新定义$F[x]$为在$[L,H]$中找$N$个数使得其$GCD=X$的方案数
我们可以得到
${\sum }_{n|d}F[d]=([\frac{H}{n}]-[\frac{L-1}{n}]{)}^N$
反演得到,
$F[n]={\sum }_{n|d}\mu (\frac{d}{n})([\frac{H}{d}]-[\frac{L-1}{d}]{)}^N$

我看到其他的博客有讲到分块的方法,但感觉复杂度有点玄学.

const ll MOD = 1e9+7;
ll N,K,L,R;
ll mod_pow(ll x,ll n){ll res = 1;while(n){if(n & 1) res = res * x % MOD;x = x * x % MOD;n >>= 1;}return res;
}
ll F[200007];
const int maxn = 200000;
int main(){cin>>N>>K>>L>>R;ll ans = K >= L && K <= R;R /= K;L --;L /= K;ll len = R - L;for(ll d = len;d >= 1;--d){ll res = R/d - L/d;F[d] = (mod_pow(res,N)-res+MOD) % MOD;for(ll k = 2*d;k <= len;k += d)F[d] = (F[d] - F[k] + MOD) % MOD;}cout<<ans+F[1]<<endl;return 0;
}

徐州网络赛 Easy Math

前置技能

1. 容斥原理
2. 线性筛
3. 杜教筛

题解

${\sum }_{i=1}^m\mu (in)$

当$n$存在平方因子的时候,直接输出$0$作为答案,这是很显然的.

否则,$n=p_1{p}_2...p_k$,不含平方因子.

考虑区间$[1,m]$中的数字$i$,如果$i$不含有因子$p_1$,则可以使用积性函数的性质.

即:$\mu (i\ast n)=\mu (i\ast \frac{n}{p_1}\ast p_1)=-\mu (i\ast \frac{n}{p_1})$

而区间$[1,m]$中的数并不全都与$p_1$互质,这就会造成$\mu (i\ast n)=0$ 而 $\mu (i\ast \frac{n}{p_1})̸=0$的情况,这就要求我们用容斥的思想来去掉这部分.

如果$i$中恰好有一个因子$p_1$,则$\mu (i\ast \frac{n}{p_1})=\mu (\frac{i}{p_1}\ast n)$

因此要减去${\sum }_{i=1}^{\lfloor \frac{n}{p_1}\rfloor }\mu (in)$

而如果$i$中有多于$1$个的$p_i$因子,那么这部分的莫比乌斯函数值无论是在那部分中都是$0$,因此不会对最终答案造成影响.

这样相当于降低了问题的规模.

我们可以列出递推式子,然后递推解决这个问题.

${\sum }_{i=1}^n\mu (in)=-{\sum }_{i=1}^n\mu (i\ast \frac{n}{p_1})+{\sum }_{i=1}^{\lfloor \frac{n}{p_1}\rfloor }\mu (in)$

递推终止的条件是

1. $n==1$ 时候,使用杜教筛得到莫比乌斯函数的前缀和
2. $\sum$的上限为$0$时候,直接返回$0$.

typedef long long ll;
const int N = 10000000;;
int mu[N+10],prime[N+10],pcnt,zhi[N+10];
void sieve(){pcnt = 0;mu[1] = zhi[1] = 1;for(int i = 2;i <= N;++i){if(!zhi[i]) mu[i] = -1,prime[pcnt++] = i;for(int j = 0;j < pcnt && prime[j]*i <= N;++j){zhi[i*prime[j]] = 1;if(i % prime[j] == 0){mu[i*prime[j]] = 0;break;}else{mu[i*prime[j]] = -mu[i];}}}for(int i = 1;i <= N;++i) mu[i] += mu[i-1];
}
map<ll,int> rec,vis;
int Mu(int x){if(x <= N) return mu[x];if(vis[x]) return rec[x];int res = 1,now = x,nxt;while(now >= 2){nxt = x/(x/now+1);res -= (now - nxt) * Mu(x/now);now = nxt;}vis[x] = 1;return rec[x] = res;
}
ll m,n;
ll np[20];int npcnt;
int Ans(ll u,ll d){if(d == 1) return Mu(u);if(u == 0) return 0;for(int i = 0;i < npcnt;++i){if(d % np[i] == 0){return -Ans(u,d/np[i]) + Ans(u/np[i],d);}}
}
bool divide(){ll x = n;for(ll i = 2;i * i <= n;++i){int cc = 0;if(x % i == 0) np[npcnt++] = i;while(x % i == 0) x /= i,cc++;if(cc >= 2) return false;}if(x > 1) np[npcnt++] = x;return true;
}
int main(){sieve();cin>>m>>n;if(!divide()) return puts("0");cout<<Ans(m,n)<<endl;return 0;
}

---

NWERC2015 Debugging

题目描述

有一份包含一个Bug的$n,(1\le n\le 10^6)$行代码,运行一次到崩溃的时间为$r,1\le r\le 10^9$.
现在你可以在任意一行花费时间$p,1\le p\le 10^9$设置一个$printf$语句来判断程序是否运行到这里.
请问在最坏情况下,最少需要多少时间可以定位到bug所在行.

题解

因为我们调试的时候往往会运行很多次程序,而每次运行完之后,我们都能定位bug在哪一块代码行中.

因此我们动态规划的思想得到递推公式.
记$dp[n]$表示调试$n$行代码所需要花的时间.

$dp[n]=mi{n}_{1\le k\le n}(dp[\frac{n}{k+1}]+k\ast p+r)$

直接递推的时间复杂度是$O(n^2)$的,而如果我们发现对于$[\frac{n}{k+1}]$相等的$k$,我们只取$k$的最小值就可以了.

所以通过底数分块,我们可以求得转移的时间复杂度为$\sqrt{n}$,而真正有效的状态数也只有$\sqrt{n}$,所以时间复杂度上限不会超过$O(n)$.