## 狄利克雷卷积+莫比乌斯反演学习笔记 ### 狄利克雷卷积 前置芝士：**卷积** 所谓的狄利克雷卷积，是定义在数论函数$(\Z^+\rightarrow C的函数)$间的一种二元运算，可这样定义： $(f*g)(n):=\sum_{xy=n}f(x)g(y)$ 也可以写成： $(f*g)(n):=\sum_{d|n}f(d)g(\frac{n}{d})$ 为了方便讨论先定义一些常用的数论函数符号： - 单位函数$\varepsilon(n):=\begin{cases} 1, \ while\ n = 1\\0,\ otherwise\end{cases}$ - 幂函数$Id_k(n):=n^k$.当$k = 1$时为恒等函数$Id(n)$，当$k = 0$时为常数函数$1(n)$ - 除数函数$\sigma_k(n):=\sum_{d|n}d^k$，当$k=1$是为因数和函数$\sigma(n)$，当$k=0$时为因数个数函数$\sigma_{0}(n)$ - 欧拉函数$\varphi(n)$ 值得注意的是，这些函数都是所谓的积性函数，即满足$f(1) = 1$，且若a，b互质，则有$f(a)f(b)=f(ab)$. 其中前两者还是完全积性函数，它们不需要满足a，b互质的条件也符合等式 积性函数之间的狄利克雷卷积还有特殊的性质： 若$f,g$为积性函数，那么$f*g$也是积性函数 首先$(f*g)(1)=f(1)*g(1)=1$，设$a,b$互质则有： $$ (f*g)(a)\cdot(f*g)(b)=\sum_{d1|a}f(d1)g(\frac{a}{d1})\cdot\sum_{d2|b}f(d2)g(\frac{b}{d2})\\ =\sum_{d1|a,d2|b}f(d1)f(d2)g(\frac{a}{d1})g(\frac{b}{d2})\\ =\sum_{d1|a,d2|b}f(d1d2)g(\frac{ab}{d1d2}) $$ 因为$a,b$互质，所以就有了$(f*g)(a) \cdot (f*g)(b)=(f*g)(ab)$该等式成立 再考虑某些数论函数之间的关系： #### 除数函数与幂函数 首先根据定义我们有$(f*1)(n)=\sum_{d|n}f(d)*1(\frac{n}{d})$ 所以可以有$(Id_k*1)(n) = \sum_{d|n}Id_k(d)1(\frac{n}{d})=\sum_{d|n}d^k=\sigma_k(n)$ 简写成$(Id_k*1)(x)=\sigma_k(x)$ #### 欧拉函数与恒等函数 $(\varphi*1)(n)=\sum_{d|n}\varphi(d)$ 当$n=p^m$时，有： $\sum_{d|n}\varphi(d)=\varphi(1)+\sum_{i=1}^m \varphi(p^i)=1+\sum_{i=1}^m(p^i-p^{i-1})=p^m=n$ 故$(\varphi * 1)(p^m)=p^m$ 现在令$n$为任意正整数，那么$n = \prod p^m$ 又因为$(\varphi*1)(n)$为积性函数，所以$(\varphi * 1)(n)=\prod p^m=n$ 所以有$(\varphi*1)=Id$ 下面对于一些卷积的法则进行证明： 交换律： $$ (f*g)(n)=\sum_{xy=n}f(x)g(y)\\ =\sum_{xy=n}f(y)g(x)\\ =(g*f)(n) $$ 结合律： $$ (f*(g*h))(n)=\sum_{xy=n}f(x)*(g*h)(y)\\ =\sum_{xy=n}f(x)*\sum_{uv=y}g(u)*h(v)\\ =\sum_{xuv=n}f(x)g(u)h(v)\\ =\sum_{yv=n}h(v)*\sum_{xu=y}f(x)g(u)\\ =((f*g)*h)(n) $$ 分配律： $$ (f*(g+h))(n)=\sum_{xy=n}f(x)*(g+h)(y)\\ =\sum_{xy=n}f(x)g(y) + f(x)h(y)\\ =\sum_{xy=n}f(x)g(y) + \sum_{xy=n}f(x)h(y)\\ =(f*g)(n)+(f*h)(n)\\ $$ 单位元的用法： $$ (f*\varepsilon)(n)=\sum_{d|n}f(d)\varepsilon(\frac{n}{d})=f(n)\\ $$ 逆元： 若$f*g=\varepsilon$时我们把$g$称为$f$的狄利克雷逆元并写作$f^{-1}$ $$ ( f * f ^ {-1})=\varepsilon $$ ### 莫比乌斯反演 我们定义数论函数$1$的狄利克雷逆元为$\mu$ 得到$\mu(x)=\begin{cases} 1 \ n = 1\\(-1)^m\ n = p_1p_2...p_m\\0\ otherwise\end{cases}$ 所以有$(f*1)=g$可推出$f=(g * \mu)$ 下面是几道经典题目的推导 #### P3455 [POI2007]ZAP-Queries ##### 简要题意 求$\sum_{x=1}^{a}\sum_{y=1}^b [gcd(x,y)=d]$ ##### 推导 $$ \sum_{x=1}^{a}\sum_{y=1}^b [gcd(x,y)=d]\\ =\sum_{x=1}^{\lfloor\frac{a}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{b}{d}\rfloor}[gcd(x,y)=1]\\ =\sum_{x=1}^{\lfloor\frac{a}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{b}{d}\rfloor}\sum_{k|gcd(x,y)}\mu(k)\\ =\sum_{k=1}^{n}\mu(k)\sum_{x=1}^{\lfloor\frac{a}{dk}\rfloor}\sum_{y=1}^{\lfloor\frac{b}{dk}\rfloor}1\\ =\sum_{k=1}^{n}\mu(k)\lfloor\frac{a}{dk}\rfloor\lfloor\frac{b}{dk}\rfloor\\ $$ 其中$n=max(a,b)$ 推导到最后一步就可以使用整除分块加速 #### P2522 [HAOI2011]Problem b ##### 简要题意 求$\sum_{i=x}^{a}\sum_{j=y}^b [gcd(i,j)=d]$ ##### 推导 这道题需要想到容斥原理。 将这个式子写成 $$ \sum_{i=1}^{a}\sum_{j=1}^b [gcd(i,j)=d]-\sum_{i=1}^{a}\sum_{j=1}^{y-1} [gcd(i,j)=d]-\sum_{i=1}^{x-1}\sum_{j=1}^b [gcd(i,j)=d]+\sum_{i=1}^{x-1}\sum_{j=1}^{y-1} [gcd(i,j)=d] $$ 就跟上面那道一样了 #### P2257 YY的GCD ##### 简要题意 求$\sum_{x=1}^{n}\sum_{y=1}^m [gcd(x,y)=p]$其中$p$为任意质数 ##### 推导 首先这里筛出质数后称作$prime$数量为$tot$ $$ \sum_{x=1}^{n}\sum_{y=1}^m [gcd(x,y)=p]\\ =\sum_{i=1}^{tot}\sum_{x=1}^{n}\sum_{y=1}^m [gcd(x,y)=prime[i]]\\ =\sum_{i=1}^{tot}\sum_{x=1}^{\lfloor\frac{n}{prime[i]}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{prime[i]}\rfloor} [gcd(x,y)=1]\\ =\sum_{i=1}^{tot}\sum_{d=1}^{max(n,m)}\mu(d)\lfloor\frac{n}{prime[i]*d}\rfloor\lfloor\frac{m}{prime[i]*d}\rfloor\\ 这里T=prime[i]*d可以将式子转化为\\ =\sum_{T=1}^{max(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{k|T,k\in prime}\mu(\frac{T}{k}) $$ 发现后面一个式子$\sum_{k|T,k\in prime}\mu(\frac{T}{k})$是可以预处理出来的，这样就减少了复杂度 #### P1829 [国家集训队]Crash的数字表格 / JZPTAB ##### 简要题意 求$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$ ##### 推导 $$ \sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)\\ =\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{gcd(i,j)}\\ =\sum_{d=1}^{max(n,m)}\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{d}[gcd(i,j)=d]\\ =\sum_{d=1}^{max(n,m)}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ijd[gcd(i,j)=1]\\ =\sum_{d=1}^{max(n,m)}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij[gcd(i,j)=1]\\ 设sum(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}ij[gcd(i,j)=1]\\ 化简这个式子\\ \sum_{i=1}^{n}\sum_{j=1}^{m}ij[gcd(i,j)=1]\\ =\sum_{d=1}^{max(n,m)}\mu(d)d^2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij\\ 设f(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}ij\\ 化简这个式子有\\ \sum_{i=1}^{n}\sum_{j=1}^{m}ij\\ =\sum_{i=1}^{n}i\sum_{j=1}^{m}j\\ =\frac{(n+1)n}{2}\cdot \frac{(m+1)m}{2}\\ 所以sum(n,m)=\sum_{d=1}^{max(n,m)}\mu(d)\cdot d^2\cdot f(\frac{n}{d},\frac{m}{d})\\ 原式=\sum_{d=1}^{max(n,m)}d\cdot sum(\frac{n}{d}, \frac{m}{d}) $$ 这里的两个式子都可以用整除分块+前缀和的方式加速 ##### 代码 ```cpp #include <bits/stdc++.h> const int N = 1e7 + 5, mod = 20101009; inline int add(int a, int b) { return (a += b) >= mod ? a - mod : a; } inline int sub(int a, int b) { return (a -= b) < 0 ? a + mod : a; } inline int mul(int a, int b) { return 1LL * a * b % mod; } std::vector < int > prime; int v[N], mu[N], sum[N]; void input(int n) { mu[1] = 1; for(int i = 2; i <= n; ++i) { if(!v[i]) prime.push_back(i), v[i] = i, mu[i] = mod - 1; for(auto x : prime) { if(x > n / i) break; v[i * x] = x; if(!(i % x)) { mu[x * i] = 0; break; } mu[i * x] = sub(mod, mu[i]); } } for(int i = 1; i <= n; ++i) sum[i] = add(sum[i - 1], mul(mu[i], mul(i, i))); } inline int sumij(int i, int j) { return mul((1LL * (i + 1) * i / 2) % mod, (1LL * (j + 1) * j / 2) % mod); } inline int f(int n, int m) { int re = 0; for(int l = 1, r; l <= m; l = r + 1) { r = std::min(n / (n / l), m / (m / l)); int cnt1 = n / l, cnt2 = m / l; re = add(re, mul(sumij(cnt1, cnt2), sub(sum[r], sum[l - 1]))); } return re; } int main() { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); int n, m; std::cin >> n >> m; if(n < m) std::swap(n, m); input(n); int ans = 0; for(int d = 1; d <= n; ++d) { ans = add(ans, mul(d, f(n / d, m / d))); } std::cout << ans << '\n'; return 0; } ``` ### 杜教筛 对于任意一个数论函数$g$，必满足： $$ \sum_{i=1}^{n}(f*g)(i)=\sum_{i=1}^{n}\sum_{d|i}g(d)f(\frac{i}{d})=\sum_{i=1}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor) $$ 重要推论： $$ g(1)S(n)=\sum_{i=1}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)-\sum_{i=2}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)\\ =\sum_{i=1}^{n}(f*g)(i)-\sum_{i=2}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor) $$ 如果可以构造一个数论函数$g$使得： - 可以快速计算$\sum_{i=1}^{n}(f*g)(i)$ - 可以快速计算$g$的前缀和，这样就能用数论分块求解$\sum_{i=2}^{n}g(i)S(\lfloor\frac{n}{i}\rfloor)$ 则我们可以取$m=\frac{2}{3}n$的复杂度求得$g(1)S(n)$ #### P3768 简单的数学题 $$ \sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j)\\ =\sum_{i=1}^{n}\sum_{j=1}^{n}ij\sum_{k|i,k|j}\varphi(k)\\ =\sum_{k=1}^{n}\varphi(k)\sum_{k|i}\sum_{k|i}ij\\ =\sum_{k=1}^{n}\varphi(k)\times k^2 \times (\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i)^2\\ =\sum_{k=1}^{n}\varphi(k)\times k^2 \times \sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}i^3 $$ 对$\varphi(k)k^2$做杜教筛 把$k$看成$Id$函数有，$f=Id^2\varphi$，卷一个$g=Id^2$可以得到$(f*g)(n)=n^3$ 这里有一个结论就是$Id^k$卷一个$Id^k$可以看成是$(Id^k*Id^k)(n)=n^k(n-\varphi(n))$ ### Powerful Number筛 Powerful Number（以下叫做PN）的定义：设$n=\prod_{i=1}^{m}p_i^{e_i}$只有所有的$e_i\gt 1$才能是PN 可以证明PN的数量级是$\Theta(\sqrt n)$的，也就是可以暴力搜索所有小于$\sqrt n$的质数处理出PN Powerful Number筛其实是杜教筛的升级版本，计算的也是$F(n)=\sum_{i=1}^n f(i)$（$f(i)$是积性函数），但是依然不适用于所有场景。 PN筛的步骤： - 找到一个积性函数$g(n)$在$n$为质数的时候$g(p)=f(p)$，并且$G(n)=\sum_{i=1}^n g(i)$ - 找到一个快速计算$G$也就是$g$的前缀和的方法 - 令$f = g*h$，此时有$h$为积性函数且$h(1) = 1$ - 最后$F(n)=\sum_{d=1,d\ is\ PN}^{n}h(d)G(\lfloor\frac{n}{d}\rfloor)$ 下面是对于这条公式的证明（可跳过）： 首先对于$h$函数有一个特殊性质： $$ f(p)=g(p)h(1)+g(1)h(p)=g(p)\\ h(p)=0 $$ 因为$h$是一个积性函数，我们知道$h(not\ PN)=0$ 于是 $$ F(n)=\sum_{i=1}^n \sum_{d|i}h(d)*g(\lfloor\frac{i}{d}\rfloor)\\ =\sum_{d=1}^n h(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}g(i)\\ =\sum_{d=1}^n h(d)G(\lfloor\frac{n}{d}\rfloor)\\ =\sum_{d=1,d\ is\ PN}^{n}h(d)G(\lfloor\frac{n}{d}\rfloor) $$ 得证 #### P5325【模板】Min_25 筛 给定$f(p^k)=p^k(p^k-1)$，求$\sum_{i=1}^n f(i)$ $$ 定义g(n)=Id(n)\varphi(n)\\ 此时G(n)可以通过杜教筛快速得出\\ G(n)=\sum_{i=1}^n i^2-\sum_{d=2}^n d*G(\lfloor\frac{n}{d}\rfloor)\\ 那么我们知道h(p^k)可以枚举计算（法1）\\ 另外也有h(p^k)直接推出的公式\\ f(p^k)=\sum_{i=0}^k g(p^{k-i})h(p^i)\\ p^k(p^k-1)=\sum_{i=0}^k p^{k-i}\varphi(p^i)h(p^i)\\ p^k(p^k-1)=\sum_{i=0}^k p^{2k-2i-1}(p-1)h(p^i)\\ h(p^k)=p^k(p^k-1)-\sum_{i=0}^{k-1} p^{2k-2i-1}h(p^i)\\ h(p^k)-p^2h(p^{k-1})=p^k(p^k-1)-p^{k+1}(p^{k-1}-1)-p(p-1)h(p^{k-1})\\ h(p^k)-ph(p^{k-1})=p^{k+1}-p^k\\ \frac{h(p^k)}{p^k}-\frac{h(p^{k-1})}{p^{k-1}}=p-1\\ 有h(p)=0推出\\ h(p^k)=(k-1)(p-1)p^k（法二）\\ 时间复杂度相似但是法二比较好写\\ $$ ```cpp void solve() { i64 n; std::cin >> n; int N1 = 10000000, N2 = sqrt(n); std::vector < int > prime, phi(N1 + 1), v(N1 + 1); std::vector < MInt > g(N1 + 1), s(N1 + 1); phi[1] = 1; for(int i = 2; i <= N1; ++i) { if(!v[i]) prime.emplace_back(i), v[i] = i, phi[i] = i - 1; for(auto x : prime) { if(x > N1 / i || x > v[i]) break; v[i * x] = x; phi[i * x] = phi[i] * (i % x == 0 ? x : x - 1); } } for(int i = 1; i <= N1; ++i) g[i] = 1LL * i * phi[i], s[i] = s[i - 1] + g[i]; MInt inv6 = MInt(6).inverse(), inv2 = 500000004; std::unordered_map < i64 , MInt > mp; auto G = [&](auto self, i64 n) -> MInt { if(n <= N1) return s[n]; if(mp.count(n)) return mp[n]; MInt re = MInt(n) * MInt(n + 1) * MInt(2 * n + 1) * inv6; for(i64 l = 2, r; l <= n; l = r + 1) { r = std::min(n, n / (n / l)); re -= (MInt(r) * MInt(r + 1) * inv2 - MInt(l - 1) * MInt(l) * inv2) * self(self, n / l); } return mp[n] = re; } ; // std::cout << g[2] << ' ' << g[3] << ' ' << g[6] << ' ' << g[9] << '\n'; MInt ans = 0; auto dfs = [&](auto self, int id, i64 now, MInt re) -> void { i64 cnt = 1LL * prime[id] * prime[id]; if(now > n / cnt) { // std::cout << now << ' ' << prime[id] << ' ' << n << '\n'; ans += re * G(G, n / now); return ; } self(self, id + 1, now, re); MInt st = cnt * g[prime[id]]; i64 las = n / now; while(cnt <= las) { MInt to = (MInt(cnt) * MInt(cnt - 1) - st); self(self, id + 1, now * cnt, re * to); if(cnt > las / prime[id]) break; st *= prime[id] * prime[id]; st += to * g[prime[id]]; cnt *= prime[id]; } } ; dfs(dfs, 0, 1, 1); std::cout << ans << '\n'; } ``` ```cpp void solve() { i64 n; std::cin >> n; int N1 = 10000000, N2 = sqrt(n); std::vector < int > prime, phi(N1 + 1), v(N1 + 1); std::vector < MInt > g(N1 + 1), s(N1 + 1); phi[1] = 1; for(int i = 2; i <= N1; ++i) { if(!v[i]) prime.emplace_back(i), v[i] = i, phi[i] = i - 1; for(auto x : prime) { if(x > N1 / i || x > v[i]) break; v[i * x] = x; phi[i * x] = phi[i] * (i % x == 0 ? x : x - 1); } } for(int i = 1; i <= N1; ++i) g[i] = 1LL * i * phi[i], s[i] = s[i - 1] + g[i]; MInt inv6 = MInt(6).inverse(), inv2 = 500000004; std::unordered_map < i64 , MInt > mp; auto G = [&](auto self, i64 n) -> MInt { if(n <= N1) return s[n]; if(mp.count(n)) return mp[n]; MInt re = MInt(n) * MInt(n + 1) * MInt(2 * n + 1) * inv6; for(i64 l = 2, r; l <= n; l = r + 1) { r = std::min(n, n / (n / l)); re -= (MInt(r) * MInt(r + 1) * inv2 - MInt(l - 1) * MInt(l) * inv2) * self(self, n / l); } return mp[n] = re; } ; // std::cout << g[2] << ' ' << g[3] << ' ' << g[6] << ' ' << g[9] << '\n'; MInt ans = 0; auto dfs = [&](auto self, int id, i64 now, MInt re) -> void { i64 cnt = 1LL * prime[id] * prime[id]; if(now > n / cnt) { // std::cout << now << ' ' << prime[id] << ' ' << n << '\n'; ans += re * G(G, n / now); return ; } self(self, id + 1, now, re); MInt i = prime[id] - 1; i64 las = n / now; while(cnt <= las) { self(self, id + 1, now * cnt, re * i * MInt(cnt)); if(cnt > las / prime[id]) break; i += MInt(prime[id] - 1); cnt *= prime[id]; } } ; dfs(dfs, 0, 1, 1); std::cout << ans << '\n'; } ``` ### Min_25筛 #### 使用要求 在$\Theta(\frac{n^{\frac{3}{4}}}{log\ n})$的时间复杂度下解决积性函数的前缀和问题。 要求$f(p)$是一个关于$p$的项数较少的多项式或可以快速求值，$f(p^c)$可以快速求值。 #### 定义说明 $p$没有特殊说明时代表全体质数 $isprime(n)$表示$n$为质数时为$1$ $p_k$表示全体质数中第$k$小的质数，特别的$p_0=1$ $lpf(n):=[1<n]min\{p:p|n\}+[1=n]$表示$n$的最小质因数。特别的，$lpf(1)=1$ $F_{prime}(n):=\sum_{2\le p \le n}f(p)$ $F_k(n):=\sum_{i=2}^n [p_k\le lpf(i)]f(i)$ #### 使用方法 观察$F_k(n)$的定义，我们要求$\sum_{i=1}^n f(i)=F_1(n)+f(1)=F_1(n)+1$ 考虑如何求出$F_k(n)$ $$ F_k(n)=\sum_{i=2}^n [p_k\le lpf(i)]f(i)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ =\sum_{k\le i,p_i^2\le n}\sum_{1\le c,p_i^c\le n}f(p^c)([c>1]+F_{i+1}(\lfloor\frac{n}{p_i^c}\rfloor))+\sum_{k\le i, p_i\le n}f(p_i)\\ =\sum_{k\le i,p_i^2\le n}(\sum_{1\le c,p_i^c\le n}f(p^c)([c>1]+F_{i+1}(\lfloor\frac{n}{p_i^c}\rfloor))+F_{prime}(n)-F_{prime}(p_{k-1}))\\ =\sum_{k\le i,p_i^2\le n}(\sum_{1\le c,p_i^{c+1}\le n}(f(p^c)F_{i+1}(\lfloor\frac{n}{p_i^c}\rfloor)+f(p_i^{c+1}))+F_{prime}(n)-F_{prime}(p_{k-1}))\\ $$ ##### 最后一步推导的说明 对于满足$p_i^c\le n< p_i^{c+1}$的$c$，有$\lfloor\frac{n}{p_i^c}\rfloor <p_i<p_{i+1}$。 故$F_{i+1}(\lfloor\frac{n}{p_i^c}\rfloor)=0$ 有$F_k(n)=0(p_k>n)$。 假设现在已经求出了所有的$F_{prime}(n)$，那么有两种方式可以求出所有的$F_k(n)$： - 直接按照递推式计算 - 从大到小枚举$p$转移，仅当$p^2<n$时转移增加值不为零，故按照递推式后缀和优化即可。 ##### 计算$F_{prime}(n)$ 一般情况下$f(p)$是一个关于$p$的低次多项式，可以表示成$f(p)=\sum a_ip^{c_i}$ 那么对于每个$p^{c_i}$，其对$F_{prime}(n)$的贡献可以表示为$a_i\sum_{2\le p\le n}p^{c_i}$ 分开考虑每个$p^{c_i}$的贡献，问题就转变为了：给定$n,s,g(p)=p^s$，对所有的$m=\lfloor\frac{n}{i}\rfloor$，求$\sum_{p\le m}g(p)$。 设： $$ G_x(n)=\sum_{2\le i \le n}[isprime(i) \or p_x<lpf(i)]g(i)\\ $$ 即埃氏筛第$k$轮晒完后剩下的数的$g$值之和。 对于一个合数$x$，必然有$lpf(x)\le\sqrt{x}$，则$F_{prime}=G_{\sqrt{n}}$故只需筛到$G_{\sqrt{n}}$即可 考虑$G$的边界值，$G_0(n)=\sum_{i=2}^n g(i)$ 对于转移，考虑埃氏筛的过程，分开讨论每部分的贡献，有： - 对于$n<p_k^2$的部分，$G$值不变，即$G_x(n)=G_{x-1}(n)$ - 对于$p_k^2\le n$部分，被筛掉的数必有质因子$p_k$即$-g(p_k)G_{k-1}(\lfloor\frac{n}{p_k}\rfloor)$ - 对于第二个部分，会有一些$lpf(i)<p_k$的$i$被减去，加回来的贡献是$g(p_k)G_{k-1}(p_{k-1})$ 则可以得到$G_k(n)=G_{k-1}(n)-[p_k^2\le n]g(p_k)(G_{k-1}(\lfloor\frac{n}{p_k}\rfloor)-G_{k-1}(p_{k-1}))$ ##### 总结 对于$F_k(n)$的计算，利用第一种方式写一个函数，这个实现难度较低，已经适用于大部分题目。 对于$F_{prime}的计算$，直接按递推式实现即可。 对于$p_k^2\le n$可以用线性筛预处理出$s_k:=F_{prime}(p_k)$来代替$F_k$递推式中的$F_{prime}(p_{k-1})$ 相应的$G_{k-1}(p_{k-1})=\sum_{i=1}^{k-1}g(p_i)$也可以预处理。 用埃氏筛求$f$的前缀和时，应当明确以下几点： - 如何快速（一般是线性复杂度）筛出前$\sqrt{n}$个$f$值 - $f(p)$的多项式表示 - 如何快速求出$f(p^c)$ 明确上述几点之后按顺序实现以下几部分即可： 1. 筛出$[1, \sqrt{n}]$内的质数与前$\sqrt{n}$个$f$值； 2. 对$f(p)$多项式表示中的每一项筛出对应的$G$，合并得到$F_{prime}$的所有$O(\sqrt{n})$有用点值 3. 按照$F_{k}$的递推式实现递归，求出$F_{1}(n)$ #### P5325【模板】Min_25 筛 ```cpp const int N = 2e5; i64 n, lim; void solve() { std::cin >> n; lim = sqrt(n); std::vector < int > prime, v(N + 10); std::vector < MInt > sum1(N + 10), sum2(N + 10); for(int i = 2; i <= N; ++i) { if(!v[i]) { prime.emplace_back(i); v[i] = i; } for(auto x : prime) { if(x > N / i || x > v[i]) break; v[i * x] = x; } } for(int i = 0; i < prime.size(); ++i) { sum1[i + 1] = sum1[i] + prime[i]; sum2[i + 1] = sum2[i] + MInt(1LL * prime[i] * prime[i]); } MInt inv2 = MInt(2).inverse(), inv6 = MInt(6).inverse(); auto f1 = [&](i64 n) -> MInt { return MInt(n) * MInt(n + 1) * inv2; } ; auto f2 = [&](i64 n) -> MInt { return MInt(n) * MInt(n + 1) * MInt(2 * n + 1) * inv6; } ; std::vector < i64 > idpre(N + 10), idsuf(N + 10), w; std::vector < std::vector < MInt > > G(2, std::vector < MInt > (N + 10)); for(i64 l = 1, r; l <= n; l = r + 1) { r = n / (n / l); w.emplace_back(n / l); G[0][w.size()] = f1(n / l) - 1, G[1][w.size()] = f2(n / l) - 1; if(n / l <= lim) idpre[n / l] = w.size(); else idsuf[r] = w.size(); } // std::cerr << '\n'; auto id = [&](i64 x) -> i64 { if(x <= lim) return idpre[x]; else return idsuf[n / x]; } ; for(int i = 0; i < prime.size(); ++i) { int x = prime[i]; for(int j = 1; j <= w.size() && 1LL * x * x <= w[j - 1]; ++j) { G[0][j] -= x * (G[0][id(w[j - 1] / x)] - sum1[i]); G[1][j] -= x * x * (G[1][id(w[j - 1] / x)] - sum2[i]); } } // std::cerr << '\n'; // std::cerr << "123" << '\n'; auto f = [&](MInt pc) -> MInt { return pc * (pc - 1); } ; auto F = [&](auto self, int k, i64 n) -> MInt { if(prime[k] > n) return 0; MInt ans = ((G[1][id(n)] - G[0][id(n)]) - (sum2[k] - sum1[k])); for(int i = k; i < prime.size() && prime[i] * prime[i] <= n; ++i) { for(i64 c = 1, sp = prime[i]; sp <= n / prime[i]; ++c, sp *= prime[i]) { ans += f(MInt(sp)) * self(self, i + 1, n / sp) + f(MInt(sp * prime[i])); } } return ans; } ; std::cout << F(F, 0, n) + 1 << '\n'; } ``` #### 2024GDCPC [另一个计数问题](https://cpc.csgrandeur.cn/csgoj/problemset/problem?pid=1286) ##### 题解 注意到，对于所有 $x \le \lfloor \frac{n + 1}{2} \rfloor$ 的结点 $x$，$x$ 与 $2x$ 间有一条边，而 $2x$ 与 $2$ 之间有一条边，因此所有 $2 \sim \lfloor \frac{n + 1}{2} \rfloor$ 范围内的结点连通。而对于 $x > \lfloor \frac{n + 1}{2} \rfloor$ 的结点 $x$，$x$ 不与 $2 \sim \lfloor \frac{n + 1}{2} \rfloor$ 中的结点连通当且仅当 $x$ 为素数。因此整个图仅包含若干个单独的素数结点与一个大的连通块。 考虑计算不连通的点对之间的贡献。不难发现，只需要求出所有 $> \lfloor \frac{n + 1}{2} \rfloor$ 的素数的和与平方和。使用 Min_25 筛或分块打表均可。 #### 2024杭电多校第2场 1004 a*b problem ##### 题意 求出$b=x_1^{t_1}x_2^{t_2}x_3^{t_3}...x_k^{t_k},1\le b\le n$在$gcd(x_1,x_2,...,x_k)=1$时， $(x_1,x_2,...,x_k)$有序对的数量。 ##### 题解 设$f(n)$表示$n=x_1^{t_1}x_2^{t_2}x_3^{t_3}...x_k^{t_k}$时$(x_1,x_2,...,x_k)$有序对的数量。 即答案为$\sum_{i=1}^n f(i)$。 如果在这里可以联想到Min_25筛的话，我们可以发现$f(p^c)$其实是与$p$无关的，即$p$的$0$次多项式。 此时我们知道$f(p)$就是$t_1,t_2,...,t_k$中$1$的数量，那么$F_{prime}(n)$就是$1$的数量与比$n$小质数个数的乘积 而$f(p^c)$我们知道每次放进$x_i^{t_i}$中，就是放一个$t_i$的倍数在其中。其实就是关于$c$的一个完全背包，但是这个完全背包有$gcd(x_1,x_2,...,x_k)=1$的限制，这个限制其实就是需要保证每个质数放的时候一定会至少有一个$0$的质量，魔改一下完全背包，加一个是否有放$0$就行。 剩下就变成板子了。 ```cpp const int N = 2e5; i64 n, lim; int m; MInt f[100001][34][2]; void solve() { std::cin >> n >> m; std::vector < int > t(m), p(34); lim = sqrt(n); for(auto &x : t) std::cin >> x, p[x]++; if(n == 1) { std::cout << "1\n"; return ; } std::vector < int > prime, v(N + 10); for(int i = 2; i <= N; ++i) { if(!v[i]) { prime.emplace_back(i); v[i] = i; } for(auto x : prime) { if(x > N / i || x > v[i]) break; v[i * x] = x; } } MInt inv2 = MInt(2).inverse(), inv6 = MInt(6).inverse(); std::vector < i64 > idpre(N + 10), idsuf(N + 10), w; std::vector < MInt > G(N + 10); for(i64 l = 1, r; l <= n; l = r + 1) { r = n / (n / l); w.emplace_back(n / l); G[w.size()] = MInt(n / l) - 1; if(n / l <= lim) idpre[n / l] = w.size(); else idsuf[r] = w.size(); } auto id = [&](i64 x) -> i64 { if(x <= lim) return idpre[x]; else return idsuf[n / x]; } ; for(int i = 0; i < prime.size(); ++i) { int x = prime[i]; for(int j = 1; j <= w.size() && 1LL * x * x <= w[j - 1]; ++j) { G[j] -= G[id(w[j - 1] / x)] - i; } } std::vector < MInt > Fprime(w.size() + 1); for(int i = 1; i <= w.size(); ++i) { Fprime[i] = G[i] * p[1]; } for(int i = 1; i <= m; ++i) { for(int j = 0; j <= 33; ++j) { f[i][j][0] = f[i][j][1] = 0; } } f[0][0][0] = 1; for(int i = 1; i <= m; ++i) { for(int j = 0; j <= 33; j += t[i - 1]) { if(j == 0) { for(int k = 0; k <= 33; ++k) { f[i][k][1] += f[i - 1][k][0] + f[i - 1][k][1]; } } else { for(int k = j; k <= 33; ++k) { f[i][k][0] += f[i - 1][k - j][0]; f[i][k][1] += f[i - 1][k - j][1]; } } } } auto dp = [&](int c) -> MInt { if(c > 33) return 0; return f[m][c][1]; } ; auto F = [&](auto self, int k, i64 n) -> MInt { if(prime[k] > n) return 0; MInt ans = Fprime[id(n)] - k * p[1]; for(int i = k; i < prime.size() && prime[i] * prime[i] <= n; ++i) { for(i64 c = 1, sp = prime[i]; sp <= n / prime[i]; ++c, sp *= prime[i]) { ans += dp(c) * self(self, i + 1, n / sp) + dp(c + 1); } } return ans; } ; std::cout << F(F, 0, n) + 1 << '\n'; } ```