狄利克雷卷积
前置
下取整的性质
$\left\lfloor\frac {\lfloor\frac{a}{b}\rfloor}{c}\right\rfloor=\lfloor \frac{a}{bc} \rfloor$。 证明略。
上取整的性质
$\left\lceil\frac {\lceil\frac{a}{b}\rceil}{c}\right\rceil= \lceil \frac{a}{bc} \rceil$
证明：
$\left\lceil\frac {\lceil\frac{a}{b}\rceil}{c}\right\rceil=\left\lfloor\frac {\lfloor\frac{a+b-1}{b}\rfloor+c-1}{c}\right\rfloor=\left\lfloor\frac {\lfloor\frac{a+b-1+bc-b}{b}\rfloor}{c}\right\rfloor=\left\lfloor\frac {a+bc-1}{bc}\right\rfloor= \lceil \frac{a}{bc} \rceil$
$g=\mu*f$，已知$f$求出$g$
这里有三种方法
方法1
for (int i = 1; i <= n; i++) g[i] = 0;
for (int i = 1; i <= n; i++) {
for (int j = 1; i * j <= n; j++) {
g[i * j] = (g[i * j] + mu[i] * f[j]) % mod;
}
}
// 依照定义，O(nlogn)方法2
for (int i = 1; i <= N; i++) g[i] = f[i];
for (int i = 1; i <= N; i++) {
for (int j = 2; i * j <= N; j++) {
g[i * j] = (g[i * j] - g[i]) % MOD;
}
}
// 类似求狄利克雷卷积逆的方式,不需要线性筛 mu ,O(nlogn)方法3
for (int i = 1; i <= n; i++) g[i] = f[i];
for (auto i : p) {
for (int j = n / i; j >= 1; j--) {
g[i * j] = (g[i * j] - g[j]) % MOD;
}
}//O(nloglogn)第三种可以理解为$dp$，设$g(i,n)=\sum\limits_{d|n,d只含前i种质因子}f(d)\mu(\frac{n}{d})$
那么$g(i,n)=\begin{cases} g(i-1,n), \;\;\;p_i\nmid n \\g(i-1,n)-g(i-1,\frac{n}{p_i}),\;\;\;p_i|n \end{cases}$
筛法
埃氏筛
枚举每个数，筛去它们的整数倍。
复杂度$O(nloglogn)$ ，证明略。
线性筛（欧拉筛）
枚举每个数$x$，筛去它们的$p_i$倍（$p_i$为质数）。当$x\%p_i=0$时，不需要再往下筛直接$break$，因为$p_i$是递增的，所以$x$乘上其他的质数的结果一定会被$p_i$的倍数筛掉。显然这样每个数只会被筛一次，故复杂度为$O(n)$。
任何积性函数都可以进行线性筛，所有线性筛积性函数都是基于线性筛质数的。
当然完全积性函数也可以线性筛，并且更加容易。不管$i$与$p$是否互质，都有$f(i\cdot p)=f(i) \cdot f(p)$，直接筛就行了。
若想线性筛出积性函数$f(x)$，就必须快速计算出以下函数值，$f(1),f(p),f(p^k)$，其中$p$为质数。
考虑筛的过程中，$i\cdot p$会被$i$乘上$p$给筛掉，将$i$唯一分解得到$p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}$，其中$p_1\le p_2\le \cdots\le p_k$
则一定有$p\leq p_1$，这是显然的，因为一旦第一次遇到$i\%p=0$也就是最小质因子时，直接$break$了。
如果$p< p_1$，那么$p$与$i$互质，可以直接得到$f(i\cdot p)=f(i)f(p)$
如果$p=p1$，这是需要对$i$记录一个$low_i$，表示$i$中最小质因子的指数次幂，即$low_i=p_1^{c_1}$
如果$i$除掉$low_i$，那么结果中最小质因子一定大于$p_1$，从而得到$gcd(\frac{i}{low_i},low_i\cdot p)=1$，那么$f(i\cdot p)=f(\frac{i}{low_i})f(low_i\cdot p)$
举例
欧拉函数$\varphi$，$\varphi(1)=1$，$\varphi(p)=p-1$。当$i\%p\neq0$时，$i$和$p$互质，$\varphi(i\cdot p)=\varphi(i)\cdot \varphi(p)=\varphi(i)\cdot(p-1)$；当$i\%p=0$时，$\varphi(i\cdot p)=\varphi(\frac{i}{p^c})\varphi(p^{c+1})=\frac{\varphi(i)}{p^c(1-\frac1{p})}p^{c+1}(1-\frac{1}{p})=p\cdot \varphi(i)$
查看代码
void init(int mod = MOD) {
np[0] = np[1] = 1;phi[1] = 1;
for (int i = 2;i < N;i++) {
if (!np[i]) {
p.push_back(i);
phi[i] = i - 1;
}
for (auto j : p) {
if (i * j < N) {
np[i * j] = 1;
if (i % j == 0) {
phi[i * j] = phi[i] * j;
break;
}
else {
phi[i * j] = phi[i] * (j - 1);
}
}
else {
break;
}
}
}
}莫比乌斯函数$\mu$，$\mu(1)=1$，$\mu(p)=-1$。当$i\%p\neq0$时，$i$和$p$互质，$\mu(i\cdot p)=\mu(i)\cdot \mu(p)=-\mu(i)$；当$i\%p=0$时，$\mu(i\cdot p)=\mu(\frac{i}{p^c})\mu(p^{c+1})=\mu(\frac{i}{p^c})\cdot 0=0$
查看代码
void init(int mod = MOD) {
np[0] = np[1] = 1;mu[1] = 1;
for (int i = 2;i < N;i++) {
if (!np[i]) {
p.push_back(i);
mu[i] = -1;
}
for (auto j : p) {
if (i * j < N) {
np[i * j] = 1;
if (i % j == 0) {
mu[i] = 0;
break;
}
else {
mu[i * j] = -mu[i];
}
}
else {
break;
}
}
}
}约数个数函数$d$，$d(1)=1$，$d(p)=2$ 。当$i\%p\neq0$时，$i$和$p$互质，$d(i\cdot p)=d(i)\cdot d(p)=2d(i)$；当$i\%p=0$时，$d(i\cdot p)=d(\frac{i}{p^c})d(p^{c+1})=\frac{d(i)}{c+1}(c+2)$
约数和函数$\sigma$，$\sigma(1)=1$，$\sigma(p)=p+1$，当$i\%p\neq0$时，$i$和$p$互质，$\sigma(i\cdot p)=\sigma(i)\cdot \sigma(p)=(p+1)\sigma(i)$；当$i\%p=0$时，$\sigma(i\cdot p)=\sigma(\frac{i}{p^c})\sigma(p^{c+1})=\frac{\sigma(i)}{\sum\limits_{i=0}^{c}p^i}\sum\limits_{i=0}^{c+1}p^i$
查看代码
void init(int mod = MOD) {
np[0] = np[1] = 1;
sigma[1] = 1;low[1] = 1;
for (int i = 2;i < N;i++) {
if (!np[i]) {
p.push_back(i);
sigma[i] = i + 1;low[i] = i;
}
for (auto j : p) {
if (i * j < N) {
np[i * j] = 1;
if (i % j == 0) {
low[i * j] = low[i] * j;
if (i == low[i]) sigma[i * j] = sigma[i] * j + 1;
else sigma[i * j] = sigma[i / low[i]] * sigma[low[i] * j];
break;
}
else {
low[i * j] = j;
sigma[i * j] = sigma[i] * (j + 1);
}
}
else {
break;
}
}
}
}下面这个不是积性函数，但是可以用来求解上面的东西，并且也可以线性筛。
最小次幂函数$low$，即$n=p_1^{c_1}p_2^{c_2}...p_k^{c_k}$，$p_1<p_2<...<p_k$，那么$low(n)=p_1^{c_1}$。
当$i\%p\neq0$时，$i$中最小质因子大于$p$，$low(i\cdot p)=p$；当$i\%p=0$时，$i$中最小质因子为$p$，$low(i\cdot p)=p\cdot low(i)$
杜教筛
模板题：P4213 - 杜教筛
在亚线性时间内求出某些特征积性函数的前缀和
所谓某些特征积性函数，指的是存在一个对应的积性函数$g$，使得$f*g=h$，且$g$，$h$的前缀和可以$O(1)$得知的函数$f$ 。
记$F,G,H$分别为$f,g,h$的前缀和，那么$H(n)=\sum\limits_{i=1}^{n}\sum\limits_{d|n}f(\frac{i}{d})g(d)=\sum\limits_{d=1}^{n}g(d)\sum\limits_{j=1}^{\lfloor \frac{n}{d} \rfloor}f(j)$
$=\sum\limits_{d=1}^{n}g(d)F(\lfloor \frac{n}{d} \rfloor)=g(1)F(n)+\sum\limits_{d=2}^{n}g(d)F(\lfloor \frac{n}{d} \rfloor)$
由于$g$为积性函数，那么$g(1)=1$，于是就有$F(n)=H(n)-\sum\limits_{d=2}^{n}g(d)F(\lfloor \frac{n}{d} \rfloor)$
注意到，如果从小到大筛的话，右边的$F$是已知的，而$G,H$都可以$O(1)$求得，则对$\lfloor\frac{n}{d} \rfloor$整除分块，单论复杂度$O(\sqrt{n})$
事实上，我们一般只关心一个点$n$的$F(n)$，故该式子一般递归求解。预处理$n^{\frac{2}{3}}$内的$f$值时，可以取得最优复杂度。
举例
求$\mu$的前缀和，注意到$\mu*1=\varepsilon$，于是有$\sum\limits_{i=1}^{n}\mu(i)=\sum\limits_{i=1}^{n}\varepsilon(i)-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}\mu(i)=1-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}\mu(i)$
求$\varphi$的前缀和，注意到$\varphi*1=id$，$\sum\limits_{i=1}^{n}\varphi(i)=\sum\limits_{i=1}^{n}id(i)-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}\varphi(i)=\frac{n(n+1)}{2}-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}\varphi(i)$
求$id\cdot\mu$的前缀和，我们需要找到合适的函数$g$，注意到$id*id=\sum\limits_{d|n}d\frac{n}{d}=\sum\limits_{d|n}n=nd(n)$
不妨令$g=id$尝试一下。于是有$(id\cdot \mu)*id=\sum\limits_{d|n}d\mu(d)\frac{n}{d}=n\sum\limits_{d|n}\mu(d)=n(\mu*1)=n\cdot \varepsilon(n)=\varepsilon(n)$
所以$\sum\limits_{i=1}^{n}id\cdot\mu(i)=\sum\limits_{i=1}^{n}\varepsilon(i)-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}id\cdot\mu(i)=1-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}id\cdot\mu(i)$
求$id\cdot\varphi$的前缀和，我们需要找到合适的函数$g$，与上面等同，$g=id$，于是有$(id\cdot \varphi)*id=\sum\limits_{d|n}d\varphi(d)\frac{n}{d}=n\sum\limits_{d|n}\varphi(d)=n(\varphi*1)=n\cdot id(n)=n^2=id_2(n)$
所以$\sum\limits_{i=1}^{n}id\cdot\varphi(i)=\sum\limits_{i=1}^{n}id_2(i)-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}id\cdot\varphi(i)=\frac{n(n+1)(2n+1)}{6}-\sum\limits_{d=2}^{n}\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor}id\cdot\varphi(i)$
杜教筛板子
int Du_sieve(int n, int pre_f[], map<int, int>& mp) {//f*g=h,求f的前缀和
if (n < N) return pre_f[n];
if (mp.count(n)) return mp[n];
int res = n * (n + 1) / 2;//h=id前缀和
for (int l = 2, r;l <= n;l = r + 1) {
r = n;
if (n / l) r = min(r, n / (n / l));
res = (res - Du_sieve(n / l, pre_f, mp) * (r - l + 1));//(r-l+1)为g=1前缀和
}
return mp[n] = res;
}
min_25筛
$Min25$筛求解前缀和本质上是一种$DP$。一般用于求解积性函数前缀和。
默认$P$为质数集，$P_i$表示第$i$个质数，默认$p\in P$。
若函数$f(x)$满足：$f(x)$是积性函数，$f(p)$存在多项式表示或可以快速求值，$f(p^k)$可以快速求值。
核心思想：将$[1,n]$分为质数，合数，$1$三个部分进行讨论。
$\sum\limits_{i=1}^{n}f(i)=\sum\limits_{i=1}^{n}[i\in P]f(i)+\sum\limits_{i=1}^{n}[i\notin P]f(i)$
$=\sum\limits_{i=1}^{n}[i\in P]f(i)+f(1)+\sum\limits_{2\le p\le \sqrt n,2\le p^e\le n}f(p^e)\left(\sum\limits_{i=2}^{\lfloor\frac{n}{p^e} \rfloor}f(i)[minp(i)>p]\right)$
首先考虑如何求解$\sum\limits_{i=1}^{n}[i\in P]f(i)$，考虑先线性筛出不大于$\sqrt n$的所有质数作为$P$。记$minp(x)$表示$x$的最小质因子，令$g(n,j)=\sum\limits_{i=1}^{n}[i\in P \; or \; minp(i)>P_j] i^k$，它表示$[1,n]$所有满足条件的数的$k$次方和，条件是要么是质数要么最小质因子大于$P_j$。这里选取$i^k$是因为它是一个完全积性函数，方便后面的公式推导。
可以发现$g(n,|P|)=\sum\limits_{i=1}^{n}[i\in P]i^k$，即$[1,n]$中所有质数的$k$次方和。为方便不妨令$g(n)=g(n,|P|)$。其中$|P|$表示$P$的大小，也就是不大于$\sqrt n$的质数数量。
考虑如何状态转移，$g(n,j)$在$g(n,j-1)$的基础上少了最小质因子恰好为$P_j$的合数。提出来一个$P_j$作为最小质因子，这样剩下的数也不能有小于它的质因子。也就是$g(\lfloor\frac{n}{P_j} \rfloor,j-1)-\sum\limits_{i=1}^{j-1}P_i^k$，后面这个是为了去掉前面满足条件的质数。这样我们可以得到递推式：$g(n,j)=g(n,j-1)- P_j^k\left[g(\lfloor\frac{n}{P_j} \rfloor,j-1)-\sum\limits_{i=1}^{j-1}P_i^k\right]$，式子中$\sum\limits_{i=1}^{j-1}P_i^k$，可以在筛素数的时候一并处理。
答案就是先求出所有质数的函数和，然后先枚举一个$p^e$，再枚举最小质因子大于$p$的数。
还是可以考虑$DP$的思想。设$S(n,x)$表示$[1,n]$中所有最小质因子大于$P_x$的函数值$f$之和。答案就是$S(n,0)$。
我们将满足条件的数分成两部分，第一部分是大于$P_x$的质数，即$g(n)-\sum\limits_{i=1}^{x}P_i^k$。另一部分则是最小质因子大于$P_x$的合数，枚举最小质因子即可。
于是$S(n,x)=g(n)-\sum\limits_{i=1}^{x}P_i^k+\sum\limits_{p_k^e\le n ,k>x }f(p_k^e)\left(S\left(\lfloor\frac{n}{p_k^e}\rfloor,k\right)+[e\ne1]\right)$
查看代码
namespace min25 {
const int SIZE = 1e5 + 10;
const int EXP_SIZE = 2;
int primes[SIZE], minFac[SIZE], pfx[SIZE][EXP_SIZE], primesPt, lim;
ll g[SIZE * EXP_SIZE][EXP_SIZE], dsc[SIZE * EXP_SIZE], inv2, inv3;
array<int, 2> indx[SIZE]; // indx[x]: index of <x, n / x>
const array<int, 2> csts[EXP_SIZE] = { {-1, 1}, {1, 2} };
void initPrimes(int siz) {
inv2 = qp(2, MOD - 2, MOD); inv3 = qp(3, MOD - 2, MOD);
fill(minFac + 0, minFac + siz + 1, 0); primesPt = 0;
for (int i = 2; i <= siz; i++) {
if (minFac[i] == 0) minFac[i] = i, primes[++primesPt] = i;
for (int j = 1; j <= primesPt && primes[j] <= min(minFac[i], siz / i); j++) minFac[i * primes[j]] = primes[j];
}
//f(p)的多项式每一项的前缀和
for (int i = 1; i <= primesPt; i++) {
for (int e = 0; e < EXP_SIZE; e++) {//展开写可能变快
pfx[i][e] = Add(pfx[i - 1][e], qp(primes[i], csts[e][1], MOD));
}
}
}
//计算f(p)
const auto f = [](ll p) {
p %= MOD;
return p * (p - 1) % MOD;
//ll res = 0;
// for (int e = 0; e < EXP_SIZE; e++) {//展开写会快很多
// res = (res + csts[e][0] * qp(p, csts[e][1], MOD)) % MOD;
// }
// return res;
};
//i^k前缀和,次数高的话需要拉插
const auto sum = [](ll n, ll exp) {
n %= MOD;
if (exp == 0) return n;
ll res = n * (n + 1) % MOD * inv2 % MOD;
if (exp == 2) return res * ((n << 1) + 1) % MOD * inv3 % MOD;
return res;
};
ll sieve(ll x, ll pt, ll n) {
if (x <= 1 || primes[pt] > x) return 0;
int k = x <= lim ? indx[x][0] : indx[n / x][1];
ll res = 0;
// for (int e = 0; e < EXP_SIZE; e++) {//展开写会快很多
// res = (res + csts[e][0] * (g[k][e] - pfx[pt - 1][e]) % MOD + MOD) % MOD;
// }
res = Add(res, Dec(pfx[pt - 1][0], g[k][0]));
res = Add(res, Dec(g[k][1], pfx[pt - 1][1]));
for (int i = pt; i <= primesPt && 1ll * primes[i] * primes[i] <= x; i++) {
ll pk = primes[i], pk1 = 1ll * primes[i] * primes[i];
for (int e = 1; pk1 <= x; pk = pk1, pk1 *= primes[i], e++) {
res = (res + f(pk) * sieve(x / pk, i + 1, n) % MOD + f(pk1)) % MOD;
}
}
return (res + MOD) % MOD;
}
ll run(ll n) {
lim = sqrt(n); initPrimes(lim); int dscPt = 0;
for (ll l = 1, r; l <= n; l = r + 1) {
r = n / (n / l); ll v = n / l; dsc[dscPt] = v;
for (int e = 0; e < EXP_SIZE; e++)
g[dscPt][e] = sum(dsc[dscPt], csts[e][1]) - 1;
v <= lim ? indx[v][0] = dscPt : indx[n / v][1] = dscPt; dscPt++;
}
for (int i = 1; i <= primesPt; i++) {
for (int j = 0; j < dscPt && 1ll * primes[i] * primes[i] <= dsc[j]; j++) {
ll v = dsc[j] / primes[i];
int k = v <= lim ? indx[v][0] : indx[n / v][1];
// for (int e = 0; e < EXP_SIZE; e++) {//展开写会快很多
// g[j][e] = (g[j][e] - qp(primes[i], csts[e][1], MOD) * (g[k][e] - pfx[i - 1][e] + MOD) % MOD + MOD) % MOD;
// }
g[j][0] = Dec(g[j][0], mul(primes[i], Dec(g[k][0], pfx[i - 1][0])));
g[j][1] = Dec(g[j][1], mul(mul(primes[i], primes[i]), Dec(g[k][1], pfx[i - 1][1])));
}
}
return Add(sieve(n, 1, n), 1);
}
整除分块
前置
给定$b,c$ ，求$a$ 。
$a\cdot b \le c$，则有$a\le \lfloor \frac{c}{b} \rfloor$
$a\cdot b < c$，则有$a< \lceil \frac{c}{b} \rceil$ ，也就是$a\cdot b\le c-1$，$a \le \lfloor\frac{c-1}{b} \rfloor$
$a\cdot b \ge c$，则有$a\ge \lceil \frac{c}{b} \rceil$
$a\cdot b > c$，则有$a> \lfloor \frac{c}{b} \rfloor$ ，也就是$a\cdot b\ge c+1$，$a \ge \lceil\frac{c+1}{b} \rceil$
对于$\lfloor\frac{n}{i}\rfloor$类似的式子，对于一个左边界$l$，其值为$k=\lfloor\frac{n}{l}\rfloor$，右边界$r$即找满足$k\le \lfloor\frac{n}{i}\rfloor$的最大的$i$ ，也就是使得$ik\le n$最大的$i$，即$r=\lfloor\frac{n}{k}\rfloor$，即$r= \left\lfloor \frac{n}{\lfloor\frac{n}{l} \rfloor} \right \rfloor$
查看代码
int lim = n;
for (int l = 1, r;l <= lim;l = r + 1) {
r = lim;
if (n / l) r = min(r, n / (n / l));
//操作
}拓展
我对整除分块的理解：
首先用$\lfloor\frac{n}{l}\rfloor$得到$k$，然后$\frac{n}{l}$应该为$k.几$，也就是$k+\varepsilon_1$，其中$\varepsilon_1\in[0,1)$
接下来用$\frac{n}{k}$，如果是$\frac{n}{k+\varepsilon_1}$，就会得到$l$，然而去掉了$\varepsilon_1$，$k$减小了，于是得到的值就变大，得到了$r+\varepsilon_2$。其中$\varepsilon_2\in[0,1)$
于是可以得出的是，在$x\in[l,r+\varepsilon_2]$范围，如果用$\frac{n}{x}$都将得到$[k,k+\varepsilon_1]$范围的值，其中$x=r+\varepsilon_2$时，$\frac{n}{x}=k$；$x=l$时，$\frac{n}{x}=k+\varepsilon_1$。于是取个整，$x\in[l,r]$的范围，$\lfloor\frac{n}{x}\rfloor=k$
下面的都可以这样理解。
给定$n$，求$\lfloor\frac{n}{ai+b}\rfloor$的整除分块
先令$i'=ai+b$，那么$l'=al+b$，$r'=ar+b$。 令$k=\lfloor \frac{n}{l'} \rfloor$，$r'$是最大的$i'$使得$ki'\le n$，那么$r'=\lfloor \frac{n}{k} \rfloor$
带入$l'=al+b$，得到$k=\lfloor \frac{n}{al+b} \rfloor$ ，$r'=\left\lfloor \frac{n}{\lfloor \frac{n}{al+b} \rfloor} \right\rfloor$
$r'=ar+b=\left\lfloor \frac{n}{\lfloor \frac{n}{al+b} \rfloor} \right\rfloor$ ，那么$r=\left\lfloor\frac{\left\lfloor \frac{n}{\lfloor \frac{n}{al+b} \rfloor} \right\rfloor-b}{a}\right\rfloor$
给定$n$，求$\lfloor\frac{n}{i^2}\rfloor$的整除分块
令$i'=i^2$，那么$l'=l^2$，$r'=r^2$。令$k=\lfloor\frac{n}{l'}\rfloor$，$r'$是最大的$i'$使得$ki'\le n$，那么$r'=\lfloor\frac{n}{k}\rfloor=\left\lfloor\frac{n}{\lfloor\frac{n}{l^2}\rfloor}\right\rfloor$
$r'=r^2=\left\lfloor\frac{n}{\lfloor\frac{n}{l^2}\rfloor}\right\rfloor$ ，得到$r=\left\lfloor\sqrt{\left\lfloor\frac{n}{\lfloor\frac{n}{l^2}\rfloor}\right\rfloor}\right\rfloor$
给定$n$，求$\lceil\frac{n}{i}\rceil$的整除分块
由于$\lceil\frac{n}{i}\rceil=\lfloor\frac{n+i-1}{i}\rfloor$ ，令$k=\lfloor\frac{n+l-1}{l}\rfloor$，$r$是最大的$i$使得$ik\le n+i-1$，那么$r=\lfloor\frac{n-1}{k-1} \rfloor=\left\lfloor\frac{n-1}{\lfloor\frac{n+l-1}{l}\rfloor-1}\right \rfloor$
$log$分块
$\sum\limits_{i=1}^{n}\lfloor log_ai \rfloor$
查看代码
for (int l = 1, r, k = (int)(log(l) / log(a));l <= n;l = r + 1, k++) {
r = min(l * a - 1, n);
res += k * (r - l + 1);
}
积性函数大赏
积性函数：一个定义域为正整数$n$的算术函数$f(n)$，有如下性质：$f(1)=1$，且当$a,b$互质时，$f(ab)=f(a)f(b)$。如果一个数字$n=p_1^{c_1} p_2^{c_2} \cdots p_k^{c_k}$，那么$f(n)=f(p_1^{c_1})f(p_2^{c_2}) \cdots f(p_k^{c_k})$，那么研究一个积性函数$f$可以转化为研究$f(p^c)$
完全积性函数：$f(1)=1$，且对任意两个正整数$a,b$都有$f(ab)=f(a)f(b)$
性质1：两个积性函数的狄利克雷卷积必定是积性函数。
性质2：狄利克雷卷积满足交换律，结合律，分配律。交换律即$f*g=g*f$，也就是$\sum\limits_{d|n}f(d)g(\frac{n}{d})=\sum\limits_{d|n}f(\frac{n}{d})g(d)$。分配律即$(f+g)*h=f*h+g*h$。结合律即$(f*g)*h=f*(g*h)$
1.
单位函数：$\varepsilon(n)=[n=1]$ （完全积性）
2.
幂函数：$id_k(n)=n^k$ （完全积性）
3.
常数函数：幂函数中，$k=0$时， $id_0(n)=1(n)=1$ ，或作$I(n)=1$ （完全积性）
4.
恒等函数：幂函数中，$k=1$时， $id_1(n)=id(n)=n$ （完全积性）
5.
除数函数：$\sigma_k(n)$，$n$所有正因数的$k$次幂之和。 这里$k$可以是任何复数和实数。
求解：
先来看$k=1$的情况，也就是约数和。
$(1+p_1+...+p_1^{c_1})...(1+p_n+...+p_n^{c_n})=\prod\limits_{i=1}^{n}(\sum\limits_{j=0}^{c_j}p_i^j)$
每个括号里面就是每一个质数的任一个次幂，就是每个括号里面选一个数字乘起来，再把所有选择加起来。显然这就是约束之和。
那么$k$次幂，显然就是$\prod\limits_{i=1}^{n}(\sum\limits_{j=0}^{c_j}(p_i^j)^k)=\prod\limits_{i=1}^{n}(\sum\limits_{j=0}^{c_j}(p_i^k)^j)$，显然可以用等比数列求和。
用PR分解longlong范围质因子，然后用上式计算。
查看代码
void Solve(int TIME) {
int n, k;cin >> n >> k;
Prime_factor(n);
map<int, int> mp;
for (auto i : res) {
while (n % i == 0) mp[i]++, n /= i;
}
int res = 1;
for (auto [i, j] : mp) {
int p = i;p = qp(p, k);
res = res * (qp(p, j + 1) - 1) % MOD * inv(p - 1) % MOD;
}
cout << res << endl;
}6.
因数个数函数：除数函数中，$k=0$时，$\sigma_o(n)=d(n)=\prod\limits_{i=1}^{k}(c_i+1)$ ，$n$的正因数个数，$n=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}$
7.
因数和函数：除数函数中，$k=1$时，$\sigma_1(n)=\sigma(n)=\prod\limits_{i=1}^{k}\sum\limits_{j=0}^{c_i}p_i^j$ ，$n$的所有正因数之和，$n=p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k}$
8.
莫比乌斯函数：$\mu(n)=\begin {cases} 1,n=1\\0,n含有平方因子\\(-1)^k ,n=\prod\limits_{i=1}^{k}p_i,p_i为质数 \end{cases} $
性质：$\sum\limits_{d|n}\mu(d)=\varepsilon(n)$ 即$\mu*1=\varepsilon$ 。这意味着在狄利克雷卷积中，$\mu$是常数函数$1$的逆元。
证明：设$n=\prod\limits_{i=1}^{k}p_i^{c_i}$，$n'=\prod\limits_{i=1}^{k}p_i$ ，那么$\sum\limits_{d|n}\mu(d)=\sum\limits_{d|n'}\mu(d)=\sum\limits_{i=0}^{k}C_{k}^{i}(-1)^i=(1+(-1))^k=[k=0]=[n=1]$
扩展：$[gcd(i,j)=1]=\sum\limits_{d|gcd(i,j)}\mu(d)$ ，将$n$直接用$gcd(i,j)$替换即可
9.
欧拉函数：$\varphi(n)=\sum\limits_{i=1}^{n}[gcd(i,n)=1]$ ，表示$[1,n]$中与$n$互质的数的个数
计算公式： $n=\prod\limits_{i=1}^{k}p_i^{c_i}$，$p_i$为质数 。 $\varphi(n)=n(1-\frac1{p_1})(1-\frac1{p_2})\cdots(1-\frac1{p_k})$
感性的证明：以$12$举例，$12$有两个质因子$2,3$。在$[1,12]$中，有$\frac1{2}$的数是$2$的倍数，所以有$1-\frac1{2}$是数不是$2$的倍数，在这不是$2$的倍数的数中，有$\frac{1}{3}$的数是3的倍数，所以既不是$2$的倍数也不是$3$的倍数的数字有$12(1-\frac1{2})(1-\frac{1}{3})$ 。
性质1：若$p$是质数，$\varphi(p)=p-1$ 。
性质2：若$p$是质数，$n=p^k$，$\varphi(n)=p^k-p^{k-1}$ 。证明：$\varphi(n)=n(1-\frac{1}{p})=p^k(1-\frac1{p})=p^k-p^{k-1}$
性质3：$\sum\limits_{d|n} \varphi(d)=n $即 $\varphi*1=id$ 。
证明：$n=\prod\limits_{i=1}^{k}p_i^{c_i}$，由于$\varphi$为积性函数，故只需要证明当$n'=p^c$时，$\varphi*1=\sum\limits_{d|n'}\varphi(d)1(\frac{n'}{d}) = id$
显然$d=p^0,p^1,p^2,...,p^c$ ，那么$\sum\limits_{i=0}^{c}\varphi(p^i)=\sum\limits_{i=0}^{c}p^i(1-\frac{1}{p})=1+(1-\frac{1}{p})\sum\limits_{i=1}^cp^i=1+(1-\frac{1}{p})\frac{p(1-p^c)}{1-p}=p^c=id(n')$
狄利克雷函数与欧拉函数的联系：对$\varphi*1=id$两边同时卷上$\mu$，得到$\varphi*1*\mu=id*\mu$
由于狄利克雷卷积满足结合律，$\varphi*\varepsilon=id*\mu$即$\varphi=\mu*id$
比较有用的式子
根据定义 $\varphi(n)=\sum\limits_{i=1}^{n}[gcd(i,n)=1]=\sum\limits_{i=1}^{n-1}[gcd(i,n)=1]+\varepsilon(n)$
$\sum\limits_{i=1}^{n}\varphi(i)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i-1}[gcd(i,j)=1]+1$，即有序对的贡献再加上$1$ 。
那么$\sum\limits_{i=1}^{n}\varphi(i)=\frac12\left(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[gcd(i,j)=1]-\sum\limits_{i=1}^{n}[gcd(i,i)=1]\right)+1=\frac12\left(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[gcd(i,j)=1]+1\right)$
即$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[gcd(i,j)=1]=2\sum\limits_{i=1}^{n}\varphi(i)-1$
10.
$\omega(n)$表示$n$的质因子数目(不可重复) ，例如$n=p_1^{c_1} p_2^{c_2} \cdots p_k^{c_k}$ ，$\omega(n)=k$
11.
$\Omega(n)$表示$n$的质因子数目(可重复) ，例如$n=p_1^{c_1} p_2^{c_2} \cdots p_k^{c_k}$ ，$\Omega(n)=\sum\limits_{i=1}^{k}c_i$
12.
刘维尔函数：$\lambda(n)=(-1)^{\Omega(n)}$
性质：$\sum\limits_{d|n}\lambda(d)=\begin{cases} 1,若n是平方数\\0,若n不是平方数 \end{cases}$ ，$\lambda(n)=\sum\limits_{d^2|n}\mu(\frac{n}{d^2})$
13.
$\gamma(n)=(-1)^{\omega(n)}$
14.
最大公约数：$gcd(n,k)$当$k$固定时，是积性函数
莫比乌斯反演
因数形式
$f(n)=\sum\limits_{d|n}g(d) ⟺ g(n)=\sum\limits_{d|n}\mu(d)f(\frac{n}{d})$
即$f=g*1 ⟺ g=\mu*f$
证明：
两边同时乘以$\mu$得到$f*\mu=g*1*\mu=g*\varepsilon=g$
莫比乌斯容斥
跟莫比乌斯反演关系不大，但是因为很有用并且用到了莫比乌斯函数，所以也放在这里
从埃筛到莫比乌斯容斥
用$cnt_i$表示$1\sim n$中，$i$出现的次数，显然$cnt_i=1$
用$k_i$表示$1\sim n$中，$i$的倍数出现的次数。
如何用$k$来表示出$cnt$呢？
这里不妨设$n=10$。
1 2 3 4 5 6 7 8 9 10
1 1 1 1 1 1 1 1 1 1
-1 -1 -1 -1 -1
-1 -1 -1
-1 -1
1
-1
1$cnt_1=k_1-k_2-k_3-k_5+k_6-k_7+k_{10}$，首先筛$1$的倍数，给每个$1$的倍数都加$1$，然后筛$2$的倍数，给每个$2$的倍数减去$1$，然后筛$3$的倍数，也给每个$3$的倍数都减去$1$。然后发现$4$已经被筛过了，跳过。然后是$5$的倍数，也都减去$1$。接下来是$6$的倍数，发现它被重复减了一次，于是枚举$6$的倍数，加回来一次。接下来筛$7$的倍数，给每个都减去$1$。接下来发现$8$的倍数已经筛过了，于是跳过。然后$9$也筛过了，跳过。最后$10$的倍数，发现多减了一次，加回来。
找一下有什么规律？$2$，$3$，$5$，$7$都包含了奇数个质因子所以是减去，$6$和$10$则包含了偶数个质因子所以是加上，而$4$，$8$，$9$包含了相同的质因子多次，也就是有平方数因数。接下来特殊定义$1$的时候是加上。
这个时候发现它恰好满足莫比乌斯函数的定义。
多推点式子看看：$cnt_2=k_2-k_4-k_6-k_{10}=k_{2\cdot 1}-k_{2\cdot 2}-k_{2\cdot 3}-k_{2\cdot 5}$
于是得到一个通式：$cnt_x=\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}k_{x\cdot i}\mu_i$
做些推广就得到了莫比乌斯容斥：
设$f(x)$表示$1\sim n$中$x$出现的次数（$x$的贡献）
$g(x)$表示$1\sim n$中有多少个是$x$的倍数（$x$的倍数的贡献）
则有$f(x)=\sum\limits_{i=1}^{\lfloor \frac{n}{x} \rfloor}g(i\cdot x)\mu(i)$
例如：
对于$g(x)=\sum\limits_{i=1}^n\sum\limits_{j=1}^{n} [gcd(a_i,a_j)=x]$的形式，可以令$f(d)=\sum\limits_{x=d,2d,...}g(x)=\sum\limits_{i=1}^n\sum\limits_{j=1}^{n}][d|gcd(a_i,a_j)]=\sum\limits_{i=1}^n\sum\limits_{j=1}^{n}[d|a_i][d|a_j]=(\sum\limits_{i=1}^{n}[d|a_i])^2$
$f(x)$可以$O(nlogn)$时间复杂度求出。
于是就有$g(d)=\sum\limits_{x=d,2d,...} f(x)\mu(\frac{x}{d})$
即$g(x)=\sum\limits_{n=x,2x,...} f(n)\mu(\frac{n}{x})=\sum\limits_{n=x,2x,...}(\sum\limits_{i=1}^{n}[n|a_i])^2\mu(\frac{n}{x})$
也可以$O(nlogn)$复杂度求出。
欧拉反演
利用$id=\varphi*1$，即$F(x)=\sum\limits_{d|F(x)}\varphi(d)$
例如：
对于$f(x)=\sum\limits_{i=1}^n\sum\limits_{j=1}^{n}gcd(a_i,a_j)$的形式，可以得到$f(x)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{d|gcd(a_i,a_j)} \varphi(d)=\sum\limits_{i=1}^{n}[d|a_i]\sum\limits_{j=1}^{n}[d|a_j]\sum\limits_{d=1}^{n} \varphi(d)$
矩阵树定理
$Laplacian$矩阵
$g(i,i)$表示$i$的度数，入度+出度。
如果$i$和$j$ $(i\neq j)$之间有边的话，就给对应的位置$g(i,j)$填上$-1$，否则为$0$。如果有重边，那么就改成对应的负几即可。
无向图生成树个数
$Laplacian$矩阵任选一个$i$后，去掉$i$行$i$列的矩阵的行列式。
生成树权值
边权的乘积。那么可以用有重边的生成树来表示所有生成树的权值之和。这个很容易理解，对于重边的图，生成树个数显然是每个边数乘起来后再乘上去重后的生成树个数。
有向图的生成树
以$r$为根的生成树，并且$r$能到达所有点。
有向图的生成树个数
以$r$为根的生成树的个数是，在$Laplacian$矩阵去掉$r$行$r$列的矩阵的行列式。
会用到的板子
辗转相除求行列式
int g[N][N];
int det(int n) {
int ans = 1;
for (int i = 1;i <= n;i++) {
for (int j = i + 1;j <= n;j++) {
//消掉g[j][i]
int x = i, y = j;
while (g[x][i]) {
int t = g[y][i] / g[x][i];
for (int k = i;k <= n;k++) {
g[y][k] = (g[y][k] - t * g[x][k]) % MOD;
}
swap(x, y);
}
//g[x][i]=0
if (x == i) {
for (int k = i;k <= n;k++) {
swap(g[i][k], g[j][k]);
}
ans = -ans;
}
}
if (g[i][i] == 0){
return 0;
}
ans = ans * g[i][i] % MOD;
}
if (ans < 0) ans += MOD;
return ans;
}
void Solve(int TIME) {
int n, m;cin >> n >> m;
for (int i = 1;i <= m;i++) {
int u, v;cin >> u >> v;
g[u][v]--, g[v][u]--;
g[u][u]++, g[v][v]++;
}
}
超几何分布
$n$个物品，其中有$m$个次品。不放回的取$k$个物品，有$x$个次品的概率为$\frac{C_{m}^{x}C_{n-m}^{k-x}}{C_{n}^{k}}$。 从组合的角度很容易理解
01分数规划
给定两个数组，$a_i$表示选取$i$的收益，$b_i$表示选取$i$的代价，定义$x[i]=1$或$0$，表示选或不选，每种物品只有选或不选两种方案，求出$\frac{\sum\limits_{i=1}^{n}a_i x_i}{\sum\limits_{i=1}^{n}b_i x_i }$的最大值/最小值
二分
假设我们求最大值，二分一个答案$mid$，当存在一组$x_i$，满足$\frac{\sum\limits_{i=1}^{n}a_i x_i}{\sum\limits_{i=1}^{n}b_i x_i }>mid$，则说明当前的$mid$可以。
$\frac{\sum\limits_{i=1}^{n}a_i x_i}{\sum\limits_{i=1}^{n}b_i x_i }>mid \;⇒\; \sum\limits_{i=1}^{n}a_ix_i-mid\sum\limits_{i=1}^{n}b_ix_i>0 \;⇒\;\sum\limits_{i=1}^{n}x_i(a_i-mid \cdot b_i)>0$
只要求出上述式子的最大值就可以了，如果最大值比$0$大，说明$mid$是可行的，否则不可行。
如果要求最小值，只要比$0$小即可。
二分求最小值
const ld eps = 1e-8;
ld l = eps, r = 1e9;
while (r - l > eps) {
ld mid = (l + r) / 2;
vc<ld> g;
for (auto [x, y, a, b] : edge) {
g.push_back(a - b * mid);
}
sort(g.begin(), g.end());
ld res = 0;
for (int i = 0;i < k;i++) {//选k个物品
res += g[i];
}
if (res < eps) r = mid;//ans
else l = mid;
}二分求最大值
const ld eps = 1e-8;
ld l = eps, r = 1e9;
while (r - l > eps) {
ld mid = (l + r) / 2;
vc<ld> g;
for (auto [x, y, a, b] : edge) {
g.push_back(-(a - b * mid));
}
sort(g.begin(), g.end());
ld res = 0;
for (int i = 0;i < k;i++) {//取k个物品
res += g[i];
}
if (res > eps) l = mid;//ans
else r = mid;
}$Dinkelbach$算法
思想是每次用上一轮的答案当做新的$mid$来输入，不断地迭代，直至答案满足精度。比二分更快。
Dinkelbach求最小值
const ld eps = 1e-8;
ld ans = 1e9;
while (1) {
DSU d(n + 1);
vc<pair<ld, ai2>> g;
for (auto [x, y, a, b] : edge) {
g.push_back({ a - b * ans,{a,b} });
}
sort(g.begin(), g.end());
for (int i = 0;i < k;i++)//取k个物品
A += g[i].second.at(0), B += g[i].second.at[1];
ld t = 1.0 * A / B;
if (abs(t - ans) < eps) break;
ans = t;
}Dinkelbach求最大值
const ld eps = 1e-8;
ld ans = eps;
while (1) {
DSU d(n + 1);
vc<pair<ld, ai2>> g;
for (auto [x, y, a, b] : edge) {
g.push_back({ -(a - b * ans),{a,b} });
}
sort(g.begin(), g.end());
for (int i = 0;i < k;i++)//取k个物品
A += g[i].second.at(0), B += g[i].second.at[1];
ld t = 1.0 * A / B;
if (abs(t - ans) < eps) break;
ans = t;
}
求二元一次不定方程
给定$a,b,c$ ，求$ax+by=c$中$x$和$y$的整数解。
首先用exgcd求出一组特解$x_0,y_0$，然后有$ax_0+by_0=gcd(a,b)=d$，如果$c \;mod\;d \neq 0$则无解。
令$base=\frac{c}{d}$，则有$a(x_0\cdot base+\frac{b}{gcd(a,b)})+b(y_0\cdot base-\frac{a}{gcd(a,b)})=c$ 。则有通解为$\begin{cases} x=x_0 \cdot \frac{c}{d} +k\cdot\frac{b}{gcd(a,b)}\\y=y_0\cdot \frac{c}{d}-k\cdot\frac{a}{gcd(a,b)} \end{cases}$ 。则很显然只要对$x_0\cdot \frac{c}{d}$关于$\frac{b}{gcd(a,b)}$取模就能得到$x$的最小非负整数解，对$y$同理，只要对$y_0\cdot \frac{c}{d}$关于$\frac{a}{gcd(a,b)}$取模就能得到$y$的最小非负整数解
模板题：P5656 二元一次不定方程
如果是非负整数，则是完全背包和同余最短路相关。具体见图论。
$Lucas$定理
求$C_n^m \;mod\; P$，$P$为质数
将$n$写成$P$进制数，获得$n_0P^0+n_1P^1+n_2P^2+...$
将$m$写成$P$进制数，获得$m_0P^0+m_1P^1+m_2P^2+...$
则$C_n^m \;mod\; P=C_{n_0}^{m_0}\cdot C_{n_1}^{m_1}\cdot C_{n_2}^{m_2}... \;mod\;P$
阶，原根，指数方程
阶：满足$a^x \equiv 1(\mod m)$ 最小的正整数$x$，读作$a$模$m$的阶。当前仅当$gcd(a,m)=1$这个方程会有解。
由欧拉公式我们知道，当$\gcd(a,m)=1$，一定有$a^{\varphi(m)}=1 (\mod m)$
阶可以理解为，从$1$开始不停的乘$a(\mod m)$，一个循环的长度。
那么如果$a^x\equiv 1(\mod m)$，那么一定有$阶|x$ 。
因此我们知道，一定有$阶|\varphi(m)$，也就是说如果要回到循环原点，步数一定是环长的倍数。
于是过程是这样的，将$\varphi(m)$用唯一分解定理表示成$\varphi(m)=p_1^{e_1}p_2^{e_2}...p_k^{e_k}$，对于每个质因子，如果去掉一个之后还是阶的倍数，即$a^{\frac{x}{p_i}}=1(\mod m)$，那么就去掉。
原根：$g\mod m$的阶为$\varphi(m)$，那么$g$为$m$的原根。
结论：只有$1$，$2$，$4$，$p^a$，$2p^a$存在原根，这里的$p$为奇素数。
原根不是唯一的。最小的原根不会很大。原根的个数为$\varphi(p-1)$个。
从小到大枚举或随机$g$，然后判断是否是原根。跟上面差不多，分解为若干质数$p_i$后，判断每个$p_i$是否都有$g^{\frac{\varphi(m)}{p_i}} \neq 1(\mod a)$即可。
需要注意的是，$p$的原根不一定是$p^2$的原根
当$g$为原根时，由于$1$到$m$之间与$m$互质的数有$\varphi(m)$个， 又因为原根的性质得到循环节长度是$\varphi(m)$，又知道循环节里面每个数之间都是互不相同的。所以这里面的任意一个数都可以用$g^i$表示。
指数方程：求$a^x \equiv b(\mod m)$，这里$\gcd(a,m)=1$
由欧拉定理，我们知道这个方程如果有解，答案一定是在$0$到$\varphi(m)-1$之间，因此最暴力的做法就是for循环累乘，但是这样复杂度就是比较高的。
这里运用meet in the middle的思想，令$T=\lceil \sqrt{m} \rceil$，那么$0$到$m-1$的数字一定可以表示成$qT+r$的形式，其中$q\in[0,T-1]$，$r\in [0,T-1]$。那么这个方程可以表示成$a^{qT+r}\equiv b(\mod m)$。
移项可得$(a^T)^q \equiv b\cdot a^{-r} (\mod m)$。就相当于在两个数组$a^{T\cdot 0},a^{T\cdot 1},...,a^{T\cdot (T-1)}$和$b\cdot a^{-0},b\cdot a^{-1},...,b\cdot a^{-(T-1)}$里面找是否存在两个相同的数字。
代码中，为了方便起见，用$x=qT-r$表示，其中$q\in[1,T]$，$r\in [1,T]$
那么就会变成$a^{qT} \equiv b \cdot a^r $，这样这个$r$就不需要求逆元，方便一点。
查看代码
int bsgs(int a, int b, int mod) {
a %= mod, b %= mod;
int T = sqrt(mod) + 2;
int v = qp(a, T, mod), d = 1;
unordered_map<int, int> hs;
for (int q = 1;q <= T;q++) {
d = d * v % mod;
if (!hs.count(d)) hs[d] = q * T;
}
int ans = mod + 1;
d = b;
for (int r = 1;r <= T;r++) {
d = d * a % mod;
if (hs.count(d)) ans = min(ans, hs[d] - r);
}
if (ans >= mod) return -1;
else return ans;
}
容斥原理
容斥原理的基本形式
$|S-\bigcup\limits_{i=1}^{n}A_i|=\sum\limits_{i=0}^{n}(-1)^{i}\sum\limits_{1\le j_1 \le j_2 \le...\le j_i \le n} |\bigcap\limits_{k=1}^{i}A_{j_k}|$
容斥原理的符号形式
设$S$是一个有限集，$a_1,...,a_n$是$n$种性质。
记$N(a_i)$为$S$中有$S$中有$a_i$性质的元素的数量。特殊的，记$N(1)=|S|$。
记$N(1-a_i)$为$S$中没有$a_i$性质的元素的数量。
$N(a_{i_1}a_{i_2}...a_{i_k})$为$S$中同时有$a_{i1},a_{i2},...,a_{i_k}$性质的元素数量。
记$N(a\pm b)=N(a) \pm N(b)$
那么容斥原理可以写成
$N((1-a_1)(1-a_2)...(1-a_n))=\sum\limits_{i=0}^{n}(-1)^i \sum\limits_{1\le j_1 <j_2<...<j_i \le n}N(a_{j_1}a_{j_2}...a_{j_i})$
$N(a_1...a_x(1-a_{x+1})...(1-a_{x+n}))=\sum\limits_{i=0}^{n}(-1)^i\sum\limits_{x \le j_1 <j_2<...<j_i \le x+n} N(a_1...a_x a_{j_1}...a_{j_i})$
而交集可以表示为$N(a_1a_2...a_n)$，并集表示为$N(1-(1-a_1)(1-a_2)...(1-a_n))$表示至少有一个满足条件的补集。
表达成了代数形式，可以简化一些问题的推导过程。
例题1:求正整数$1$到$n$中，既不是2的倍数而不是3的倍数，但是是5的倍数的数的数量。
设$a_1:是2的倍数 ,a_2:是3的倍数,a_3:是5的倍数$
$N((1-a_1)(1-a_2)a_3)=N(a_3-a_1a_3-a_2a_3+a_1a_2a_3)$
$=N(a_3)-N(a_1a_3)-N(a_2a_3)+N(a_1a_2a_3)=\lfloor\frac{n}{5}\rfloor-\lfloor\frac{n}{10}\rfloor-\lfloor\frac{n}{15}\rfloor+\lfloor\frac{n}{30}\rfloor$
例题2：刚好存在$k$个$i$满足$p_i=i$的排列个数。
设$a_i: \; p_i=i$
$N(a_{1}...a_{k}(1-a_{k+1})...(1-a_n))=\sum\limits_{i=0}^{n-k}(-1)^iC_{n-k}^{i}N(a_1...a_ka_{j_i}...a_{j_i})$
$=\sum\limits_{i=0}^{n-k}(-1)^iC_{n-k}^{i}(n-(k+i))!$
又因为前面这$k$个是任取的，所以要总体再乘上$C_{n}^{k}$ 。答案就是$C_n^k\sum\limits_{i=0}^{n-k}(-1)^iC_{n-k}^{i}(n-(k+i))!$
这个形式可以窥见一些二项式反演的迹象，有时候使用比二项式反演还要方便，但有时候还是二项式反演方便。建议把两个方法都熟练掌握。
位运算常见结论
$a|b + a\&b =a+b$
$a|b-a\&b=a\text{^} b$
$a-b=2a-(a+b)=2a-a|b-a\&b$
可以从集合的角度画个图理解：两个集合的总和是两个集合的交集加上并集，两个集合的异或是两个集合的并集减去两个集合的交集。
裴蜀定理
$a,b$是不全为$0$的整数，对任意整数$x,y$，满足$gcd(a,b) \;|\; ax+by$，且存在整数$x,y$，使得$ax+by=gcd(a,b)$。
推论：在长为$n$的循环数组（即数组中的最后一个元素的下一个元素是数组中的第一个元素，数组中第一个元素的前一个元素是数组中的最后一个元素）中，如果有周期$n$和$k$（$n$是循环数组自带的性质导致的），有$a_i=a_{i+nx+ky}=a_{i+gcd(n,k)}$ 。
勒让德定理：
$n!$做质因子分解后，质因子$p$的次数为$\sum\limits_{k=1}\lfloor\frac{n}{p^k} \rfloor$。
证明：这$n$个数的连乘积中，只有$p、2p、...、\lfloor \frac{n}{p} \rfloor \cdot p$能被$p$整除。从而有了$\lfloor \frac{n}{p} \rfloor$个质因子$p$ 。只有$p^2、2p^2、...、\lfloor \frac{n}{p^2} \rfloor\cdot p^2$，考虑到每个能被$p^2$整除的数都被$p$整除，所以只需要加上$\lfloor\frac{n}{p^2}\rfloor$个即可，以此类推即是结果。
$[n,2n]$之间必然存在质数。
证明很复杂，但是可以用哥德巴赫猜想做伪证明。显然$2n$是个偶数，那么右哥猜得到存在两个质数$p+q=2n$，如果$[n,2n]$不存在质数，那么$p<n,q<n$，不成立。
加强：伯特兰-切比雪夫定理，若整数$n>3$，则至少存在一个质数$p$，符合$p\in(n,2n-2)$。
线性规划
最大化目标函数：
需要最大化$\boldsymbol {C^{T}X} $，即$c_1x_1+...+c_mx_m$
其中$\boldsymbol X\ge 0$ $\boldsymbol {AX}\le \boldsymbol B$
其中$\boldsymbol A$是一个$n\times m$的矩阵，$\boldsymbol B$是一个$n$维向量，$\boldsymbol X$是一个$m$维向量。
$n$和$m$表示约束的个数和变量的个数。参数$\boldsymbol a$和$\boldsymbol b$表示线性规划问题的系数矩阵和约束向量。
线性代数
向量空间
向量空间是由向量组成的线性空间，定义$(F,V,+,\cdot)$为向量空间，其中$F$为域，$V$为集合，$V$中元素称为向量，$+$为向量加法，$\cdot $为标量乘法。
线性相关
对于向量空间中$V$上$n$个元素的向量组$(\boldsymbol{v}_{1},...,\boldsymbol{v}_{n})$，若存在不全为$0$的数$a_i\in F$，满足$a_1\boldsymbol v_1+a_2\boldsymbol v_2+...+a_n\boldsymbol v_n=0$，则称这$n$个向量线性相关，否则称为线性无关。
线性组合
对于向量空间中$V$上$n$个元素的向量组$(\boldsymbol{v}_{1},...,\boldsymbol{v}_{n})$，其线性组合是如下形式的向量，$a_1\boldsymbol v_1+a_2\boldsymbol v_2+...+a_n\boldsymbol v_n$，其中$a_1,...,a_n\in F$。一组向量线性无关$\Longleftrightarrow$没有向量可用有限个其他向量线性组合所表示
张成
对于向量空间中$V$上$n$个元素的向量组$(\boldsymbol{v}_{1},...,\boldsymbol{v}_{n})$，其所有线性组合所构成的集合称为$(\boldsymbol{v}_1,...,\boldsymbol{v}_{n})$的张成，记为$span(\boldsymbol{v}_{1},...,\boldsymbol{v}_{n})$。
基
若向量空间$V$中向量组$\mathcal B$既是线性无关的又可以张成$V$，则称其为$V$的基。
$\mathcal B$中的元素称为基向量，如果基中的元素个数有限，就称向量空间为有限维向量空间，将元素的个数称作向量空间的维数。
性质：设$\mathcal B$是向量空间$V$的基，则$\mathcal B$具有以下性质：
1.$V$是$\mathcal B$的最小生成集，就是说只有$\mathcal B$能张成$V$，而它的任何真子集都不张成全部的向量空间。
2.$\mathcal B$是$V$中线性无关向量的极大集合，就是说$\mathcal B$在$V$中是线性无关集合，而且$V$中没有其他线性无关集合包含它作为真子集。
3.$V$中所有向量都可以按唯一的方式被表达为$\mathcal B$中向量的线性组合。
线性相关性引理
如果$(\boldsymbol v_1,...,\boldsymbol v_n)$在$V$中是线性相关的，并且$\boldsymbol v_1\ne 0$，则有至少一个$j\in \{2,...,n\}$使得下列成立：
1. $\boldsymbol v_j\in span(\boldsymbol v_1,...,\boldsymbol v_{j-1})$
2. 如果从$(\boldsymbol v_1,...,\boldsymbol v_n)$去掉第$j$项，则剩余向量组的张成仍然等于$span(\boldsymbol v_1,...,\boldsymbol v_n)$
下面介绍算法竞赛中的线性基
异或线性基
对于数$a_0,a_1,...,a_n$将$a_i$的二进制表示$(b_m...b_0)_2$看作一个向量$\boldsymbol a_i=(b_m,...,b_0)$。
下文称向量$\boldsymbol a_i$的第$j$位为$\boldsymbol a_{i,j}$。
向量组$\boldsymbol a_1,...,\boldsymbol a_n$可以张成一个向量集合$span(\boldsymbol a_1,...,\boldsymbol a_n)$，加上异或运算和乘法运算，即可形成一个向量空间$V=(\{0,1\},span(\boldsymbol a_1,...,\boldsymbol a_n),\text{^} ,\cdot )$。
性质：对于任意存在于线性基的二进制位$i$，至多只有一个$\boldsymbol a$满足第$i$位为$1$。
模拟这个过程，$n=5$，$a={7,1,4,3,5}$。
一开始矩阵是$\left[\begin{array}{l} 0&0&0 \\0&0&0\\0&0&0 \end{array}\right]$，
加入$7=(111)_2$，矩阵变成：$\left[\begin{array}{l} 1&1&1 \\0&0&0\\0&0&0 \end{array}\right]$
加入$1=(001)_2$，添加到最后一行，同时消去第一行的最低位，矩阵变成：$\left[\begin{array}{l} 1&1&0 \\0&0&0\\0&0&1 \end{array}\right]$，
加入$4=(100)_2$，由于第一行已经有数字，它被异或为$(010)_2$，加入第二行，同时消去第一行的第二位，矩阵变成：$\left[\begin{array}{l} 1&0&0 \\0&1&0\\0&0&1 \end{array}\right]$.
剩下的数字都加不上了。
威尔逊定理
$(p-1)! \equiv-1(mod \;p)$是$p$为质数的充分必要条件。
验证充分性：
非质数的情况，带入$1$发现结果为$0$，带入$4$结果为$2$。
当$p>4$时，如果$p$是完全平方数，那么$p=k^2$，由于$p>4$，可得$k>2$，即$2k-p=2k-k^2=-(k-1)^2+1<0$，于是$2k<p$成立。
故$(p-1)!=1\times2\times ...\times k\times...\times 2k\times ...\times (p-1)=2k^2x=2px$，于是$(p-1)! \equiv 0(mod\; p)$。
当$p$不是完全平方数，那么$p$就必然等于两个完全不相等的数$a,b$的乘积，不妨设$a<b$，则有$1<a<b<p$，于是$(p-1)!=1\times ...\times a\times...\times b \times...\times (p-1)=abx=px$，于是也有$(p-1)!\equiv 0(mod \;p)$。
验证必要性：
当$p$为质数，考虑二次剩余式，$x^2\equiv 1(mod\;p)$ ，移项$x^2-1\equiv 0(mod\;p)$，因式分解$(x+1)(x-1)\equiv 0(mod\;p)$，可以得到$x\equiv p-1(mod\;p)$或者$x\equiv 1(mod\;p)$。抛开这两个数，对于$a\in [2,p-2]$，必然存在一个和它不相等的逆元$a^{-1}\in [2,p-2]$满足$aa^{-1}\equiv -1(mod\;p)$。于是必然有$\frac{p-3}{2}$对数相乘的乘积是$1$，即$(p-2)!\equiv 1(mod \;p)$，再乘上$(p-1)$即得到$(p-1)!\equiv p-1(mod\;p)$即$(p-1)!\equiv -1(mod\;p)$。
狄利克雷卷积
1.P2522 - Problem b
知识点
整除分块，狄利克雷卷积
题意
$\sum\limits_{i=a}^{b}\sum\limits_{j=c}^{d}[gcd(i,j)=k]$
题解
$\sum\limits_{i=a}^{b}\sum\limits_{j=c}^{d}[gcd(i,j)=k]=\sum\limits_{ki=a}^{b}\sum\limits_{kj=c}^{d}[gcd(ki,kj)=k]=\sum\limits_{i=\lceil \frac{a}{k}\rceil}^{\lfloor \frac{b}{k}\rfloor}\sum\limits_{j=\lceil \frac{c}{k}\rceil}^{\lfloor \frac{d}{k}\rfloor}[gcd(i,j)=1]$ （这里上限是下取整，下限是上取整。这样不会多枚举。如果反过来的话，就可能多枚举了。）
$=\sum\limits_{i=\lceil \frac{a}{k}\rceil}^{\lfloor \frac{b}{k}\rfloor}\sum\limits_{j=\lceil \frac{c}{k}\rceil}^{\lfloor \frac{d}{k}\rfloor}\sum\limits_{D|gcd(i,j)}\mu(D)=\sum\limits_{D=1}^{min(\lfloor \frac{b}{k}\rfloor,\lfloor \frac{d}{k}\rfloor)}\mu(D)\sum\limits_{i=\lceil \frac{a}{k}\rceil}^{\lfloor \frac{b}{k}\rfloor}\sum\limits_{j=\lceil \frac{c}{k}\rceil}^{\lfloor \frac{d}{k}\rfloor}[D|gcd(i,j)]$
$=\sum\limits_{D=1}^{min(\lfloor \frac{b}{k}\rfloor,\lfloor \frac{d}{k}\rfloor)}\mu(D)\sum\limits_{Di=\lceil \frac{a}{k}\rceil}^{\lfloor \frac{b}{k}\rfloor}\sum\limits_{Dj=\lceil \frac{c}{k}\rceil}^{\lfloor \frac{d}{k}\rfloor}[D|gcd(Di,Dj)]=\sum\limits_{D=1}^{min(\lfloor \frac{b}{k}\rfloor,\lfloor \frac{d}{k}\rfloor)}\mu(D)\sum\limits_{Di=\lceil \frac{a}{k}\rceil}^{\lfloor \frac{b}{k}\rfloor}\sum\limits_{Dj=\lceil \frac{c}{k}\rceil}^{\lfloor \frac{d}{k}\rfloor}1$
$=\sum\limits_{D=1}^{min(\lfloor \frac{b}{k}\rfloor,\lfloor \frac{d}{k}\rfloor)}\mu(D)\sum\limits_{i=\lceil \frac{a}{Dk}\rceil}^{\lfloor \frac{b}{Dk}\rfloor}\sum\limits_{j=\lceil \frac{c}{Dk}\rceil}^{\lfloor \frac{d}{Dk}\rfloor}1=\sum\limits_{D=1}^{min(\lfloor \frac{b}{k}\rfloor,\lfloor \frac{d}{k}\rfloor)}\mu(D)({\lfloor \frac{b}{Dk}\rfloor}-{\lceil \frac{a}{Dk}\rceil}+1)({\lfloor \frac{d}{Dk}\rfloor}-{\lceil \frac{c}{Dk}\rceil}+1)$
整除分块即可，简单证明一下两个整除分块。
$\lfloor \frac{n}{ik} \rfloor$，令$t=\lfloor \frac{n}{lk} \rfloor$ ，$r$为使得$itk\le n$最大的$i$，那么$r=\lfloor \frac{n}{tk} \rfloor=\lfloor\frac{n}{k\lfloor \frac{n}{lk} \rfloor} \rfloor$
$\lceil \frac{n}{ik} \rceil=\lfloor\frac{n+ik-1}{ik} \rfloor$，令$t=\lfloor \frac{n+lk-1}{lk} \rfloor$ ，$r$为使得$itk\le n+ik-1$最大的$i$，那么$r=\lfloor \frac{n-1}{tk-k} \rfloor=\lfloor\frac{n-1}{k\lfloor \frac{n+lk-1}{lk} \rfloor-k} \rfloor$
查看代码
void Solve(int TIME) {
int res = 0;
int t;
int a, b, c, d, k;cin >> a >> b >> c >> d >> k;
for (int l = 1, r;l <= min(b / k, d / k);l = r + 1) {
r = min(b / k, d / k);
t = (a + l * k - 1) / (l * k); if (k * t - k) r = min(r, (a - 1) / (k * t - k));
t = b / (l * k);if (k * t) r = min(r, b / (k * t));
t = (c + l * k - 1) / (l * k); if (k * t - k) r = min(r, (c - 1) / (k * t - k));
t = d / (l * k);if (k * t) r = min(r, d / (k * t));
res += (pre_mu[r] - pre_mu[l - 1]) * (1 + b / (l * k) - (a + l * k - 1) / (l * k)) * (1 + d / (l * k) - (c + l * k - 1) / (l * k));
}
cout << res << endl;
}2.P2257 - YY的GCD
知识点
狄利克雷卷积，整除分块
题意
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(x,y)\in Prime]$
题解
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[gcd(x,y)\in Prime]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d=1}^{min(n,m)}[gcd(i,j)=d][d \in Prime]$
$=\sum\limits_{d=1}^{min(n,m)}\sum\limits_{di=1}^{n}\sum\limits_{dj=1}^{m}[gcd(di,dj)=d][d \in Prime]=\sum\limits_{d=1}^{min(n,m)}\sum\limits_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d} \rfloor}[gcd(i,j)=1][d \in Prime]$
$=\sum\limits_{d=1}^{min(n,m)}\sum\limits_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d} \rfloor}\sum\limits_{D|gcd(i,j)}\mu(D)[d \in Prime]$
$=\sum\limits_{d=1}^{min(n,m)}\sum\limits_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d} \rfloor}\sum\limits_{D=1}^{min(\lfloor\frac{n}{d} \rfloor,\lfloor\frac{m}{d} \rfloor)}\mu(D)[d \in Prime][D|gcd(i,j)]$
$=\sum\limits_{D=1}^{min(\lfloor\frac{n}{d} \rfloor,\lfloor\frac{m}{d} \rfloor)}\mu(D)\sum\limits_{d=1}^{min(n,m)}[d \in Prime]\sum\limits_{Di=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{Dj=1}^{\lfloor\frac{m}{d} \rfloor}[D|gcd(Di,Dj)]$
$=\sum\limits_{d=1}^{min(n,m)}[d \in Prime]\sum\limits_{D=1}^{min(\lfloor\frac{n}{d} \rfloor,\lfloor\frac{m}{d} \rfloor)}\mu(D){\lfloor\frac{n}{dD} \rfloor}{\lfloor\frac{m}{dD} \rfloor}$ （到这一步已经可以整除分块，但是复杂度不够）
令$T=dD$，那么$\sum\limits_{T=1}^{min(n ,m)}{\lfloor\frac{n}{T} \rfloor}{\lfloor\frac{m}{T} \rfloor}\sum\limits_{d|T}[d \in Prime]\mu(\frac{T}{d})$
令$g(d)=[d\in Prime]$，$f(T)=\sum\limits_{d|T}g(d)\mu(\frac{T}{d})$，那么$\sum\limits_{T=1}^{min(n ,m)}{\lfloor\frac{n}{T} \rfloor}{\lfloor\frac{m}{T} \rfloor}\sum\limits_{d|T}g(d)\mu(\frac{T}{d})=\sum\limits_{T=1}^{min(n ,m)}{\lfloor\frac{n}{T} \rfloor}{\lfloor\frac{m}{T} \rfloor}f(T)$
这个$f(T)$显然是可以预处理的，然后便做完了。由于只用对质数计算贡献，复杂度约为$O(NloglogN)$，这是因为质数在自然数的占比约为$logn$，这样就对原本的$NlogN$进一步做了优化。
总时间复杂度$O(NloglogN+T(\sqrt{n}+\sqrt{m}))$ ，这里的$log$均以$e$为底
查看代码
int f[N];
void init(int mod = MOD) {
np[0] = np[1] = 1;mu[1] = 1;
for (int i = 2;i < N;i++) {
if (!np[i]) {
p.push_back(i);
mu[i] = -1;
}
for (auto j : p) {
if (i * j < N) {
np[i * j] = 1;
if (i % j == 0) {
break;
}
else {
mu[i * j] = -mu[i];
}
}
else {
break;
}
}
}
int t = 0;
for (int i = 1;i < N;i++) {
if (np[i]) continue;
for (int j = 1;i * j < N;j++) {
f[i * j] += mu[j];t++;
}
}
for (int i = 1;i < N;i++) {
f[i] += f[i - 1];
}
}
void Solve(int TIME) {
int n, m;cin >> n >> m;int res = 0;
int lim = min(n, m);
for (int l = 1, r;l <= min(n, m);l = r + 1) {
r = lim;
if (n / l) r = min(r, n / (n / l));
if (m / l) r = min(r, m / (m / l));
res += (n / l) * (m / l) * (f[r] - f[l - 1]);
}
cout << res << endl;
}3.CF1884D - Counting Rhyme
知识点
狄利克雷卷积
题解
注意$a_i\in[1,n]$，这个时候倍数的处理就显得轻松很多。只需要记录$[1,n]$哪些数字出现，然后在循环中枚举这个数字的倍数，打上标记。最坏时间复杂度是$nln(n)$，即$\frac{n}{1}+\frac{n}{2}+...+\frac{n}{n}=n(1+\frac{1}{2}+...+\frac{1}{n})=nln(n)$
设一个集合$S$，里面包含$[1,n]$中所有不为$a$数组中任何一个数的倍数的数。
$\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}\sum\limits_{x\in S}[gcd(a_i,a_j)=x]$
$=\frac{1}{2}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{x\in S}[gcd(a_i,a_j)=x]=\frac{1}{2}\sum\limits_{i=1}^{n}cnt_i\sum\limits_{j=1}^{n}cnt_j\sum\limits_{x\in S}[gcd(i,j)=x]$
$=\frac{1}{2}\sum\limits_{x\in S}[gcd(xi,xj)=x]\sum\limits_{xi=1}^{n}cnt_{xi}\sum\limits_{xj=1}^{n}cnt_{xj}=\frac{1}{2}\sum\limits_{x\in S}\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{xi}\sum\limits_{j=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{xj}[gcd(i,j)=1]$
$=\frac{1}{2}\sum\limits_{x\in S}\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{xi}\sum\limits_{j=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{xj}\sum\limits_{d|gcd(i,j)}\mu(d)=\frac{1}{2}\sum\limits_{x\in S}\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{xi}\sum\limits_{j=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{xj}\sum\limits_{d=1}^{\lfloor\frac{n}{x} \rfloor}\mu(d)[d|gcd(i,j)]$
$=\frac{1}{2}\sum\limits_{x\in S}\sum\limits_{di=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{dxi}\sum\limits_{dj=1}^{\lfloor\frac{n}{x}\rfloor}cnt_{dxj}\sum\limits_{d=1}^{\lfloor\frac{n}{x} \rfloor}\mu(d)[d|gcd(di,dj)]=\frac{1}{2}\sum\limits_{x\in S}\sum\limits_{d=1}^{\lfloor\frac{n}{x} \rfloor}\mu(d)\sum\limits_{i=1}^{\lfloor\frac{n}{dx}\rfloor}cnt_{dxi}\sum\limits_{j=1}^{\lfloor\frac{n}{dx}\rfloor}cnt_{dxj}$
令$T=dx$，那么$\frac{1}{2}\sum\limits_{x\in S}\sum\limits_{T=1}^{ n}\mu(\frac{T}{x})\sum\limits_{i=1}^{\lfloor\frac{n}{T}\rfloor}cnt_{Ti}\sum\limits_{j=1}^{\lfloor\frac{n}{T}\rfloor}cnt_{Tj}=\frac{1}{2}\sum\limits_{T=1}^{ n}\sum\limits_{ x|T}\mu(\frac{T}{x})[x\in S](\sum\limits_{i=1}^{\lfloor\frac{n}{T}\rfloor}cnt_{Ti})^2$
令$f(T)=\sum\limits_{ x|T}\mu(\frac{T}{x})[x\in S]$，那么就有$\frac{1}{2}\sum\limits_{T=1}^{n}f(T)(\sum\limits_{i=1}^{\lfloor\frac{n}{T}\rfloor}cnt_{Ti})^2$，$f(T)$可以暴力$nlogn$预处理，而后面的整体暴力求也是$O(nlogn)$复杂度
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi cnt(n + 1);
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i], cnt[a[i]] += 1;
vi not_in_s(n + 1);
for (int i = 1;i <= n;i++) {
if (cnt[i] == 0) continue;
for (int j = i;j <= n;j += i) {
not_in_s[j] = 1;
}
}
vi f(n + 1);
int res = 0;
for (int i = 1;i <= n;i++) {
if (not_in_s[i]) continue;
for (int j = 1;j * i <= n;j++) {
f[i * j] += mu[j];
}
}
for (int i = 1;i <= n;i++) {
int t = 0;
for (int j = 1;j <= n / i;j++) {
t += cnt[i * j];
}
res += t * t * f[i];
}
cout << res / 2 << endl;
}法2：利用莫比乌斯反演的第二形式
P.S.这种$[gcd(a_i,a_j)=x]$用这个方法会方便很多
$\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}\sum\limits_{x\in S}[ gcd(a_i,a_j)=x]$
$=\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}\sum\limits_{x\in S}\sum\limits_{x|D}[D|gcd(a_i,a_j)]\mu(\frac{D}{x})$
$=\sum\limits_{i=1}^{n}[D|a_i]\sum\limits_{j=i+1}^{n}[D|a_j]\sum\limits_{x\in S}\sum\limits_{x|D}\mu(\frac{D}{x})$
$=\sum\limits_{x\in S}\sum\limits_{x|D}\sum\limits_{i=1}^{n}[D|a_i]\sum\limits_{j=i+1}^{n}[D|a_j]\mu(\frac{D}{x})$
$F(D)=\sum\limits_{i=1}^{n}[D|a_i]\sum\limits_{j=i+1}^{n}[D|a_j]=\frac{1}{2}\left((\sum\limits_{i=1}^{n}[D|a_i])^2-\sum\limits_{i=1}^n[D|a_i]\right)$
原式$=\sum\limits_{x\in S}\sum\limits_{x|D}F(D)\mu(\frac{D}{x})$
查看代码
int n;cin >> n;
vi cnt(n + 1);
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i], cnt[a[i]] += 1;
vi not_in_s(n + 1);
for (int i = 1;i <= n;i++) {
if (cnt[i] == 0) continue;
for (int j = i;j <= n;j += i) {
not_in_s[j] = 1;
}
}
vi f(n + 1);
for (int i = 1;i <= n;i++) {
for (int j = 1;j * i <= n;j++) {
f[i] += cnt[i * j];
}
}
for (int i = 1;i <= n;i++) f[i] = (f[i] * f[i] - f[i]) / 2;
int res = 0;
for (int i = 1;i <= n;i++) {
if (not_in_s[i]) continue;
for (int j = 1;j * i <= n;j++) {
res += f[i * j] * mu[j];
}
}
cout << res << endl;
}
4.P3312 - 数表
知识点
狄利克雷卷积，整除分块
题解
设约数和函数$\sigma$
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sigma(gcd(i,j))[\sigma(gcd(i,j))\le a]=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d=1}^{min(n,m)}\sigma(d)[gcd(i,j)=d][\sigma(d)\le a]$
$=\sum\limits_{di=1}^{n}\sum\limits_{dj=1}^{m}\sum\limits_{d=1}^{min(n,m)}\sigma(d)[gcd(di,dj)=d][\sigma(d)\le a]=\sum\limits_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)[gcd(i,j)=1][\sigma(d)\le a]$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)\sum\limits_{D|gcd(i,j)}\mu(D)[\sigma(d)\le a]$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)\sum\limits_{D=1}^{min(\lfloor\frac{n}{d} \rfloor,\lfloor\frac{m}{d} \rfloor)}\mu(D)[D|gcd(i,j)][\sigma(d)\le a]$
$=\sum\limits_{Di=1}^{\lfloor\frac{n}{d} \rfloor}\sum\limits_{Dj=1}^{\lfloor\frac{m}{d} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)\sum\limits_{D=1}^{min(\lfloor\frac{n}{d} \rfloor,\lfloor\frac{m}{d} \rfloor)}\mu(D)[D|gcd(Di,Dj)][\sigma(d)\le a]$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{Dd} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{Dd} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)\sum\limits_{D=1}^{min(\lfloor\frac{n}{d} \rfloor,\lfloor\frac{m}{d} \rfloor)}\mu(D)[\sigma(d)\le a]$
令$T=Dd$，则有 $\sum\limits_{i=1}^{\lfloor\frac{n}{T} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{T} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)\sum\limits_{T=1}^{min(n,m)}\mu(\frac{T}{d})[\sigma(d)\le a]=\sum\limits_{T=1}^{min(n,m)}\mu(\frac{T}{d})\sum\limits_{i=1}^{\lfloor\frac{n}{T} \rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{T} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)[\sigma(d)\le a]$
$=\sum\limits_{T=1}^{min(n,m)}{\lfloor\frac{n}{T} \rfloor}{\lfloor\frac{m}{T} \rfloor}\sum\limits_{d=1}^{min(n,m)}\sigma(d)\mu(\frac{T}{d})[\sigma(d)\le a]$
令$g(T)=\sum\limits_{d=1}^{min(n,m)}\sigma(d)\mu(\frac{T}{d})[\sigma(d)\le a]$
那么有$\sum\limits_{T=1}^{min(n,m)}{\lfloor\frac{n}{T} \rfloor}{\lfloor\frac{m}{T} \rfloor}g(T)$ ，只有$\sigma(d)\le a$会对其产生贡献
多组询问，考虑将询问的$a$从小到大排序处理。再将$\sigma$按大小排序，进一步优化。
然后使用树状数组每次加入新的贡献，再计算这个询问的答案。
树状数组最坏情况需要加入$nlogn$次，每次复杂度是$logn$ ，复杂度$O(nlog^2n)$。
总时间复杂度$O(nlog^2n+q(\sqrt{n}logn))$
查看代码
void init(int mod = MOD) {
np[0] = np[1] = 1;mu[1] = 1;
sigma[1] = 1;low[1] = 1;
for (int i = 2;i < N;i++) {
if (!np[i]) {
p.push_back(i);
mu[i] = -1;
sigma[i] = i + 1;low[i] = i;
}
for (auto j : p) {
if (i * j < N) {
np[i * j] = 1;
if (i % j == 0) {
low[i * j] = low[i] * j;
if (i == low[i]) sigma[i * j] = sigma[i] * j + 1;
else sigma[i * j] = sigma[i / low[i]] * sigma[low[i] * j];
break;
}
else {
low[i * j] = j;
mu[i * j] = -mu[i];
sigma[i * j] = sigma[i] * (j + 1);
}
}
else {
break;
}
}
}
}
struct Q {
int a, id;
int n, m;
bool friend operator<(Q a, Q b) {
return a.a < b.a;
}
};
void Solve(int TIME) {
int q;cin >> q;
vc<Q> qu(q + 1);
for (int i = 1;i <= q;i++) {
int n, m, a;cin >> n >> m >> a;
qu[i] = { a,i,n,m };
}
BIT tr(1e5 + 10);
vi ans(q + 1);
vc<pi> Sigma(1e5 + 1);
for (int i = 1;i <= 1e5;i++) {
Sigma[i] = { sigma[i],i };
}
sort(Sigma.begin() + 1, Sigma.end());
sort(qu.begin() + 1, qu.end());
int now = 0;
for (int i = 1;i <= q;i++) {
auto [a, id, n, m] = qu[i];
int lim = min(n, m);
for (int j = now + 1;j <= 1e5;j++) {
if (Sigma[j].first > a) {
now = j - 1;
break;
}
for (int k = 1;k * Sigma[j].second <= 1e5;k++) {
tr.add(k * Sigma[j].second, Sigma[j].first * mu[k] % MOD);
}
}
for (int l = 1, r;l <= lim;l = r + 1) {
r = lim;
if (n / l) r = min(r, n / (n / l));
if (m / l) r = min(r, m / (m / l));
ans[id] = (ans[id] + tr.query(l, r) * (n / l) % MOD * (m / l) % MOD) % MOD;
}
}
for (int i = 1;i <= q;i++) cout << (ans[i] % MOD + MOD) % MOD << endl;
}5.P2398 - GCD SUM
知识点
莫比乌斯反演，整除分块
题解
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}gcd(i,j)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{k=1}^{min(n,m)}k[gcd(i,j)=k]$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor \frac{m}{k}\rfloor}\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d|gcd(i,j)}\mu(d)=\sum\limits_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum\limits_{j=1}^{\lfloor \frac{m}{k}\rfloor}\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d=1}^{min(\lfloor \frac{n}{k}\rfloor,\lfloor \frac{m}{k}\rfloor)}\mu(d)[d|gcd(i,j)]$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum\limits_{j=1}^{\lfloor \frac{m}{dk}\rfloor}\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d=1}^{min(\lfloor \frac{n}{k}\rfloor,\lfloor \frac{m}{k}\rfloor)}\mu(d)$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\sum\limits_{j=1}^{\lfloor \frac{m}{dk}\rfloor}\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d=1}^{min(\lfloor \frac{n}{k}\rfloor,\lfloor \frac{m}{k}\rfloor)}\mu(d)=\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d=1}^{min(\lfloor \frac{n}{k}\rfloor,\lfloor \frac{m}{k}\rfloor)}\mu(d)\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor$
$=\sum\limits_{T=1}^{min(n,m)}\mu(\frac{T}{k})\sum\limits_{k|T}k\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor$
$=\sum\limits_{T=1}^{min(n,m)}\varphi(T)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor$ ，（由$\mu*id=\varphi$得出）
推到这里才发现$n=m$，那么直接就是$\sum\limits_{T=1}^{n}\varphi(T)(\lfloor\frac{n}{T}\rfloor)^2$
查看代码
void Solve(int TIME) {
int n;cin >> n;
int res = 0;
for (int l = 1, r;l <= n;l = r + 1) {
r = n;
if (n / l) r = min(r, n / (n / l));
res += (pre_phi[r] - pre_phi[l - 1]) * (n / l) * (n / l);
}
cout << res << endl;
}
6.AGC038C - LCMs
知识点
狄利克雷卷积
题意
$\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}lcm(A_i,A_j)$
题解
$\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}lcm(A_i,A_j)=\frac12(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}lcm(A_i,A_j)-\sum\limits_{i=1}^{n}A_i)$
$=\frac12(\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}cnt_i \cdot cnt_j\cdot lcm(i,j)-\sum\limits_{i=1}^{N}i\cdot cnt_i)$
考虑计算$res1=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}cnt_i \cdot cnt_j\cdot lcm(i,j)=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}cnt_i \cdot cnt_j\cdot \frac{ij}{gcd(i,j)}$
$=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}\sum\limits_{k=1}^{N}cnt_i \cdot cnt_j\cdot \frac{ij}{k}[gcd(i,j)=k]=\sum\limits_{i=1}^{\lfloor \frac{N}{k} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{N}{k} \rfloor}\sum\limits_{k=1}^{N}cnt_{ki} \cdot cnt_{kj}\cdot \frac{ki\cdot kj}{k}[gcd(i,j)=1]$
$=\sum\limits_{i=1}^{\lfloor \frac{N}{k} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{N}{k} \rfloor}\sum\limits_{k=1}^{N}cnt_{ki} \cdot cnt_{kj}\cdot kij\sum\limits_{d|gcd(i,j)}\mu(d)=\sum\limits_{i=1}^{\lfloor \frac{N}{k} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{N}{k} \rfloor}\sum\limits_{k=1}^{N}cnt_{ki} \cdot cnt_{kj}\cdot kij\sum\limits_{d=1}^{\lfloor \frac{N}{k} \rfloor}\mu(d)[d|gcd(i,j)]$
$=\sum\limits_{i=1}^{\lfloor \frac{N}{dk} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{N}{dk} \rfloor}\sum\limits_{k=1}^{N}cnt_{dki} \cdot cnt_{dkj}\cdot d^2kij\sum\limits_{d=1}^{\lfloor \frac{N}{k} \rfloor}\mu(d)$
令$T=dk$，则有$\sum\limits_{i=1}^{\lfloor \frac{N}{T} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{N}{T} \rfloor}\sum\limits_{k|T}cnt_{Ti} \cdot cnt_{Tj}\cdot T\frac{T}{k}ij\sum\limits_{T=1}^{N}\mu(\frac{T}{k})$
$=\sum\limits_{T=1}^{N}T\sum\limits_{k|T}\mu(\frac{T}{k})\frac{T}{k}\sum\limits_{i=1}^{\lfloor \frac{N}{T} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{N}{T} \rfloor}cnt_{Ti} \cdot cnt_{Tj}\cdot ij=\sum\limits_{T=1}^{N}T\sum\limits_{k|T}\mu(\frac{T}{k})\frac{T}{k}(\sum\limits_{i=1}^{\lfloor \frac{N}{T} \rfloor}cnt_{Ti} \cdot i)^2$
$=\sum\limits_{T=1}^{N}T\sum\limits_{k|T}\mu(k)k(\sum\limits_{i=1}^{\lfloor \frac{N}{T} \rfloor}cnt_{Ti} \cdot i)^2$
最后的答案$\frac12(\sum\limits_{T=1}^{N}T\sum\limits_{k|T}\mu(k)k(\sum\limits_{i=1}^{\lfloor \frac{N}{T} \rfloor}cnt_{Ti} \cdot i)^2-\sum\limits_{i=1}^{N}i\cdot cnt_i)$
查看代码
void Solve(int TIME) {
int n;cin >> n;
for (int i = 1;i <= n;i++) {
int x;cin >> x;cnt[x] += 1;
}
for (int i = 1;i <= 1e6;i++) {
for (int j = 1;j * i <= 1e6;j++) {
f[i * j] = (f[i * j] + mu[i] * i % MOD) % MOD;
}
}
int res = 0;
for (int T = 1;T <= 1e6;T++) {
int t = 0;
for (int j = 1;j <= 1e6 / T;j++) {
t = (t + cnt[T * j] * j % MOD) % MOD;
}
res = (res + T * f[T] % MOD * t % MOD * t % MOD) % MOD;
}
for (int i = 1;i <= 1e6;i++) {
res = (res - 1ll * i * cnt[i] % MOD) % MOD;
}
res = res * inv(2) % MOD;
cout << (res % MOD + MOD) % MOD << endl;
}
或者改成下一题用莫反第二形式，代码如下
查看代码
void Solve(int TIME) {
int n;cin >> n;
const int M = 1e6;
vi cnt(M + 1);
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i], cnt[a[i]] += 1;
vi f(M + 1);
for (int i = 1;i <= M;i++) {
for (int j = i;j <= M;j += i) {
f[i] = (f[i] + cnt[j] * j % MOD) % MOD;
}
}
for (int i = 1;i <= M;i++) f[i] = f[i] * f[i] % MOD;
int res = 0;
for (int i = 1;i <= M;i++) {
int t = 0;
for (int j = 1;j * i <= M;j++) {
t = (t + f[i * j] * mu[j] % MOD) % MOD;
}
res = (res + t * inv(i) % MOD) % MOD;
}
for (int i = 1;i <= n;i++) res = (res - a[i] + MOD) % MOD;
cout << res * inv(2) % MOD << endl;
}
7.P3911 最小公倍数之和
知识点
莫比乌斯容斥
题意
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}lcm(a_i,a_j)$
题解
两种方法，一种是将式子的$a_i$转化为枚举值域然后记录数组$cnt_i$。
另一种方法是使用莫比乌斯反演的第二种形式。
上一题跟这个差不多，用的是法一。这里用第二种方法写一下。
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}lcm(a_i,a_j)=\sum\limits_{i=1}^{n}a_i\sum\limits_{j=1}^{n}a_j \frac{1}{gcd(a_i,a_j)}$
$=\sum\limits_{i=1}^{n}a_i\sum\limits_{j=1}^{n}a_j\sum\limits_{d=1}^{n} \frac{1}{d}[gcd(a_i,a_j)=d]$
$=\sum\limits_{i=1}^{n}a_i\sum\limits_{j=1}^{n}a_j\sum\limits_{d=1}^{n} \frac{1}{d}\sum\limits_{D=d,2d,...}[D|gcd(a_i,a_j)]\mu(\frac{D}{d})$
$=\sum\limits_{D=d,2d,...}\sum\limits_{i=1}^{n}a_i[D|a_i]\sum\limits_{j=1}^{n}a_j[D|a_j]\sum\limits_{d=1}^{n} \frac{1}{d}\mu(\frac{D}{d})$
$=\sum\limits_{D=d,2d,...}(\sum\limits_{i=1}^{n}a_i[D|a_i])^2\sum\limits_{d=1}^{n} \frac{1}{d}\mu(\frac{D}{d})$
令$F(D)=(\sum\limits_{i=1}^{n}a_i[D|a_i])^2$
原式$=\sum\limits_{d=1}^{n}\frac{1}{d} \sum\limits_{D=d,2d,...}F(D)\mu(\frac{D}{d})$
其中$F(x)$可以$O(nlogn)$预处理，而后面这个也可以$O(nlogn)$解决。
查看代码
void Solve(int TIME) {
int n;cin >> n;
const int M = 5e4;
vi cnt(M + 1);
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i], cnt[a[i]] += 1;
vi f(M + 1);
for (int i = 1;i <= M;i++) {
for (int j = i;j <= M;j += i) {
f[i] += cnt[j] * j;
}
}
for (int i = 1;i <= M;i++) f[i] = f[i] * f[i];
int res = 0;
for (int i = 1;i <= M;i++) {
int t = 0;
for (int j = 1;j * i <= M;j++) {
t += f[i * j] * mu[j];
}
res += t / i;
}
cout << res << endl;
}
8.HDU 7325 - GCD Magic
知识点
狄利克雷卷积，杜教筛
题意
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[\gcd(2^i-1,2^j-1)]^k$
题解
推公式或打表可以得出$gcd(2^{i}-1,2^{j}-1)=gcd(2^i-1,2^j-2^i)=gcd(2^i-1,2^i(2^{j-i}-1))$
$=gcd(2^i-1,2^{j-i}-1)=...=2^{gcd(i,j)}-1$
法1.
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[\gcd(2^i-1,2^j-1)]^k=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}(2^{gcd(i,j)}-1)^k$
$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{d=1}^{n}(2^{d}-1)^k[gcd(i,j)=d]=\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{d=1}^{n}(2^{d}-1)^k[gcd(i,j)=1]$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{d=1}^{n}(2^{d}-1)^k\sum\limits_{D|gcd(i,j)}\mu(D)=\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{d=1}^{n}(2^{d}-1)^k\sum\limits_{D=1}^{\lfloor\frac{n}{d} \rfloor}\mu(D)[D|gcd(i,j)]$
$=\sum\limits_{i=1}^{\lfloor\frac{n}{dD}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{dD}\rfloor}\sum\limits_{d=1}^{n}(2^{d}-1)^k\sum\limits_{D=1}^{\lfloor\frac{n}{d} \rfloor}\mu(D)[D|gcd(Di,Dj)]$
令$T=dD$，那么有$\sum\limits_{i=1}^{\lfloor\frac{n}{T}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{T}\rfloor}\sum\limits_{T=1}^{n}\sum\limits_{d|T}(2^{d}-1)^k\mu(\frac{T}{d})=\sum\limits_{T=1}^{n}(\lfloor\frac{n}{T}\rfloor)^2\sum\limits_{d|T}(2^{d}-1)^k\mu(\frac{T}{d})$
令$f(d)=(2^d-1)^k$，$g(T)=\sum\limits_{d|T}f(d)\mu(\frac{T}{d})$， 那么有$\sum\limits_{T=1}^{n}(\lfloor\frac{n}{T}\rfloor)^2g(T)$
现在即求$g=f*\mu$的前缀和$G$，考虑将$g$卷上$1$，则$\sum\limits_{i=1}^{n}f(i)=\sum\limits_{i=1}^{n}(g*1)(i)=\sum\limits_{i=1}^{n}\sum\limits_{j|i}g(j)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{\lfloor \frac{n}{i} \rfloor}g(j)=\sum\limits_{i=1}^{n}G(\lfloor \frac{n}{i} \rfloor)$
即有$G(n)=\sum\limits_{i=1}^{n}f(i)-\sum\limits_{i=2}^{n}G(\lfloor \frac{n}{i} \rfloor)$，这一步也就是杜教筛的结论。
然后计算$f$的前缀和，$F(n)=\sum\limits_{i=1}^{n}(2^i-1)^k=\sum\limits_{i=1}^{n}\sum\limits_{j=0}^{k}C_{k}^{j}2^{ij}(-1)^{k-j}=\sum\limits_{j=0}^{k}C_{k}^{j}(-1)^{k-j}\sum\limits_{i=1}^{n}(2^{j})^i$
后面是个等比数列，然后可以每次$klogn$查询$f$的前缀和
预处理出$1e6$范围的 $G$，然后杜教筛即可。复杂度$O(T(n^{\frac23}logn^\frac32+k\sqrt{n}logn))$
查看代码
int g[N];
int sum_f(int n, int k) {
int res = 0;
for (int j = 0;j <= k;j++) {
int sgn = (k - j) % 2 == 1 ? MOD - 1 : 1;
int q = (1 << j);
if (j == 0) {
res = (res + n * sgn % MOD) % MOD;
}
else {
res = add(res, sgn * comb(k, j) % MOD * (q * (qp(q, n) - 1) % MOD * inv(q - 1) % MOD) % MOD);
}
}
return norm(res);
}
int Du_sum_g(int n, int pre_f[], umap<int, int>& mp, int k) {
if (n < N) return pre_f[n];
if (mp.count(n)) return mp[n];
int res = sum_f(n, k);
for (int l = 2, r;l <= n;l = r + 1) {
r = n;
if (n / l) r = min(r, n / (n / l));
res = (res - Du_sum_g(n / l, pre_f, mp, k) * (r - l + 1) % MOD) % MOD;
}
return mp[n] = res;
}
void Solve(int TIME) {
int n, k;cin >> n >> k;
umap<int, int> mp;
for (int i = 0;i <= 1e6;i++) g[i] = 0;
for (int i = 1;i <= 1e6;i++) {
int t = qp((qp(2, i) - 1), k);
for (int j = 1;j * i <= 1e6;j++) {
g[i * j] = (g[i * j] + mu[j] * t % MOD) % MOD;
}
}
for (int i = 1;i <= 1e6;i++) {
g[i] = (g[i - 1] + g[i]) % MOD;
}
int res = 0;
for (int l = 1, r;l <= n;l = r + 1) {
r = n;
if (n / l) r = min(r, n / (n / l));
res = add(res, (n / l) * (n / l) % MOD * (Du_sum_g(r, g, mp, k) - Du_sum_g(l - 1, g, mp, k)) % MOD);
}
cout << norm(res) << endl;
}法2.
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[\gcd(2^i-1,2^j-1)]^k=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}(2^{gcd(i,j)}-1)^k$
$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{d=1}^{n}(2^{d}-1)^k[gcd(i,j)=d]=\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{d=1}^{n}(2^{d}-1)^k[gcd(i,j)=1]$
$=\sum\limits_{d=1}^{n}(2^{d}-1)^k(\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}2\varphi(i)-1)$
考虑对后面的$\lfloor\frac{n}{d}\rfloor$整除分块，杜教筛处理后面的$\varphi$的前缀和，前面的也可以通过上面所说的变成等比数列，每次$klogn$复杂度的查询。
根据$\varphi*1=id$做杜教筛
$\sum\limits_{i=1}^{n}id(i)=\sum\limits_{i=1}^{n}\varphi(i)\sum\limits_{j=1}^{\lfloor \frac{n}{i} \rfloor}1=\sum\limits_{i=1}^{n}1\sum\limits_{j=1}^{\lfloor \frac{n}{i} \rfloor}\varphi(j)=\sum\limits_{i=1}^{n}\varphi(i)+\sum\limits_{i=2}^{n}\sum\limits_{j=1}^{\lfloor \frac{n}{i} \rfloor}\varphi(j)$
即$\sum\limits_{i=1}^{n}\varphi(i)=\sum\limits_{i=1}^{n}id(i)-\sum\limits_{i=2}^{n}\sum\limits_{j=1}^{\lfloor \frac{n}{i} \rfloor}\varphi(j)=\frac{n(n+1)}{2}-\sum\limits_{i=2}^{n}\sum\limits_{j=1}^{\lfloor \frac{n}{i} \rfloor}\varphi(j)$
查看代码
int g[N];
int sum_f(int n, int k) {
int res = 0;
for (int j = 0;j <= k;j++) {
int sgn = (k - j) % 2 == 1 ? MOD - 1 : 1;
int q = (1 << j);
if (j == 0) {
res = (res + n * sgn % MOD) % MOD;
}
else {
res = add(res, sgn * comb(k, j) % MOD * (q * (qp(q, n) - 1) % MOD * inv(q - 1) % MOD) % MOD);
}
}
return norm(res);
}
int Du_sum_phi(int n, int pre_f[], umap<int, int>& mp) {
if (n < N) return pre_f[n];
if (mp.count(n)) return mp[n];
int res = n * (n + 1) / 2;
for (int l = 2, r;l <= n;l = r + 1) {
r = n;
if (n / l) r = min(r, n / (n / l));
res = (res - Du_sum_phi(n / l, pre_f, mp) * (r - l + 1) % MOD) % MOD;
}
return mp[n] = res;
}
umap<int, int> mp_sum_phi;
void Solve(int TIME) {
int n, k;cin >> n >> k;
int res = 0;
for (int l = 1, r;l <= n;l = r + 1) {
r = n;
if (n / l) r = min(r, n / (n / l));
res = add(res, (sum_f(r, k) - sum_f(l - 1, k)) % MOD * (2 * Du_sum_phi(n / l, pre_phi, mp_sum_phi) % MOD - 1) % MOD);
}
cout << norm(res) << endl;
}
9.[2021CCPC威海热身赛] Number Theory
知识点
打表
题意
$\sum\limits_{k=1}^{n} \sum\limits_{i|k} \sum\limits_{j|i}\lambda(i)\lambda(j)$，其中$\lambda(x)=(-1)^{\sum\limits_{i=1}^{k}c_i}$，$x=\prod\limits_{i=1}^{k}p_i^{c_i}$ 。$1\le n\le 10^{12}$
题解
打表易知$\sum\limits_{d|n}\lambda(d)=[n=a^2,a\in N^+]$
$\sum\limits_{k=1}^{n} \sum\limits_{i|k} \sum\limits_{j|i}\lambda(i)\lambda(j)=\sum\limits_{k=1}^{n} \sum\limits_{i|k}\lambda(i) [i=a^2,a\in N^{+}]$
设$f(i)=\lambda(i) [i=a^2,a\in N^{+}]=[i=a^2,a\in N^{+}]$，那么就有$\sum\limits_{k=1}^{n} \sum\limits_{i|k}f(i)=\sum\limits_{i=1}^{n}f(i)\lfloor \frac{n}{i} \rfloor$
这样枚举平方数即可，复杂度$O(\sqrt{n})$
查看代码
void Solve(int TIME) {
int n;cin >> n;int res = 0;
for (int i = 1;i * i <= n;i++) {
res = (res + n / (i * i) % MOD) % MOD;
}
cout << res << endl;
}
10.2019ICPC西安区域赛 D - Easy Problem
知识点
欧拉降幂，莫比乌斯反演
题解
$\sum\limits_{a_1=1}^{m}\sum\limits_{a_2=1}^{m}...\sum\limits_{a_n=1}^{m}(a_1a_2...a_n)^k[gcd(a_1,a_2,...,a_n)=d]$
$=d^{nk}\sum\limits_{a_1=1}^{\lfloor \frac{m}{d }\rfloor}\sum\limits_{a_2=1}^{\lfloor \frac{m}{d }\rfloor}...\sum\limits_{a_n=1}^{\lfloor \frac{m}{d }\rfloor}(a_1a_2...a_n)^k[gcd(da_1,da_2,...,da_n)=d]$
$=d^{nk}\sum\limits_{a_1=1}^{\lfloor \frac{m}{d }\rfloor}\sum\limits_{a_2=1}^{\lfloor \frac{m}{d }\rfloor}...\sum\limits_{a_n=1}^{\lfloor \frac{m}{d }\rfloor}(a_1a_2...a_n)^k[gcd(a_1,a_2,...,a_n)=1]$
$=d^{nk}\sum\limits_{a_1=1}^{\lfloor \frac{m}{d }\rfloor}\sum\limits_{a_2=1}^{\lfloor \frac{m}{d }\rfloor}...\sum\limits_{a_n=1}^{\lfloor \frac{m}{d }\rfloor}(a_1a_2...a_n)^k\sum\limits_{D|gcd(a_1,a_2,...,a_n)}\mu(D)$
$=d^{nk}\sum\limits_{a_1=1}^{\lfloor \frac{m}{d }\rfloor}\sum\limits_{a_2=1}^{\lfloor \frac{m}{d }\rfloor}...\sum\limits_{a_n=1}^{\lfloor \frac{m}{d }\rfloor}(a_1a_2...a_n)^k\sum\limits_{D=1}^{\lfloor \frac{m}{d }\rfloor}\mu(D)[D|gcd(a_1,a_2,...,a_n)]$
$=d^{nk}D^{nk}\sum\limits_{a_1=1}^{\lfloor \frac{m}{dD }\rfloor}\sum\limits_{a_2=1}^{\lfloor \frac{m}{dD }\rfloor}...\sum\limits_{a_n=1}^{\lfloor \frac{m}{dD }\rfloor}(a_1a_2...a_n)^k\sum\limits_{D=1}^{\lfloor \frac{m}{d }\rfloor}\mu(D)[D|gcd(Da_1,Da_2,...,Da_n)]$
$=(dD)^{nk}\sum\limits_{a_1=1}^{\lfloor \frac{m}{dD }\rfloor}\sum\limits_{a_2=1}^{\lfloor \frac{m}{d D}\rfloor}...\sum\limits_{a_n=1}^{\lfloor \frac{m}{dD }\rfloor}(a_1a_2...a_n)^k\sum\limits_{D=1}^{\lfloor \frac{m}{d }\rfloor}\mu(D)$
$=\sum\limits_{D=1}^{\lfloor \frac{m}{d }\rfloor}\mu(D)(dD)^{nk}\sum\limits_{a_1=1}^{\lfloor \frac{m}{d D}\rfloor}a_1^k\sum\limits_{a_2=1}^{\lfloor \frac{m}{d D}\rfloor}a_2^k...\sum\limits_{a_n=1}^{\lfloor \frac{m}{d D}\rfloor}a_n^k$
$=d^{nk}\sum\limits_{D=1}^{\lfloor \frac{m}{d }\rfloor}\mu(D)D^{nk}(\sum\limits_{i=1}^{\lfloor \frac{m}{d D}\rfloor}i^k)^n$
然后用欧拉降幂即可。
可以线性筛预处理$i^k$，与$i^k$的前缀和，然后$O(mlog(mod))$扫过去。
查看代码
int k;
void init(int mod = MOD) {
np[0] = np[1] = 1;mu[1] = 1;idk[1] = 1;
for (int i = 2;i < N;i++) {
if (!np[i]) {
p.push_back(i);
mu[i] = -1;
idk[i] = qp(i, k);
}
for (auto j : p) {
if (i * j < N) {
np[i * j] = 1;
if (i % j == 0) {
idk[i * j] = idk[i] * idk[j] % MOD;
break;
}
else {
mu[i * j] = -mu[i];
idk[i * j] = idk[i] * idk[j] % MOD;
}
}
else {
break;
}
}
}
for (int i = 1;i < N;i++) {
sum_idk[i] = (idk[i] + sum_idk[i - 1]) % MOD;
}
}
int phi_mod = getphi(MOD);
void Solve(int TIME) {
string n;cin >> n;
int m, d;cin >> m >> d >> k;
init();
int res = qp(d, k);
int nn = exp_euler(n, MOD);
res = qp(res, nn);
int t = 0;
for (int i = 1;i <= m / d;i++) {
t = (t + mu[i] * qp(idk[i], nn) % MOD * qp(sum_idk[m / d / i], nn) % MOD) % MOD;
}
cout << norm(res * t) << endl;
}11.CF1900D - Small GCD
知识点
欧拉反演，gcd的调和级数容斥，莫比乌斯容斥
题意
$\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}gcd(a_i,a_j) \cdot(n-j)$
题解
考虑右端点递增，然后思考如何去维护左侧所有区间的信息。
$f(t)$表示$gcd(a_i,a_j)=t$的对数
考虑容斥：$g(t)$表示$t|gcd(a_i,a_j)$的对数
$g(t)=f(t)+f(2t)+f(3t)+...$
那么$f(t)=g(t)-f(2t)-f(3t)-...$
而$g(t)$可以通过预处理得到，对每个数的因子进行判断记录贡献即可。
答案是$\sum\limits_{i=1}^{N}i\cdot f(i)$
查看代码
vvc<int> factor(N + 10);
void Solve(int TIME) {
int T;cin >> T;
for (int i = 1;i <= N;i++) {
for (int j = 1;j * i <= N;j++) {
factor[i * j].push_back(i);
}
}
while (T--) {
int n;cin >> n;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
sort(a.begin() + 1, a.end());
vi cnt(N + 10);
vi f(N + 10);
for (int i = 1;i <= n;i++) {
for (auto x : factor[a[i]]) {
f[x] += cnt[x] * (n - i);
cnt[x] += 1;
}
}
vi F(N + 10);
int res = 0;
for (int i = N;i >= 1;i--) {
F[i] = f[i];
for (int j = i + i;j <= N;j += i) {
F[i] -= F[j];
}
res += F[i] * i;
}
cout << res << endl;
}
}法2.
利用欧拉反演：$gcd(a_i,a_j)=\sum\limits_{d|gcd(a_i,a_j)}\varphi(d)$
$\sum\limits_{i=1}^{j-1}\sum\limits_{j=1}^{n}gcd(a_i,a_j) \cdot (n-j)=\sum\limits_{i=1}^{j-1}[d|a_i]\sum\limits_{j=1}^{n}[d|a_j](n-j)\sum\limits_{d=1}^{n}\varphi(d) $
$=\sum\limits_{d=1}^{N}\sum\limits_{j=1}^{n}[d|a_j](n-j)\sum\limits_{i=1}^{j-1}[d|a_i]$
枚举右端点，枚举当前$a_j$的约数记录贡献即可。
查看代码
void Solve(int TIME) {
int T;cin >> T;
vvc<int> factor(N + 10);
for (int i = 1;i <= N;i++) {
for (int j = 1;j * i <= N;j++) {
factor[i * j].push_back(i);
}
}
while (T--) {
int n;cin >> n;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
sort(a.begin() + 1, a.end());
int res = 0;
vi cnt(N + 10);
for (int i = 1;i <= n;i++) {
for (auto x : factor[a[i]]) {
res += cnt[x] * (n - i) * phi[x];
cnt[x] += 1;
}
}
cout << res << endl;
}
}法3.利用莫比乌斯反演的倍数形式
$\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}\sum\limits_{x=1}^{N}x[gcd(a_i,a_j)=x]\cdot(n-j)$
$=\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}\sum\limits_{x=1}^{N}x\sum\limits_{d=x,2x,...}[d|gcd(a_i,a_j)]\mu(\frac{d}{x})\cdot(n-j)$
$=\sum\limits_{x=1}^{N}x\sum\limits_{d=x,2x,...}\sum\limits_{i=1}^{n}[d|a_i]\sum\limits_{j=i+1}^{n}[d|a_j]\cdot(n-j)\mu(\frac{d}{x})$
令$f(d)=\sum\limits_{i=1}^{n}[d|a_i]\sum\limits_{j=i+1}^{n}[d|a_j]\cdot(n-j)$
原式$=\sum\limits_{x=1}^{N}x\sum\limits_{d=x,2x,...}f(d)\mu(\frac{d}{x})$
查看代码
vvc<int> factor(N + 10);
void Solve(int TIME) {
int T;cin >> T;
for (int i = 1;i <= N;i++) {
for (int j = 1;j * i <= N;j++) {
factor[i * j].push_back(i);
}
}
while (T--) {
int n;cin >> n;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
sort(a.begin() + 1, a.end());
int res = 0;
vi cnt(N + 10);
vi f(N + 10);
for (int i = 1;i <= n;i++) {
for (auto x : factor[a[i]]) {
f[x] += cnt[x] * (n - i);
cnt[x] += 1;
}
}
for (int i = 1;i <= N;i++) {
for (int j = 1;i * j <= N;j++) {
res += i * f[i * j] * mu[j];
}
}
cout << res << endl;
}
}
12.I-hx的数列
知识点
数论分块，莫比乌斯反演
题解
$(a,a\cdot \frac{c}{b},a\cdot{(\frac{c}{b})}^2)$，令$t=a\cdot (\frac{c}{b})^2$，需要$t\le n$
保证$gcd(b,c)=1$且$b<c$，于是就有$b^2|a$，这样才能有$t$为整数
于是令$k=\frac{a}{b^2}$，则有$k\cdot c^2 \le n$，于是有$k\le \lfloor\frac{n}{c^2}\rfloor$
$\sum\limits_{c=1}^{\lfloor\sqrt{n}\rfloor}\lfloor\frac{n}{c^2}\rfloor \sum\limits_{b=1}^{c-1}[\gcd(b,c)=1]=\sum\limits_{c=1}^{\lfloor\sqrt{n}\rfloor}\lfloor\frac{n}{c^2}\rfloor \sum\limits_{b=1}^{c-1}\sum\limits_{D|gcd(b,c)}\mu(D)$
$=\sum\limits_{Dc=1}^{\lfloor\sqrt{n}\rfloor}\lfloor\frac{n}{D^2c^2}\rfloor \sum\limits_{Db=1}^{Dc-1}\sum\limits_{D=1}^{\lfloor\sqrt{n}\rfloor}\mu(D)[D|gcd(Db,Dc)]$
$=\sum\limits_{Dc=1}^{\lfloor\sqrt{n}\rfloor}\lfloor\frac{n}{D^2c^2}\rfloor \sum\limits_{Db=1}^{Dc-1}\sum\limits_{D=1}^{\lfloor\sqrt{n}\rfloor}\mu(D)=\sum\limits_{D=1}^{\lfloor\sqrt{n}\rfloor}\mu(D)\sum\limits_{c=1}^{\lfloor\frac{\lfloor\sqrt{n}\rfloor}{D}\rfloor}\lfloor\frac{n}{D^2c^2}\rfloor \sum\limits_{b=1}^{c-1}1$
$=\sum\limits_{D=1}^{\lfloor\sqrt{n}\rfloor}\mu(D)\sum\limits_{c=1}^{\lfloor\frac{\lfloor\sqrt{n}\rfloor}{D}\rfloor}\lfloor\frac{n}{D^2c^2}\rfloor (c-1)$
复杂度$O(\sqrt1+\sqrt2+...+\sqrt{\sqrt{n}})$
大抵是不会很大的。
查看代码
void Solve(int TIME) {
int n;cin >> n;
int res = 0;
int sqn = sqrt(n);
for (int i = 1;i <= sqn - 1;i++) {
int t = 0;
for (int l = 1, r;l <= sqn / i;l = r + 1) {
int m = n / i / i;
r = sqn / i;
if (m / (l * l)) r = min((int)sqrt(m / (m / (l * l))), r);
t += m / (l * l) * (l + r - 2) * (r - l + 1) / 2;
}
res += t * mu[i];
}
cout << res << endl;
}
13.牢大与GCD
知识点
莫比乌斯容斥
题意
给出长度为$n$的数组$a$和$num$，求下面的式子。
$\sum\limits_{x=1}^{n}\sum\limits_{k_1=1}^{n}\sum\limits_{k_2=1}^{n}...\sum\limits_{k_{num[x]}=1}^{n}[\gcd\limits_{i=1}^{num[x]}(a_{k_i})=x]$
题解
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{k=1}^{n}[gcd(a_i,a_j,a_k)=x]$
$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{k=1}^{n}\sum\limits_{d=x,2x,...}[d|gcd(a_i,a_j,a_k)]\mu(\frac{d}{x})$
$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\sum\limits_{k=1}^{n}\sum\limits_{d=x,2x,...}[d|a_i][d|a_j][d|a_k]\mu(\frac{d}{x})$
$=\sum\limits_{d=x,2x,...}\mu(\frac{d}{x})\sum\limits_{i=1}^{n}[d|a_i]\sum\limits_{j=1}^{n}[d|a_j]\sum\limits_{k=1}^{n}[d|a_k]$
$=\sum\limits_{d=x,2x,...}\mu(\frac{d}{x})(\sum\limits_{i=1}^{n}[d|a_i])^3$
可以看出，当需要gcd的个数扩大至$X$，变的仅仅是后面的幂次。
于是令$f(d)=\sum\limits_{i=1}^{n}[d|a_i]$
结果为$\sum\limits_{x=1}^{n}\sum\limits_{d=x,2x,...}\mu(\frac{d}{x})f(d)^{pw[x]}$
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
int M = 1e5;
vi cnt(M + 1);for (int i = 1;i <= n;i++) cnt[a[i]] += 1;
//正解
vi pw(n + 1);for (int i = 1;i <= n;i++) cin >> pw[i];
vi f(M + 1);
for (int i = 1;i <= M;i++) {
for (int j = 1;j * i <= M;j++) {
f[i] += cnt[i * j];
}
}
int res = 0;
for (int i = 1;i <= n;i++) {
for (int j = 1;j * i <= M;j++) {
res += mu[j] * qp(f[i * j], pw[i]);
}
}
cout << res << endl;
//暴力对拍
// int res = 0;
// for (int i = 1;i <= n;i++) {
// map<int, int> mp;mp[0] = 1;
// for (int j = 1;j <= pw[i];j++) {
// map<int,int> nmp;
// for (int k = 1;k <= n;k++) {
// for (auto [o1, o2] : mp) {
// nmp[gcd(o1, a[k])] += o2;
// }
// }
// mp = nmp;
// }
// res += mp[i];
// }
// cout << res << endl;
}
14.G-小苯的逆序对_牛客小白月赛87 (nowcoder.com)
知识点
莫比乌斯容斥及其dp方法
题解
令$f(x)$表示$gcd(a_i,a_j)=x$并且$i和j$是逆序对的对数。
令$g(x)=f(x)+f(2x)+f(3x)...$，那么$f(x)=g(x)-(f(2x)+f(3x)+...)$ 。这里可以用dp做。
$f(x)=g(x)\mu(1)+g(2x)\mu(2)+...+g(kx)\mu(k)+...$ ，这里转换成莫比乌斯容斥。
$g(x)$是好算的，树状数组维护即可，做完了。
dp:
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi a(n + 1);
vi pos(n + 1);
for (int i = 1;i <= n;i++) cin >> a[i], pos[a[i]] = i;
BIT tr(n + 1);
vi f(n + 1);
for (int i = n;i >= 1;i--) {
for (int j = i;j <= n;j += i) {
f[i] += tr.query(pos[j], n);
tr.add(pos[j], 1);
}
for (int j = i;j <= n;j += i) {
tr.add(pos[j], -1);
}
for (int j = 2 * i;j <= n;j += i) {
f[i] -= f[j];
}
}
cout << f[1] << endl;
}莫比乌斯容斥：
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi a(n + 1);
vi pos(n + 1);
for (int i = 1;i <= n;i++) cin >> a[i], pos[a[i]] = i;
BIT tr(n + 1);
int res = 0;
for (int i = 1;i <= n;i++) {
int t = 0;
for (int j = i;j <= n;j += i) {
t += tr.query(pos[j], n);
tr.add(pos[j], 1);
}
for (int j = i;j <= n;j += i) {
tr.add(pos[j], -1);
}
res += t * mu[i];
}
cout << res << endl;
}
15.P2398 GCD SUM
知识点
莫比乌斯反演
题解
莫比乌斯反演
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}gcd(i,j)=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[gcd(i,j)=d]=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(di,dj)=d]$
$=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(i,j)=1]=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{D|gcd(i,j)}\mu(D)$
$=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{D=1}^{\lfloor\frac{n}{d}\rfloor}\mu(D)[D|gcd(i,j)]=\sum\limits_{d=1}^{n}d\sum\limits_{D=1}^{\lfloor\frac{n}{d}\rfloor}\mu(D)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}[D|i]\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}[D|j]$
$=\sum\limits_{d=1}^{n}d\sum\limits_{D=1}^{\lfloor\frac{n}{d}\rfloor}\mu(D)(\lfloor\frac{n}{Dd}\rfloor)^2$
查看代码
void Solve(int TIME) {
int n;cin >> n;
int res = 0;
for (int i = 1;i <= n;i++) {
int m = n / i;
for (int l = 1, r;l <= m;l = r + 1) {
r = m;
if (m / l) r = min(m / (m / l), r);
res += (pre_mu[r] - pre_mu[l - 1]) * i * (m / l) * (m / l);
}
}
cout << res << endl;
}
莫比乌斯容斥
$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}gcd(i,j)=\sum\limits_{d=1}^{n}d\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[gcd(i,j)=d]$
令$g(d)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[gcd(i,j)=d]+\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}[gcd(i,j)=2d]+...$
得到$\sum\limits_{d=1}^{n}d\sum\limits_{k=1}^{n}\mu(k)g(kd)$
$g(x)$很好算，即每个$x$的倍数都可以选，选两个可以重复。即$(\lfloor\frac{n}{x} \rfloor)^2$
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi g(n + 1);
for (int i = 1;i <= n;i++) {
g[i] = (n / i) * (n / i);
}
int res = 0;
for (int i = 1;i <= n;i++) {
for (int j = 1;j * i <= n;j++) {
res += i * mu[j] * g[j * i];
}
}
cout << res << endl;
}
数论与杂项
1. 51nod1236 - 序列求和 V3
知识点
特征方程，二次剩余
题解
$F[n]=F[n-1]+F[n-2]$ ，由特征方程 $q^2=q+1$，得到特征根$q_1=\frac{1+\sqrt5}{2}$ ,$q_2=\frac{1-\sqrt5}{2}$
则有$F[0]=k_1q_1^0+k_2q_2^0=0$ ，$F[1]=k_1q_1^1+k_2q_2^1=1$，
解得$k_1=\frac{1}{\sqrt5}$ , $k_2=-\frac{1}{\sqrt5}$ , 得到 $F[n]=\frac{1}{\sqrt5}(\frac{1+\sqrt5}{2})^n-\frac{1}{\sqrt5}(\frac{1-\sqrt5}{2})^n=\frac{1}{\sqrt5}[(\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n]$
$F[n]^k=(\frac{1}{\sqrt5})^k[(\frac{1+\sqrt5}{2})^n-(\frac{1-\sqrt5}{2})^n]^k =(\frac{1}{\sqrt5})^k\sum\limits_{i=0}^{k}C_{k}^{i}(\frac{1+\sqrt5}{2})^{ni}(\frac{1-\sqrt5}{2})^{n(k-i)}(-1)^{k-i}$
$S(n,k)=\sum\limits_{i=1}^{n}F[i]^k=\sum\limits_{i=1}^{n}(\frac{1}{\sqrt5})^k\sum\limits_{j=0}^{k}C_{k}^{j}(\frac{1+\sqrt5}{2})^{ij}(\frac{1-\sqrt5}{2})^{i(k-j)}(-1)^{k-j}$
$=(\frac{1}{\sqrt5})^k\sum\limits_{i=1}^{n}\sum\limits_{j=0}^{k}C_{k}^{j}[(\frac{1+\sqrt5}{2})^{j}(\frac{1-\sqrt5}{2})^{(k-j)}]^i(-1)^{k-j}$
$=(\frac{1}{\sqrt5})^k\sum\limits_{j=0}^{k}(-1)^{k-j}C_{k}^{j}\sum\limits_{i=1}^{n}[(\frac{1+\sqrt5}{2})^{j}(\frac{1-\sqrt5}{2})^{(k-j)}]^i$
后面的循环是个等比数列，快速幂处理即可，总时间复杂度$O(klogn)$
$\sqrt5$在模意义下的值用二次剩余得到 ，可以直接打表找或使用Cipolla算法。在本题中可以得到 383008016是5的一个二次剩余。
查看代码
void Solve(int TIME) {
int n, k;cin >> n >> k;
int res = 0;
int a = (1 + sqrt5) % MOD * inv(2) % MOD;
int b = (MOD + 1 - sqrt5) % MOD * inv(2) % MOD;
for (int j = 0;j <= k;j++) {
int t = qp(a, j) * qp(b, k - j) % MOD;
int tt;
if (t != 1) tt = t * (qp(t, n) - 1) % MOD * inv(t - 1) % MOD;
else tt = n % MOD;
if ((k - j) & 1) {
res = sub(res, mul(comb(k, j), tt));
}
else {
res = add(res, mul(comb(k, j), tt));
}
}
res = mul(res, qp(inv(sqrt5), k));
cout << norm(res) << endl;
}
2.CF1891D - Suspicious logarithms
知识点
数论分块
题解
$2^{f(x)}\le x$，则有$f(x)=\lfloor log_2x \rfloor$
$f(x)^{g(x)}\le x$，则有$g(x)\cdot log_2f(x)\le log_2x$，即$g(x)\cdot log_2{\lfloor log_2x \rfloor}\le log_2x$，即$g(x)=\lfloor \frac{log_2x}{log_2{\lfloor log_2x \rfloor}} \rfloor$，即求$\sum\limits_{i=l}^{r}g(i)$
对于相等的$\lfloor log_2x \rfloor$，当$x$增大，$g(x)$也增大。而$\lfloor log_2x \rfloor$的范围比较小，考虑枚举。之后在该范围，将其划分为多个不同的连续相同数字子段。
注意：$log_2(x)$这类函数默认使用$double$类型，使用前要强制转换为$long\; double$，否则会爆。
查看代码
int f(int x) {
return 63 - __builtin_clzll(x);
}
int g(int x) {
return (ll)(log2(ld(x)) / log2((ld)f(x)));
}
struct qwq {
int l, r, v;
};
void Solve(int TIME) {
vc<qwq> a;int L = 4, R = 1E18;
for (int l = L, r;l <= R;l = r + 1) {
int fx = f(l);int lim = min((1ll << (fx + 1)) - 1, R);
int lo = l + 1, hi = lim;
int gx = g(l);
while (lo <= hi) {
int mid = (lo + hi) >> 1;
if (g(mid) <= gx) {
lo = mid + 1;
}
else {
hi = mid - 1;
}
}
r = lo - 1;
a.push_back({ l,r,gx });
}
int q;cin >> q;
while (q--) {
cin >> L >> R;
int res = 0;
auto pt = partition_point(a.begin(), a.end(), [&](qwq A) {return A.r < L;}) - a.begin();
for (int i = pt;i < a.size();i++) {
auto [l, r, v] = a[i];
if (l > R) break;
else {
l = max(l, L), r = min(r, R);
res = (res + (r - l + 1) * v % MOD) % MOD;
}
}
cout << res << endl;
}
}
3.CFgym104076 E - Identical Parity
知识点
exgcd
题解
发现最后$1,1+k,1+2k...$ ， $2,2+k,2+2k...$, ... $k-1,k-1+k,k-1+2k,....$这些都是同奇偶性。
也就是分成$k$组，其中$n \;mod \;k$组是$\lfloor\frac{n}{k}\rfloor+1$个元素，剩下$k-n \;mod\;k$组是$\lfloor\frac{n}{k}\rfloor$个元素。令$a=\lfloor\frac{n}{k}\rfloor$
然后要使得$ax+(a+1)y=n/2$有解$x,y$使得$0\le x\le (n \;mod\;k)$ 并且$0\le y \le (k-n \;mod\;k) $
二分求出每个边界的可行范围看有无交集即可。
查看代码
void Solve(int TIME) {
int n, k;cin >> n >> k;
int a = n / k;//人数
int t = n % k;//人数为a+1的组数t，人数为a的组数k-t
int x, y;
int d = exgcd(a, a + 1, x, y);
int m = n / 2;
if (m % d) return NO, void();
int base = m / d;
x *= base, y *= base;
int norm_1 = a + 1;
int norm_2 = a;
int l, r;
l = -1e9, r = 1e9;
while (l <= r) {
int mid = l + r >> 1;
if (x + mid * (a + 1) <= k - t) l = mid + 1;
else r = mid - 1;
}
int r_a = l - 1;
l = -1e9, r = 1e9;
while (l <= r) {
int mid = l + r >> 1;
if (x + mid * (a + 1) >= 0) r = mid - 1;
else l = mid + 1;
}
int l_a = r + 1;
l = -1e9, r = 1e9;
while (l <= r) {
int mid = l + r >> 1;
if (y - mid * (a) <= t) r = mid - 1;
else l = mid + 1;
}
int l_b = r + 1;
l = -1e9, r = 1e9;
while (l <= r) {
int mid = l + r >> 1;
if (y - mid * (a) >= 0) l = mid + 1;
else r = mid - 1;
}
int r_b = l - 1;
if (min(r_a, r_b) >= max(l_a, l_b)) YES;
else NO;
}4.2023ICPC区域赛南京站C
知识点
异或的性质
题意
求满足 $ i \text{^} (P-1)\equiv 1 (mod\; P)$ 并且$0\le i\le m$的$i$的个数
范围$1\le m \le 10^{18} $，$1\le P \le 10^{18}$
题解
即 $i\text{^} (P-1)= kP+1$ ，$(k\ge 0)$ ，即$i=(kP+1) \text{^} (P-1)\le m$的数量
我们知道异或的性质$a-b \le $ $a \text{^} b$ $\le a+b$
所以$(k-1)P+2\le(kP+1) \text{^} (P-1) \le(k+1)P$
当$k\ge \lceil\frac{m}{p}\rceil+1$ ，$(kP+1) \text{^} (P-1) \ge(k-1)P+2 > m$
当$0\le k\le \lfloor\frac{m}{p}\rfloor-1$ ，$(kP+1) \text{^} (P-1) \le(k+1)P \le m$
所以只需检查$\lfloor\frac{m}{p}\rfloor \le k\le \lceil\frac{m}{p} \rceil$范围的$k$是否满足$(kP+1)\text{^} (P-1)\le m$即可
查看代码
void Solve(int TIME) {
int P, m;cin >> P >> m;
int res = m / P;
for (int i = m / P;i <= (m + P - 1) / P;i++) {
int t = (i * P + 1) ^ (P - 1);
if (t <= m) res += 1;
}
cout << res << endl;
}好像不需要卡这么死，知道大致上界后往后枚举个10位就行了。
查看代码
void Solve(int TIME) {
const __int128 ONE = 1;
int p, m;cin >> p >> m;
int res = m / p;
int n = m / p - 1;
for (int i = n + 1;i <= n + 10;i++) {
__int128 t = ONE * (ONE * i * p + 1) ^ (ONE * p - 1);
if (t <= m) res += 1;
}
cout << res << endl;
}
解2
根据上面的式子，将$m$拆位枚举，枚举公共前缀，在$1$的位置替换为$0$，然后后面的可以为任意。之后用类似前缀和的方式算一下答案。
查看代码
void Solve(int TIME) {
int p, m;cin >> p >> m;
int mask = 0;
int res = 0;
auto calc = [&](int x) {//[0,x]中满足kp+1=x的个数
x--;
if (x < 0) return 0ll;
return x / p + 1;
};
auto calc2 = [&](int l, int r) {//[l,r]中满足
return calc(r) - calc(l - 1);
};
for (int i = 60;i >= 0;i--) {
if (~m >> i & 1) continue;
int t = mask ^ (p - 1);
int l = t >> i << i, r = l + (1ll << i) - 1;
res += calc2(l, r);
mask += 1ll << i;
}
res += calc2(mask ^ (p - 1), mask ^ (p - 1));
cout << res << endl;
}
5.2019ICPC西安站 F - Function!
知识点
数论分块
题意
$\sum\limits_{a=2}^{n}(a\sum\limits_{b=a}^{n}\lfloor log_ab \rfloor)$
题解
$\sum\limits_{a=2}^{n}(a\sum\limits_{b=a}^{n}\lfloor log_ab \rfloor \lceil log_ba\rceil)=\sum\limits_{a=2}^{n}(a\sum\limits_{b=a}^{n}\lfloor log_ab \rfloor)$
发现当$a\cdot a > n$，$\lfloor log_ab \rfloor =1$ ，此时有$\sum\limits_{a=\lfloor\sqrt{n}\rfloor+1}^{n}a(n-a+1)=\sum\limits_{a=\lfloor\sqrt{n}\rfloor+1}^{n}\left(a(n+1)-a^2\right)$
于是在$a\cdot a \le n$时，分块计算，$\sum\limits_{a=2}^{\lfloor \sqrt{n} \rfloor}\sum\limits_{b=a}^{n}\lfloor log_ab \rfloor$
复杂度$O(log1+log2+...+log\sqrt{n})=O(log{\frac{(1+\sqrt{n})\sqrt{n}}{2}})=O(logn)$
查看代码
int inv6 = inv(6);
int inv2 = inv(2);
int calc_n(int x) {
x %= MOD;
return (1 + x) % MOD * x % MOD * inv2 % MOD;
}
int calc_n_2(int x) {
x %= MOD;
return (2 * x + 1) % MOD * (x + 1) % MOD * x % MOD * inv6 % MOD;
}
void Solve(int TIME) {
int n;cin >> n;
int res = 0;
int a = 2;
for (a = 2;a * a <= n;a++) {
for (int l = a, r, k = 1;l <= n;l = r + 1, k++) {
r = min(l * a - 1, n);
res = (res + (r - l + 1) % MOD * k % MOD * a % MOD) % MOD;
}
}
res = (res + (n + 1) % MOD * (calc_n(n) - calc_n(a - 1)) % MOD) % MOD;
res = (res - (calc_n_2(n) - calc_n_2(a - 1)) % MOD) % MOD;
res = norm(res);
cout << res << endl;
}
6.P5110 块速递推
知识点
二次剩余，特征方程，光速幂
题解
$a_n=233a_{n-1}+666a_{n-2}$ ，特征方程$q^2=233q+666$
得到$\begin{cases} q_1=\frac{233+\sqrt{56953}}{2} \\q_2= \frac{233-\sqrt{56953}}{2} \end{cases}$
通项$a_n=k_1\cdot q_1^n+k_2\cdot q_2^n$，带入$a_0=k_1+k_2=0$ ，$a_1=k_1\cdot q_1+k_2\cdot q_2=1$
得到$k_1=\frac{1}{q_1-q_2}=\frac{1}{\sqrt{56953}}$，$k_2=-\frac{1}{\sqrt{56953}}$。则有$a_n=\frac {(\frac{233+\sqrt{56953}}{2})^n-(\frac{233-\sqrt{56953}}{2})^n}{\sqrt{56953}}$
发现$\sqrt{56953}$在$10^9+7$下的二次剩余为$188305837$，那么直接带入进去即可。
$a_n=\frac {(\frac{233+188305837}{2})^n-(\frac{233-188305837}{2})^n}{\sqrt{56953}}=\frac {94153035^n-(-94152802)^n}{188305837}=\frac {94153035^n-905847205^n}{188305837}(mod \;10^9+7) $
$=233230706\cdot(94153035^n-905847205^n)$
然后光速幂$O(\sqrt{MOD})$预处理，$O(1)$查询即可。总复杂度$O(T+\sqrt{MOD})$
查看代码
struct LP {//光速幂,处理同底数同模数的幂
int base1[N], basesqrt[N];
int Block_len;
int Phi;
ll maxn = 1e10;//模数的最大值
LP(int x) { init(x); }
void init(int x) {//初始化底数为x
Phi = MOD - 1;
Block_len = sqrt(maxn) + 1;
base1[0] = 1;for (int i = 1;i <= Block_len;i++) base1[i] = 1ll * base1[i - 1] * x % MOD;
basesqrt[0] = 1;for (int i = 1;i <= Block_len;i++) basesqrt[i] = 1ll * basesqrt[i - 1] * base1[Block_len] % MOD;
}
int qp(ull x) {
x %= Phi;
return 1ll * basesqrt[x / Block_len] * base1[x % Block_len] % MOD;
}
};
namespace Mker
{
unsigned long long SA, SB, SC;
void init() { cin >> SA >> SB >> SC; }
unsigned long long rand()
{
SA ^= SA << 32, SA ^= SA >> 13, SA ^= SA << 1;
unsigned long long t = SA;
SA = SB, SB = SC, SC ^= t ^ SA;return SC;
}
}
void Solve(int TIME) {
int T;cin >> T;
Mker::init();
LP a(94153035);
LP b(905847205);
int res = 0;
while (T--) {
int n = Mker::rand();
res ^= 233230706 * (a.qp(n) - b.qp(n) + MOD) % MOD;
}
cout << res << endl;
}
7.abc331D - Tile Pattern
知识点
容斥
题解
考虑对一个右下角$(x,y)$计算它到$(1,1)$的$1$的个数。
发现可以分为四个部分来求：
1.左上角被$(x/n)\cdot (y/n)$个$n\cdot n$覆盖的正方形区域；
2.左下角被$(y/n)$个$(x \mod n )\cdot n$覆盖的矩形区域
3.右上角被$(x/n)$个$x\cdot (y\mod n)$覆盖的矩形区域
4.右下角被$1$个$(x\mod n)\cdot (y\mod n)$覆盖的矩形区域
画图大点的图手动模拟一下可以容易得出这个结论，我用了$8\cdot 7$的矩形，$n=3$
然后容斥一下，用类似二位前缀和的方式，得到答案$res=calc(x_2,y_2)-calc(x_2,y_1-1)-calc(x_1-1,y_2)+calc(x_1-1,y_1-1)$
查看代码
void Solve(int TIME) {
int n, q;cin >> n >> q;
vvc<int> g(n + 1, vc<int>(n + 1));
for (int i = 1;i <= n;i++) {
for (int j = 1;j <= n;j++) {
char ch;cin >> ch;
if (ch == 'B') g[i][j] = 1;
}
}
vvc<int> sumG(n + 1, vc<int>(n + 1));
for (int i = 1;i <= n;i++) {
for (int j = 1;j <= n;j++) {
sumG[i][j] = sumG[i][j - 1] + sumG[i - 1][j] - sumG[i - 1][j - 1] + g[i][j];
}
}
auto calc = [&](int x, int y) {
int res = 0;
res += (x / n) * (y / n) * sumG[n][n];
res += (y / n) * sumG[x % n][n];
res += (x / n) * sumG[n][y % n];
res += sumG[x % n][y % n];
return res;
};
while (q--) {
int x1, y1, x2, y2;cin >> x1 >> y1 >> x2 >> y2;
x1 += 1, y1 += 1, x2 += 1, y2 += 1;
cout << calc(x2, y2) - calc(x2, y1 - 1) - calc(x1 - 1, y2) + calc(x1 - 1, y1 - 1) << endl;
}
}
8.青蛙的约会
知识点
exgcd
题解
$x+am=y+an+kL$
$a(m-n)-kL=y-x$
查看代码
void Solve(int TIME) {
int x, y, m, n, l;cin >> x >> y >> m >> n >> l;
int a, k;
int d = exgcd(m - n, -l, a, k);
y -= x;
if (y % d) return cout << "Impossible" << endl, void();
int base = y / d;
a *= base, k *= base;
a = norm(a, l / gcd(l, m - n));
cout << a << endl;
}
9.C-Prime Distance
知识点
线性筛，埃筛
题解
可以发现只需要用先线性筛预处理出$\sqrt{值域}$范围内的质数，然后对指定区间埃筛即可。
复杂度$O(nloglog{10^6})$
查看代码
void Solve(int TIME) {
int l, r;cin >> l >> r;
vi v(r - l + 10);
for (auto i : p) {
if (i > r) continue;
for (int j = l / i;j <= r / i;j++) {
if (j == 1) continue;
if (i * j >= l) v[j * i - l] = 1;
}
}
if (l == 1) v[0] = 1;
vi ans;
for (int i = 0;i <= r - l;i++) {
if (v[i] == 0) ans.push_back(i + l);
}
if (ans.size() < 2) cout << "There are no adjacent primes." << endl;
else {
int mn = inf;
int res = 1;
int mx = 0;
int res2 = 1;
for (int i = 1;i < ans.size();i++) {
if (ans[i] - ans[i - 1] < mn) {
mn = ans[i] - ans[i - 1];
res = i;
}
if (ans[i] - ans[i - 1] > mx) {
mx = ans[i] - ans[i - 1];
res2 = i;
}
}
cout << ans[res - 1] << "," << ans[res] << " are closest, ";
cout << ans[res2 - 1] << "," << ans[res2] << " are most distant." << endl;
}
}
10.X-factor Chains
知识点
多重集排列
题解
按唯一分解定理，将这个数字分解为若干质数的乘积，然后发现相同质因子之间可以交换顺序，其实就是个多重集排列。
查看代码
void Solve(int TIME) {
int x;cin >> x;
int res = 0;
vi ans;
for (int i = 2;i <= x / i;i++) {
int cnt = 0;
while (x % i == 0) {
res += 1;
cnt += 1;
x /= i;
}
if (cnt) ans.push_back(cnt);
}
if (x > 1) res += 1, ans.push_back(1);
vi fact(22);fact[0] = 1;
for (int i = 1;i <= 20;i++) {
fact[i] = fact[i - 1] * i;
}
debug(fact[20]);
cout << res << " ";
res = fact[res];
for (auto i : ans) {
res /= fact[i];
}
cout << res << endl;
}11.CF1840G2 - In Search of Truth
知识点
meet in the middle
题解
G1
用类似BSGS的思想，其实也是meet in the middle的思想。或者更类似根号分治？
令$T=\lceil\sqrt{n}\rceil$，$i \in[1,T]$，$j\in[1,T]$ ，这样$T\cdot i-j$可以表示所有$[0,n-1]$的数。
想要得到答案，可以令$T\cdot i-j\equiv 0(\mod n)$，那么有$T\cdot i \equiv j (\mod n)$
先预处理出前$T$个数放入哈希表，然后再每次走$T$步看有没有和前$T$个数相同的数。
查看代码
void Solve(int TIME) {
vi vis(1e6 + 1, -1);
int now;cin >> now;
vis[now] = 0;
auto add = [&](int x) {
cout << "+ " << x << endl;
cin >> now;
};
auto answ = [&](int x) {
cout << "! " << x << endl;
};
int res = now;
for (int i = 1;i <= 1000;i++) {
add(1);
if (vis[now] != -1) {
answ(res);
return;
}
res = max(res, now);
vis[now] = i;
}
for (int i = 1;;i++) {
add(1000);
if (vis[now] != -1) {
answ((i + 1) * 1000 - vis[now]);
return;
}
}
}G2
还是看不懂啊，明年再来看吧
类似题 2022区域赛杭州站I - Guess Cycle Length
12.K-寂寞沙洲冷
知识点
01分数规划
题解
$k$个生成树，也就是一开始是$n$棵树（连通块），每次连边都会减少一个连通块，于是最后变成k个的时候退出就可以了。
然后边权就是典型的01分数规划求比率生成树。注意这里求最大生成树。
Dinkelbach方法
void Solve(int TIME) {
int n, m, k;cin >> n >> m >> k;
vc<ai4> edge;
for (int i = 1;i <= m;i++) {
int x, y, a, b;cin >> x >> y >> a >> b;
edge.push_back({ x,y,a,b });
}
const ld eps = 1e-8;
ld ans = eps;
while (1) {
DSU d(n + 1);
vc<pr<ld, ai4>> g;
for (auto [x, y, a, b] : edge) {
g.push_back({ -a + b * ans,{x,y,a,b} });
}
sort(g.begin(), g.end());
ld res = 0;int now = n;
int A = 0, B = 0;
for (auto [w, u] : g) {
if (d.find(u[0]) != d.find(u[1])) {
d.merge(u[0], u[1]);
A += u[2];
B += u[3];
now -= 1;
if (now == k) break;
}
}
ld t = 1.0 * A / B;
if (abs(t - ans) < eps) break;
ans = t;
}
cout << ans << endl;
}
二分方法
void Solve(int TIME) {
int n, m, k;cin >> n >> m >> k;
vc<ai4> edge;
for (int i = 1;i <= m;i++) {
int x, y, a, b;cin >> x >> y >> a >> b;
edge.push_back({ x,y,a,b });
}
const ld eps = 1e-8;
ld l = eps, r = 1e4 + 10;
while (r - l > eps) {
DSU d(n + 1);
ld mid = (l + r) / 2;
vc<pr<ld, ai2>> g;
for (auto [x, y, a, b] : edge) {
g.push_back({ -a + b * mid,{x,y} });
}
sort(g.begin(), g.end());
ld res = 0;int now = n;
for (auto [w, u] : g) {
if (d.find(u[0]) != d.find(u[1])) {
d.merge(u[0], u[1]);
res += -w;
now -= 1;
if (now == k) break;
}
}
if (res > eps) l = mid;
else r = mid;
}
cout << l << endl;
}二分慢了4~5倍。
13.G-scx的记忆碎片
知识点
dp，组合数
题解
dp求出用1~9组成w位，总和为k的方案数。然后组合数搞一下。
注意这里m很大，组合数需要特殊处理。
注意特判k=1的情况。
查看代码
int comb(int a, int b, int mod = MOD) {//C(a,b)=A(a,b)/b!
if (b < 0 || a < 0 || b > a) return 0;
int res = 1;
for (int i = a;i >= a - b + 1;i--) res = i % mod * res % mod;
res = res * infact[b] % mod;
return res;
}
void Solve(int TIME) {
int m, k;cin >> m >> k;
int W = min(100ll, m);
vvc<int> f(W + 1, vi(k + 1));
f[0][0] = 1;
for (int w = 1;w <= W;w++) {
for (int i = 0;i <= k;i++) {
for (int j = 1;j <= 9;j++) {
if (i >= j) f[w][i] = (f[w][i] + f[w - 1][i - j]) % MOD;
}
}
}
int res = 0;
for (int i = 0;i <= W;i++) {
int cnt = f[i][k];
res = (res + cnt * comb(m, i) % MOD) % MOD;
}
if (k == 1) res = (res + 1) % MOD;
cout << res << endl;
}
14.5.数学尖子生
知识点
容斥
题解
容易推知，如果一个数$x$被$lcm(1,2,...n-1)$整除，但是不被$lcm(1,2,...n)$整除，那么这个数的$f(x)=n$。
于是预处理$Lcm$，然后容斥即可。
查看代码
void Solve(int TIME) {
int lc = 1;
vi Lcm(50);
for (int i = 1;i <= 42;i++) {
lc = lcm(lc, i);
Lcm[i] = lc;
}
int t;cin >> t;
while (t--) {
int a, n;cin >> a >> n;
if (a >= 42 || a == 1) {
cout << 0 << endl;
continue;
}
int res = n / Lcm[a - 1] - n / Lcm[a];
if (res <= 0) cout << 0 << endl;
else cout << res << endl;
}
}
15.4.取余
知识点
无
题解
如果范围小一点，可以直接数论分块$O(a\sqrt {b})$做，但是这里范围比较大所以考虑换其他方法。
考虑枚举$b$，这样每个合法的数的范围是$[bx+l,bx+r]$，其中$l=s$，$r=min(b-1,t)$，x=0，1，2...
然后根据这个范围计算一下有几个数符合条件，需要特判$l=0$和最后一个块的东西。不得不说调特判的能力是非常重要的！
查看代码
void Solve(int TIME) {
int A, B, s, t;cin >> A >> B >> s >> t;
int res = 0;
for (int b = 1;b <= B;b++) {
if (b <= s) continue;
int l = s, r = min(b - 1, t);
int x = A / b;
res += (r - l + 1) * x;
if (l == 0) res -= 1;
int t = A - b * x;
if (t < l) continue;
res += min(t, r) - l + 1;
}
cout << res << endl;
}
16.CF1861E - Non-Intersecting Subpermutations
知识点
容斥，划分型dp
题解
法1.
显然有个贪心的策略，应该尽量选择靠前的合法串。
考虑贡献的思想，以每个位置为结尾的合法段都求出贡献然后累加。
设$f_i$表示前$i$个位置中，最后一个长度为$k$的段以第$i$个位置为结尾的贡献次数。
显然可以先将$f_i=k!\cdot k^{i-k}$，然后减去不合法的方案。
于是不合法等价于，存在一个$j\in[i-k+1,i-1]$，使得$j$位置为结尾存在一个合法段。
于是方案数就是$f_j \cdot (i-j)!$，$(i-j)!$表示剩下的数字任意次序。
于是有$f_i=k!\cdot k^{i-k}-\sum\limits_{j=i-k+1}^{i-1}f_j \cdot(i-j)!$，后面的$n-i$个数随便选。
于是答案为$\sum\limits_{i=1}^{n}f_i\cdot k^{n-i}$
查看代码
void Solve(int TIME) {
int n, k;cin >> n >> k;
vi f(n + 1);
int res = 0;
for (int i = k;i <= n;i++) {
int t = 0;
for (int j = i - k + 1;j <= i - 1;j++) {
t = Add(t, mul(f[j], fact[i - j]));
}
f[i] = Dec(mul(fact[k], qp(k, i - k)), t);
res = Add(res, mul(f[i], qp(k, n - i)));
}
cout << res << endl;
}
法2.
从前往后扫，如果找到了一个长度为$k$的且满足条件的段就选择。
所以可以想到一个$dp$，设$f_{i,j}$表示前$i$个位置，末尾最长的，值只出现一次的子数组。
$f_{i+1,j+1} += f_{i,j}\cdot(k-j)$，$(i+1不是[i-j+1,i]出现的值)$
$f_{i+1,j} += f_{i,t}$ ，$t\in[j,k-1] $ ，$(i+1是[i-j+1，i]出现的值)$
$f_{i+1,0}=f_{i+1,k}$，插入一个值后所有的值都出现了，让长度归$0$重新开始。
于是答案就是$\sum\limits_{i=1}^{n}f_{i,k}\cdot k^{n-i}$
其实也用了贡献的思想。不过这里直接dp求出了上面法1容斥求出的东西。
查看代码
void Solve(int TIME) {
int n, k;cin >> n >> k;
vvc<int> f(n + 1, vi(k + 1));
f[0][0] = 1;
for (int i = 0;i < n;i++) {
for (int j = 0;j < k;j++) {
f[i + 1][j + 1] = Add(f[i + 1][j + 1], mul(f[i][j], (k - j)));
}
int s = 0;
for (int j = k - 1;j >= 1;j--) {
s = Add(s, f[i][j]);
f[i + 1][j] = Add(f[i + 1][j], s);
}
f[i + 1][0] = f[i + 1][k];
}
int res = 0;
for (int i = 1;i <= n;i++) {
res = Add(res, mul(f[i][k], qp(k, n - i)));
}
cout << res << endl;
}
17.CF1790E - Vlad and a Pair of Numbers
考虑把$a+b$换成位运算符号，才方便构造。
$a\text{^} =a|b-a\&b=x$
$a+b=a|b+a\&b=2x$
得 $a\&b=\frac{x}{2}$ ，又根据$a\text{^} b=x$
可以很容易构造解，分别看$x$和$y=\frac{x}{2}$的每一位数，如果$\{x,y\}=\{0,1\}$则$\{a,b\}=\{1,1\}$；如果$\{x,y\}=\{1,0\}$则$\{a,b\}=\{1,0\}或\{0,1\}$；如果$\{x,y\}=\{0,0\}$则$\{a,b\}=\{0,0\}$；如果$\{x,y\}=\{1,1\}$则无解。根据这个写一下即可。
查看代码
void Solve(int TIME) {
int x;cin >> x;
if (x & 1) return cout << -1 << endl, void();
int y = x / 2;
int a = 0, b = 0;
for (int i = 30;i >= 0;i--) {
int u = x >> i & 1;
int v = y >> i & 1;
if (u == 1 && v == 1) return cout << -1 << endl, void();
else if (u == 1 && v == 0) a += 1 << i;
else if (u == 0 && v == 1) a += 1 << i, b += 1 << i;
}
cout << a << " " << b << endl;
}
期望与概率
1. UVA12004 Bubble Sort
知识点
贡献计算
题解
考虑每个有序对对答案的贡献$f(i,j)$,其中$i<j$
$res=\sum\limits_{i=1}^{n}\sum\limits_{j=i+1}^{n}f(i,j)$。而$f(i,j)=1时$，$i$在$j$的后面。 $f(i,j)=0$时，i在j的前面。这两种情况的概率都是$\frac{1}{2}$
即$f(i,j)=\frac{1}{2}$， 那么$res=\frac{n(n+1)}{2}\cdot \frac{1}{2}=\frac{n(n+1)}{4}$
查看代码
void Solve(int TIME) {
int n;
cin >> n;
ll res = 1ll * n * (n - 1) / 2;
cout << "Case " << TIME << ": ";
if (res & 1) {
cout << res << "/" << 2 << endl;
}
else {
cout << res / 2 << endl;
}
}
2.ABC326E - Revenge of "The Salary of AtCoder Inc."
知识点 期望dp 题解 期望dp考虑逆推。设状态$f_i$表示当前的$x$为$i$时，终止时获得钱的期望数量
显然有$f_n=0$，因为此时无论抽中什么都会中终止。
$f_{n-1}=\frac{a_n}{n}$，即抽中$n$有贡献，否则终止。
$f_{n-2}=\frac{a_{n-1}+f_{n-1}}{n}+\frac{a_{n}}{n}$
$f_{n-3}=\frac{a_{n-2}+f_{n-2}}{n}+\frac{a_{n-1}+f_{n-1}}{n}+\frac{a_n}{n}$
然后发现$f_n=f_{n+1}+\frac{a_{n+1}+f_{n+1}}{n}$，倒着递推即可。
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
int invn = inv(n);
vi f(n + 1);//目前x=i,至结束获得钱的期望
//f[n-1]=a[n]/n+0,f[n-2]=(a[n-1]+f[n-1])/n+a[n]/n,f[n-3]=(a[n-2]+f[n-2])/n+(a[n-1]+f[n-1])n+(a[n])/n
f[n] = 0;
for (int i = n - 1;i >= 0;i--) {
f[i] = (f[i + 1] + (a[i + 1] + f[i + 1]) % MOD * invn % MOD) % MOD;
}
cout << f[0] << endl;
}
3.nc68504M - 让二追三
知识点
贡献计算期望
题解
考虑计算每个长度为5的连续子段的对答案贡献的期望，然后答案就是所有的这些单个的累加，即$(n-4)$个的累加，即$(n-4)\cdot (1-p)(1-p)(1-p)\cdot p\cdot p $
查看代码
void Solve(int TIME) {
int a, b, n;cin >> a >> b >> n;
int p = a * inv(b) % MOD;
if (n <= 4) return cout << 0 << endl, void();
int res = (n - 4) * (1 - p + MOD) % MOD * (1 - p + MOD) % MOD * p % MOD * p % MOD * p % MOD;
cout << res << endl;
}
4.牛客挑战赛71 - A 和的期望
知识点
贡献计算期望
题解
仅考虑单个数字对答案的贡献，如果要选出$i$个数字，那么就需要在选出这个数字后，再在剩下的$n-1$个数字选出$i-1$个。而总方案数是$C_{n}^{i}$，所以答案是$\frac{C_{n-1}^{i-1}S}{C_{n}^{i}}=\frac{(n-1)!}{(i-1)!(n-i)!} \frac{i!(n-i)!}{n!}S=\frac{i}{n}S$，其中$S=\sum\limits_{i=1}^{n}a_i$
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
int s = accumulate(a.begin() + 1, a.end(), 0ll);
s %= MOD;
for (int i = 1;i <= n;i++) {
cout << s * i % MOD * inv(n) % MOD << " ";
}
cout << endl;
}
5.abc332F - Random Update Query
知识点
计算期望
题解
发现每次对于一个区间里的任意一个数，设它的之前的期望是$E$。它的新的期望是$\frac{(r-l)\cdot E}{r-l+1}+\frac{x}{r-l+1}$
维护一个区间加乘线段树，每次给区间整体乘上$\frac{r-l}{r-l+1}$，然后再加上$\frac{x}{r-l+1}$
查看代码
struct Segtree {
#define ls (x << 1)
#define rs (x << 1 | 1)
struct Tag {
ll k = 1;
ll b = 0;
};
struct Info {
int sz;
ll sum = 0;
};
struct node {
Info info;
Tag tag;
};
Info friend operator +(const Info& l, const Info& r) {
return { l.sz + r.sz,l.sum + r.sum };
}
Info friend operator +(const Info& info, const Tag& tag) {
return { info.sz,(tag.k * info.sum + tag.b * info.sz) % MOD };
}
Tag friend operator+(const Tag& tag1, const Tag& tag2) {
return { tag1.k * tag2.k % MOD,(tag1.b * tag2.k + tag2.b) % MOD };
}
int n;
vector<node> tr;
Segtree(const vector<int>& a, int n) :n(n) {
tr.resize(n << 2);
build(a, 1, 1, n);
}
void build(const vector<int>& a, int x, int l, int r) {
tr[x].info.sz = r - l + 1;
if (l == r) {
tr[x].info = { 1,a[l] };
tr[x].tag = { 1,0 };
return;
}
else {
int mid = l + r >> 1;
build(a, ls, l, mid);
build(a, rs, mid + 1, r);
pushup(x);
}
}
void pushup(int x) {//从下往上更新
tr[x].info = tr[ls].info + tr[rs].info;
}
void settag(int x, Tag tag) {
tr[x].info = tr[x].info + tag;
tr[x].tag = tr[x].tag + tag;
}
void pushdown(int x) {//下传标记
settag(ls, tr[x].tag);
settag(rs, tr[x].tag);
tr[x].tag.k = 1;
tr[x].tag.b = 0;
}
//l,r代表操作区间
void update(int x, int l, int r, int ql, int qr, Tag tag) {
if (l == ql && r == qr) {
settag(x, tag);
return;
}
int mid = l + r >> 1;
pushdown(x);
if (qr <= mid) update(ls, l, mid, ql, qr, tag);
else if (mid < ql) update(rs, mid + 1, r, ql, qr, tag);
else {
update(ls, l, mid, ql, mid, tag);
update(rs, mid + 1, r, mid + 1, qr, tag);
}
pushup(x);
}
Info query(int x, int l, int r, int ql, int qr) {
if (l == ql && r == qr) return tr[x].info;
int mid = l + r >> 1;
pushdown(x);
if (qr <= mid) return query(ls, l, mid, ql, qr);
else if (mid < ql) return query(rs, mid + 1, r, ql, qr);
else return query(ls, l, mid, ql, mid) + query(rs, mid + 1, r, mid + 1, qr);
}
};
void Solve(int TIME) {
int n, m;cin >> n >> m;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
Segtree tr(a, n);
for (int i = 1;i <= m;i++) {
int l, r, x;cin >> l >> r >> x;
int len = r - l + 1;
tr.update(1, 1, n, l, r, { (MOD + 1 - inv(len)) % MOD, inv(len) * x % MOD });
}
for (int i = 1;i <= n;i++) {
cout << tr.query(1, 1, n, i, i).sum << " \n"[i == n];
}
}6.CF103637A - Agile permutation
知识点
置换环，斯特林排列数
题解
肯定是先一直洗牌，直到某个大于等于决策点不再洗牌，然后一直使用交换操作。
枚举这个决策点。决策点即需要交换的次数，在这里即为置换环的个数。
考虑置换环的个数，显然能用斯特林排列数表示。用$SA_{n,k}$表示$1\sim n$的排列中，环的个数为$k$的方案数。
每次同时计算出大于等于决策点的总值和次数，两者相除就是交换操作的期望代价。同时由于是几何分布，期望等于概率的倒数，计算出随机操作的期望代价。两者相加就是该决策点是期望。对所有决策点取最小值即为答案。
还有注意不随机也要计算一下答案，dfs得到置换环数量之后计入即可。
查看代码
void Solve(int TIME) {
int n, a, b;cin >> n >> a >> b;
vi p(n + 1);for (int i = 1;i <= n;i++) cin >> p[i];
vc<ld> fact(n + 1);fact[0] = 1.0;
for (int i = 1;i <= n;i++) {
fact[i] = fact[i - 1] * i;
}
vvc<int> SA(n + 1, vi(n + 1));
SA[0][0] = 1;
for (int i = 1;i <= n;i++)
for (int j = 1;j <= i;j++)
SA[i][j] = SA[i - 1][j - 1] + (i - 1) * SA[i - 1][j];
vi vis(n + 1);
auto dfs = [&](auto&& dfs, int u) {
if (vis[u]) return;
vis[u] = 1;
dfs(dfs, p[u]);
};
int tot = 0;
for (int i = 1;i <= n;i++) {
if (!vis[i]) dfs(dfs, i), tot++;
}
ld res = (n - tot) * a;
ld cnt = 0, val = 0;
for (int i = n;i >= 1;i--) {
cnt += SA[n][i], val += 1.0 * SA[n][i] * a * (n - i);
res = min(res, fact[n] / cnt * b + val / cnt);
}
cout << res << endl;
}
7.机器人
知识点
无穷级数$\frac{1}{1-x}$
题解
发现要停在一个地方，首先需要先到这个地方，然后转任意圈（可能为0）之后停在这个地方。预处理出转一圈的概率$P$。
设到达某地的概率为$q$，那么停在这个地方的概率是$q+qP+qP^2+...=\frac{q}{1-P}$。
排名的话，手玩一下发现排名依次降低（除了0），额外把为0的项放到最后。
查看代码
void Solve(int TIME) {
int n;cin >> n;
vi k(n + 1);for (int i = 1;i <= n;i++) cin >> k[i];
int inv2 = inv(2);
vi INV(5e6 + 1);INV[0] = 1;
for (int i = 1;i <= 5e6;i++) INV[i] = INV[i - 1] * inv2 % MOD;
vi p(n + 1);for (int i = 1;i <= n;i++) p[i] = INV[k[i]];
int pi = 1;
for (int i = 1;i <= n;i++) {
pi = pi * p[i] % MOD;
}
vi f(n + 1);
int now = 0;
int nnow = 0;
for (int i = 1;i <= n;i++) {
if (k[i]) f[i] = ++now;
else f[i] = ++nnow;
}
for (int i = 1;i <= n;i++) if (k[i] == 0) f[i] = f[i] + now;
int t = 1;
int res = 0, res2 = 0;
for (int i = 1;i <= n;i++) {
res = (res + i * t % MOD * (MOD + 1 - p[i]) % MOD * inv(MOD + 1 - pi) % MOD) % MOD;
res2 = (res2 + i * f[i] % MOD) % MOD;
t = t * p[i] % MOD;
}
cout << res << endl;
cout << res2 << endl;
}8.P9777 Fujisaki 讨厌数学
知识点
矩阵快速幂
题解
令$a_n=x^n+\frac{1}{x^n}$
$(x+\frac{1}{x})a_{n}=x^{n+1}+x^{n-1}+\frac{1}{x^{n-1}}+\frac{1}{x^{n+1}}=a_{n+1}+a_{n-1}$
即$a_n=(x+\frac1x)a_{n-1}-a_{n-2}=ka_{n-1}-a_{n-2}$
$\left[\begin{array}\ a_1 &a_2 \end{array}\right] \times\left[\begin{array}\ 0&-1\\1&k \end{array}\right]=\left[\begin{array}\ a_2&a_3 \end{array}\right]$
矩阵快速幂加速即可，$a_1=k,a_2=k^2-2$
查看代码
void Solve(int TIME) {
int k, n;cin >> MOD >> k >> n;
if (n == 0) return cout << 2 << endl, void();
mat<int> b;b(1, 1) = 0;b(1, 2) = MOD - 1, b(2, 1) = 1, b(2, 2) = k;
b = qp(b, n - 1);
mat<int> a;a(1, 1) = k, a(1, 2) = Dec(mul(k, k), 2);
a = a * b;
cout << a(1, 1) << endl;
}
9.D. Lonely Mountain Dungeons (codeforces.com)
知识点
凸函数
题解
固定了队伍数量，那么不同小队里同一物种的分配越均匀越好。
假设共$m$只队伍，$n$个人，设每个队的人数$x_1,x_2,...,x_m$，即$\frac12(x_1(n-x_1)+x_2\cdot(n-x_2)+...+x_m(n-x_m))=\frac12(n(x_1+...+x_m)-(x_1^2+...+x_m^2))$最大，即$\frac12(n^2-(x_1^2+...+x_m^2))$，显然平均分配是最好的。
数据范围很弱，随便模拟一下即可。
查看代码
void Solve(int TIME) {
int n, b, x;cin >> n >> b >> x;
map<int, int> mp;
vi c(n + 1);for (int i = 1;i <= n;i++) cin >> c[i], mp[c[i]]++;
sort(c.begin() + 1, c.end());
int mx = c[n];
auto calc = [&](int x, int gr) {
int avg = x / gr;
int re = x % gr;
int c1 = gr - re, c2 = re;
return (x * x - c1 * avg * avg - c2 * (avg + 1) * (avg + 1)) / 2;
};
int res = 0, now = 0;
for (int i = 1;i <= mx;i++) {
if (mp.count(i)) {
now += mp[i] * calc(i, i);
}
int tnow = now;
for (auto it = mp.upper_bound(i);it != mp.end();it++) {
auto [j, k] = *it;
tnow += k * calc(j, i);
}
res = max(res, tnow * b - (i - 1) * x);
}
cout << res << endl;
}打个表看出这是个凸函数，三分也可以做，不过不太会证明为什么是凸函数？
查看代码
void Solve(int TIME) {
int n, b, x;cin >> n >> b >> x;
map<int, int> mp;
vi c(n + 1);for (int i = 1;i <= n;i++) cin >> c[i], mp[c[i]]++;
sort(c.begin() + 1, c.end());
int mx = c[n];
auto calc = [&](int x, int gr) {
int avg = x / gr;
int re = x % gr;
int c1 = gr - re, c2 = re;
return (x * x - c1 * avg * avg - c2 * (avg + 1) * (avg + 1)) / 2;
};
auto check = [&](int k) {
if (k <= 0) return -inf;
int res = 0;
for (auto [i, j] : mp) {
res += j * calc(i, k);
}
return res * b - (k - 1) * x;
};
int l = 1, r = mx;
while (l <= r) {
int mid = l + r >> 1;
if (check(mid - 1) < check(mid)) l = mid + 1;
else r = mid - 1;
}
cout << check(l - 1) << endl;
}
10.6.集合划分【算法赛】 - 蓝桥云课 (lanqiao.cn)
知识点
容斥
题解
每个相同位上有$1$的，不能分到同一个集合。但是这样不好算，考虑它的反面，所有该位有$1$的，都分到同一个集合。然后容斥算一下即可，即考虑哪些位需要满足都分到一个集合。这个用并查集算即可。
查看代码
void Solve(int TIME) {
int n;cin >> n;
DSU d(n + 1);
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
int u = 0;for (int i = 1;i <= n;i++) u |= a[i];
auto check = [&](int x) {
vvc<int> del(15);
int one = 0;
for (int i = 0;i < 15;i++) {
if (x >> i & 1) {
one += 1;
for (int j = 1;j <= n;j++) {
if (a[j] >> i & 1) {
if (del[i].size()) d.merge(del[i].back(), j);
del[i].push_back(j);
}
}
}
}
int cnt = 0;
for (int i = 1;i <= n;i++) {
cnt += d.find(i) == i;
}
for (auto i : del) for (auto j : i) d.p[j] = j;
if (one & 1) return MOD - qp(2, cnt);
else return qp(2, cnt);
};
for (int i = 0;i < 15;i++) {
if (u >> i & 1) {
int cnt = 0;
for (int j = 1;j <= n;j++) {
if (a[j] >> i & 1) cnt++;
}
if (cnt == 1) {
return cout << 0 << endl, void();
}
}
}
int res = 0;
for (int mask = u;;mask = (mask - 1) & u) {
res = (res + check(mask)) % MOD;
if (mask == 0) break;
}
cout << res << endl;
}
由于本题数据很弱导致复杂度$2^n$的解也能过。
错解代码
void Solve(int TIME) {
int n;cin >> n;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
map<ai2, int> mp;
mp[{0, 0}] = 1;
int mask = 0;
int u = (1 << 15);
for (int i = 1;i <= n;i++) {
mask |= a[i];
map<ai2, int> nmp;
for (auto [c, k] : mp) {
auto [x, y] = c;
(nmp[{x | a[i], y}] += k) %= MOD;
(nmp[{x, y | a[i]}] += k) %= MOD;
}
swap(nmp, mp);
}
cout << mp[{mask, mask}] << endl;
}11.E - 7x7x7 (atcoder.jp)
知识点
容斥，公共部分经典trick
首先对于两个在$x$轴上的线段$(x_1,x_2),(x_3,x_4)$来说，公共部分的范围应该是$max(x_1,x_2)\le min(x_3,x_4)$。推广到二维的矩形和三维的立方体也是一样的道理。
由容斥原理，
$v3=V(a\cap b\cap c)$，$v_2=V(a\cap b)+V(b\cap c)+V(a\cap c)-3v_3,v_1=V_{总}-2v_2-3v_3$
暴力枚举第二个和第三个的位置，复杂度$14^6$，看看是否有位置符合。
查看代码
void Solve(int TIME) {
int V1, V2, V3;cin >> V1 >> V2 >> V3;
int V = 3 * 7 * 7 * 7;
auto func = [&](int x1, int y1, int z1, int x2, int y2, int z2) {
int res = 1;
res *= max(0ll, min(x1, x2) + 7 - max(x1, x2));
res *= max(0ll, min(y1, y2) + 7 - max(y1, y2));
res *= max(0ll, min(z1, z2) + 7 - max(z1, z2));
return res;
};
auto func2 = [&](int x1, int y1, int z1, int x2, int y2, int z2, int x3, int y3, int z3) {
int res = 1;
res *= max(0ll, min({ x1, x2,x3 }) + 7 - max({ x1, x2,x3 }));
res *= max(0ll, min({ y1, y2,y3 }) + 7 - max({ y1, y2,y3 }));
res *= max(0ll, min({ z1, z2,z3 }) + 7 - max({ z1, z2,z3 }));
return res;
};
for (int a1 = -7;a1 <= 7;a1++) {
for (int b1 = -7;b1 <= 7;b1++) {
for (int c1 = -7;c1 <= 7;c1++) {
for (int a2 = -7;a2 <= 7;a2++) {
for (int b2 = -7;b2 <= 7;b2++) {
for (int c2 = -7;c2 <= 7;c2++) {
int v3 = func2(0, 0, 0, a1, b1, c1, a2, b2, c2);
int v2 = func(a1, b1, c1, a2, b2, c2) + func(0, 0, 0, a1, b1, c1) + func(0, 0, 0, a2, b2, c2) - 3 * v3;
int v1 = V - 2 * v2 - 3 * v3;
if (v1 == V1 && v2 == V2 && v3 == V3) {
YES;
cout << 0 << " " << 0 << " " << 0 << " ";
cout << a1 << " " << b1 << " " << c1 << " ";
cout << a2 << " " << b2 << " " << c2 << endl;
return;
}
}
}
}
}
}
}
NO;
}
D. Exam in MAC (codeforces.com)
知识点
容斥
题解
对于$a_i$，有$\lfloor\frac{a_i}{2}\rfloor+1$对$(x,y)$，满足$x+y=a_i$。
有$c-a_i+1$对$(x,y)$满足$y-x=a_i$。
之后算一下同时满足两个条件的对数，分奇偶算一下即可。
然后容斥一下。
查看代码
void Solve(int TIME) {
int n, c;cin >> n >> c;
vi a(n + 1);for (int i = 1;i <= n;i++) cin >> a[i];
int res = (c + 1) * (c + 2) / 2;
ai2 par = { 0,0 };
for (int i = 1;i <= n;i++) {
res -= a[i] / 2 + 1;
res -= c - a[i] + 1;
par[a[i] & 1]++;
}
res += par[0] * (par[0] + 1) / 2;
res += par[1] * (par[1] + 1) / 2;
cout << res << endl;
}17届浙江省赛F
知识点
线性规划Simplex
题解
原本$x\in[-c,c]$，那么首先给所有的$x$减去一个$c$，使得$x$可以取$[0,2c]$。
接下来直接套板子即可。
$w_1x_1+...+w_mx_m+b_1>0$，$w_1'x_1+...+w_m'x_m+b_2<0$
加个$eps$给$AX\le B$变成$>0$和$<0$即可。
最后的答案记得减去$c$。
查看代码
const ld eps = 1e-10;
constexpr int MAXN = 405, MAXM = 405;
struct Simplex {
inline static ld a[MAXN][MAXM], b[MAXN], c[MAXM];//矩阵a,约束条件b,
inline static ld d[MAXN][MAXM], x[MAXM];
inline static int ix[MAXN + MAXM];
Simplex() {}
ld run(int n, int m) {
m++;
for (int i = 0; i < m + 1; i++) d[n][i] = d[n + 1][i] = 0;
for (int i = 0; i < n + m; i++) ix[i] = i;
int r = n, s = m - 1;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < m - 1; ++j) d[i][j] = -a[i][j];
d[i][m - 1] = 1;
d[i][m] = b[i];
if (d[r][m] > d[i][m]) r = i;
}
for (int i = 0;i < m - 1;i++) d[n][i] = c[i];
d[n + 1][m - 1] = -1;
for (ld dd; ; ) {
if (r < n) {
swap(ix[s], ix[r + m]);
d[r][s] = 1.0 / d[r][s];
for (int j = 0;j <= m;j++) if (j != s) d[r][j] *= -d[r][s];
for (int i = 0;i <= n + 1;i++) {
if (i != r) {
for (int j = 0;j <= m;j++) {
if (j != s) {
d[i][j] += d[r][j] * d[i][s];
}
}
d[i][s] *= d[r][s];
}
}
}
r = s = -1;
for (int j = 0;j < m;j++) {
if (s < 0 || ix[s] > ix[j]) {
if (d[n + 1][j] > eps || (d[n + 1][j] > -eps && d[n][j] > eps)) s = j;
}
}
if (s < 0) break;
for (int i = 0; i < n; ++i) {
if (d[i][s] < -eps) {
if (r < 0 || (dd = d[r][m] / d[r][s] - d[i][m] / d[i][s]) < -eps || (dd < eps && ix[r + m] > ix[i + m])) r = i;
}
}
if (r < 0) return -1;//无解
}
if (d[n + 1][m] < -eps) return -1;
ld ans = 0;
for (int i = 0;i < m;i++) x[i] = 0;
for (int i = m; i < n + m; i++) {
if (ix[i] < m - 1) {
ans += d[i - m][m] * c[ix[i]];
x[ix[i]] = d[i - m][m];
}
}
return ans;
}
};
const int C = 9000;
const ld eps2 = 1e-7;
void Solve(int TIME) {
int n;cin >> n;
vvc<int> A(2, vi(n));
vi B(2);
for (int i = 0;i < 2;i++) {
int sum = 0;
for (int j = 0;j < n;j++) {
cin >> A[i][j];sum += A[i][j];
}
cin >> B[i];
B[i] -= C * sum;
}
Simplex sim;
for (int cas = 0;cas < 2;cas++) {
for (int i = 0;i < n;i++) sim.c[i] = 1;
for (int i = 0;i < n;i++) for (int j = 0;j < n;j++) sim.a[i][j] = 0;
for (int j = 0;j < n;j++) {
sim.a[j][j] = 1;sim.b[j] = 2 * C;
}
for (int i = 0;i < 2;i++) {
for (int j = 0;j < n;j++) sim.a[n + i][j] = 0;
for (int j = 0;j < n;j++) {
if (cas == i) sim.a[n + i][j] = -A[i][j];
else sim.a[n + i][j] = A[i][j];
}
if (cas == i) sim.b[n + i] = B[i] - eps2;
else sim.b[n + i] = -B[i] - eps2;
}
ld res = sim.run(n + 2, n);
if (res == -1) continue;
for (int j = 0;j < n;j++) cout << sim.x[j] - C << " ";cout << endl;
return;
}
cout << "No" << endl;
}6.最后的加k【算法赛】 - 蓝桥云课 (lanqiao.cn)
知识点
矩阵光速幂
题解
注意到每一位是独立的，可以分别计算。列出$dp$方程，然后用矩阵光速幂预处理同底数同模数的矩阵，之后左乘初始矩阵得到答案。
（初始矩阵可以看作向量，这样常数会更快一点）
查看代码
template<typename T>
struct mat {
int m[11][11];
mat() {
memset(m, 0, sizeof m);
}
mat(int epsilon) {
memset(m, 0, sizeof m);
for (int i = 1;i <= 10;i++) m[i][i] = 1;
}
T& operator()(int i, int j) { return m[i][j]; }
T operator()(int i, int j)const { return m[i][j]; }
int is_I() {
for (int i = 1;i <= 10;i++) {
for (int j = 1;j <= 10;j++) {
if (i == j && m[i][j] != 1) return 0;
if (i != j && m[i][j] != 0) return 0;
}
}
return 1;
}
};
//矩阵乘法
template<typename T>
mat<T> operator *(const mat<T>& x, const mat<T>& y) {
mat<T> t;
for (int k = 1;k <= 10;k++) {
for (int i = 1;i <= 10;i++) {
for (int j = 1;j <= 10;j++) {
(t.m[i][j] += x.m[i][k] * y.m[k][j] % MOD) %= MOD;
}
}
}
return t;
}
//向量乘矩阵
template<typename T>
mat<T> operator &(const mat<T>& x, const mat<T>& y) {
mat<T> t;
for (int k = 1;k <= 10;k++) {
for (int i = 1;i <= 1;i++) {
for (int j = 1;j <= 10;j++) {
(t.m[i][j] += x.m[i][k] * y.m[k][j] % MOD) %= MOD;
}
}
}
return t;
}
//矩阵快速幂
mat<int> qp(mat<int> a, int k) {
mat<int> res(1);
while (k) {
if (k & 1) res = res * a;
a = a * a;
k >>= 1;
}
return res;
}
//光速幂,处理同底数同模数的幂
namespace LP_Mat {
ll getphi(ll x) {
ll res = x;
for (int i = 2; i * i <= x; i++) {
if (x % i == 0) {
res -= res / i;
while (x % i == 0) x /= i;
}
}
if (x > 1) res -= res / x;
return res;
}
mat<int> base1[N], basesqrt[N];
int Block_len;
int Phi;
ll maxn = 1e9 + 7;//模数的最大值
void init(mat<int> x) {//初始化底数为x
Phi = getphi(MOD);
Block_len = sqrt(maxn) + 1;
base1[0] = mat<int>(1);for (int i = 1;i <= Block_len;i++) base1[i] = base1[i - 1] * x;
basesqrt[0] = mat<int>(1);for (int i = 1;i <= Block_len;i++) basesqrt[i] = basesqrt[i - 1] * base1[Block_len];
}
mat<int> qp(ull x) {
x %= Phi;
return basesqrt[x / Block_len] * base1[x % Block_len];
}
}
void Solve(int TIME) {
int T, x;cin >> T >> x;
mat<int> base;
for (int i = 0;i <= 9;i++) {
int nx = x + i;
auto s = to_string(nx);
for (auto j : s) base(i + 1, j - '0' + 1) += 1;
}
LP_Mat::init(base);
while (T--) {
int n, k;cin >> n >> k;
mat<int> m;
for (auto i : to_string(n)) m(1, i - '0' + 1) += 1;
m = m & LP_Mat::qp(k);
int res = 0;
for (int i = 1;i <= 10;i++) res = (res + m(1, i)) % MOD;
cout << res << endl;
}
}P3746 [六省联考 2017] 组合数问题
知识点
矩阵快速幂
题解
$\sum\limits_{i=0}^{n}\left(\begin{array}\ nk\\ik+r \end{array}\right)=\sum\limits_{i=0}^{in+r}\left(\begin{array}\ nk\\i \end{array}\right)[i\equiv r(mod \; k)]$
即从$nk$个物品中选$i$个物品，且$i\equiv r(mod\;k)$的方案数。
于是考虑dp，$f_{i,j}$表示前$i$个物品，选择$j$个（注意需要在模$k$意义下）的方案数。
那么有$f_{i,j}=f_{i-1,j}+f_{i-1,(j-1)\;mod\; k}$，矩阵快速幂即可。注意特判k=1。
查看代码
void Solve(int TIME) {
mat<int> ma;
int n, k, r;cin >> n >> MOD >> k >> r;
if (k == 1) return cout << qp(2, n) << endl, void();
n *= k;
for (int i = 0;i < k;i++) {
ma(i, i) = ma((i - 1 + k) % k, i) = 1;
//ma(i, i) = ma(i, (i - 1 + k) % k) = 1;原本写反成这个也能通过！？思考一下数据水了还是为啥？
}
mat<int> res;
res(0, 0) = 1;
res = res * qp(ma, n);
cout << res(0, r) << endl;
}6.组合数【算法赛】 - 蓝桥云课 (lanqiao.cn)
双倍经验
查看代码
void Solve(int TIME) {
mat<int> ma;
int a, b, n;cin >> a >> b >> n;
if (a == 1) return cout << qp(2, n) << endl, void();
for (int i = 0;i < a;i++) {
ma(i, i) = ma((i - 1 + a) % a, i) = 1;
}
mat<int> res;
res(0, 0) = 1;
res = res * qp(ma, n);
int ans = res(0, b);
cout << ans << endl;
}
标签：lfloor,frac,limits,数论,sum,rfloor,int,数学,杂项 From： https://www.cnblogs.com/mrsuns/p/18133592