狄利克雷卷積
前置
下取整的性質
$\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."
設狀態$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;
}