@总结 - 13@ 筛
@0 - 参考资料@
zsy的线性筛blog
唐老师's杜教筛(orz)
jiry's杜教筛(orz)
fjzzq's杜教筛(orz)
fjzzq'sMin_25筛(orz)
一篇关于Min_25筛拓展功能的blog(orz)
@1 - 线性筛@
基本的线性筛咕且不提，我们在这里只讨论使用线性筛筛出任意的积性函数值（一般是两个积性函数的卷积）。
假如我们要筛 f(n)，将 n 唯一因数分解得到 \(n = p_1^{a_1}*p_2^{a_2}*...*p_k^{a_k}\)，且 \(p_1<p_2<...<p_k\)。
则：\(f(n) = f(p_1^{a_1})*f(p_2^{a_2})*...*f(p_k^{a_k})\)。
在线性筛中，我们通过一个数的最小值因子去筛掉 n。要是我们在筛的时候同时得到 \(p_1^{a_1}\) 的值，就可以愉快地筛出 f(n)。
首先：我们处理 f(1) = 1。
当我们筛出 p 为质数时，我们根据定义得到 f(p), f(p^2), ... f(p^c)，其中 c 是满足 p^c ≤ MAX 的最大值。这时往往可以利用前一个 f(p^(k-1)) 得到后一个 f(p^k)（比如欧拉函数）。
然后我们同时处理出 low(n) 表示 n 对应的 \(p_1^{a_1}\)。根据线性筛的具体过程，可以判断 a1 是否等于 1。如果等于 1，low(n) 就等于 p1；否则 low(n) = low(n/p1)*p1。
这个时候时间复杂度依然控制在线性复杂度。
@2 - 杜教筛@
显然上面那个不是我们研究的重点。我们不妨来看个问题进行引入：
给定 n（n ≤ 10^11），求 \(\sum_{i=1}^{n}\phi(i)\)。
乍一看，还以为是线性筛模板题。结果定睛一看，发现 n ≤ 10^11。这。。。可把人整懵逼了。
这类问题称作积性函数前缀和问题，而这类问题往往是近年来数论毒瘤的考察点。
像这时候就可以使用杜教筛（orz dyh）。
我们知道，关于 \(\phi\) 有一个很经典的狄利克雷卷积式：\(\phi * I = id\)。
而我们又知道，\(id\) 的前缀和是很好求的。
所以我们尝试将 \(\phi\) 转成 \(id\) 来求解，加快我们的速度：
不妨记 \(S(n) = \sum_{i=1}^{n}\phi(i)\)，则上式可以进一步转化：
这样就建立了 S(n) 与其他 S 之间的关系。
假如使用记忆化搜索，则用到的 \(\lfloor\frac{n}{i}\rfloor\) 只有 \(2*\sqrt{n}\) 种，且其中有 \(\sqrt{n}\) 种 \(\le \sqrt{n}\)。
我们分析一下时间复杂度，设 T(n) 表示求解 S(n) 的时间复杂度。因为用的是记忆化搜索，所以有：
用积分去拟合一下，得到 \(T(n) = O(n^{\frac{3}{4}})\)。
如果预处理前 k 个正整数的 S，且 \(k \ge \sqrt{n}\)，则可以得到：
当 \(k = O(n^{\frac{2}{3}})\) 时，有 \(O(k) = O(\frac{n}{\sqrt{k}})\)，此时复杂度平衡至最优。
一般地，假如我们要求某数论函数 \(f(n)\) （注：这里不一定是积性函数。只是因为积性函数方便预处理。）的前缀和 \(S(n)\)，我们可以尝试构造函数 \(g(n)\) 与 \(h(n)\) 使得 \(f*g = h\)，且 \(h\) 的前缀和是很好求的。
则仿照上面的欧拉函数前缀和，可以推导出如下结果：
预处理出 S 的前 \(O(n^{\frac{2}{3}})\) 项即可实行杜教筛。
当然，你问我这个 \(g\) 怎么构造。。。emmm这我也不好说，自己感受一下就可以构造出来了。
（因为 \(g\) 的构造往往是难点所在啊喂）
一个例题。
一份参考代码：
#include<cstdio>
#include<cstring>
typedef long long ll;
const int MAXN = 5000000;
const int MOD = 1000000007;
ll sphi[MAXN + 5]; int phi[MAXN + 5];
int prm[MAXN + 5], pcnt = 0;
bool is_prm[MAXN + 5];
void sieve() {
phi[1] = sphi[1] = 1;
for(int i=2;i<=MAXN;i++) {
if( !is_prm[i] ) {
prm[++pcnt] = i;
phi[i] = i-1;
}
for(int j=1;i*prm[j]<=MAXN;j++) {
is_prm[i*prm[j]] = true;
if( i % prm[j] == 0 ) {
phi[i*prm[j]] = phi[i]*prm[j];
break;
}
else phi[i*prm[j]] = phi[i]*phi[prm[j]];
}
sphi[i] = (sphi[i-1] + phi[i])%MOD;
}
}
bool tag_phi[MAXN + 5]; ll dp_phi[MAXN + 5];
ll solve_phi(ll n, int k) {
if( n <= MAXN ) return sphi[n];
if( tag_phi[k] ) return dp_phi[k];
ll ret = (n%MOD)*((n+1)%MOD)%MOD*500000004%MOD; ll l = 2;
while( l <= n ) {
ll r = n/(n/l);
ret = ret - solve_phi(n/l, k*l)*(r-l+1)%MOD;
l = r + 1;
}
tag_phi[k] = true;
return dp_phi[k] = ret;
}
int main() {
sieve(); ll N;
scanf("%lld", &N);
printf("%lld\n", (solve_phi(N, 1)%MOD + MOD)%MOD);
}
@3 - Min-25 筛@
（P.S：目前博主只会基本的 Min-25 操作，更高深的拓展还不会。。。）
杜教筛太依赖于函数本身的性质进行构造了。
假如给定任意的积性函数 f(n)，我们应该怎么求 f 的前缀和呢？
曾经有一个叫洲阁筛的能够实现 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的复杂度。
然后被 Min-25 筛的玄学复杂度爆踩，成为时代的眼泪。
事实上，zzt dalao 证明了在 \(n \leq 10^{13}\) 内洲阁筛跑得过的数据，Min-25 筛都跑得过。
该算法的优化思路为：合数总存在一个 \(\le \sqrt{n}\) 的质因子。
不过这个算法要是高效，前提是：
（1）f(p) 是一个多项式（等会儿你就知道为什么要求这个了）。
（2）f(p^k) 能够容易算出（k > 1）。
首先我们需要处理出 g(n)，表示 n 以内的素数的积性函数和。为了得到 g(n)，我们可以逐步递推 h(i, n)：
注：这个以及下面的 f(k) 并不是它的真实积性函数值，而是根据 f(p) 的多项式表达式算出来。
实现时，必须 f(k) 拆成若干个单项式分别求和，这样才能同时满足积性函数与多项式两条性质。
其实就是一个埃氏筛法的过程。h 的转移分为以下两类：
（1）当 \(prime[i]^2 \le n\) 时：
因为 h(i, n) 包含两个部分：质数、最小值质因子 \(>prime[i]\) 的，所以要减去质数的情况。
我们通过这个操作就可以筛去 \(prime[i]\) 的倍数中还未筛掉的数。
（2）当 \(prime[i]^2 > n\) 时：
这个很显然。根据埃筛的过程，这个情况根本不会筛去任何数。
可以发现第（2）种情况可以直接滚动数组滚掉。
设 \(\le \sqrt{n}\) 的质数个数为 pcnt，则最后的 \(h(pcnt, n)\) 即为 \(g(n)\)
同时，一样的套路，\(\lfloor\frac{n}{prime[i]}\rfloor\) 只有 \(O(\sqrt{n})\) 种可能的取值。
这个过程的复杂度？是 \(O(\frac{n^{\frac{3}{4}}}{\log n})\) 的。证明？抱歉我太菜了我不会。。。
现在的问题是，我们怎么预处理出 \(h(0, n)\)。
因为 1 实在太特殊了，我们仅令 \(h(0, n) = \sum_{i=2}^{n}f(i)\)，最后的最后再特殊处理 f(1)。
现在可以发现：假如说 f 的表达式是多项式，我们就可以愉快地一波自然数幂和。
好的。现在我们已经处理出 \(g(n)\) 表示 n 以内的质数积性函数之和。
现在我们来处理 \(S(i, n)\)，类似地，定义为：
（注：这里的 f 就是它的真实积性函数值了）
首先处理质数部分，即：
这个比较容易理解。
然后合数部分，我们去枚举合数的最小质因子及其幂：
看起来有点长，其实因为合数也分为两种：一种是形如 p^k (k > 1) 式的，另一种则包含两种以上的质因子。
因为我们 f(1) 是特殊计算的，所以这种分类是必要的。
最后就可以得到 \(S(i, n) = A + B\)。计算合数的复杂度就是 O(玄学)。
一个例题。
可以发现除了 f(2) = 3 以外，f(p) = p - 1（都是奇数）。
所以只需要最后特殊处理一下 f(2)。
参考代码：
#include<cstdio>
#include<cmath>
typedef long long ll;
const int MOD = int(1E9) + 7;
const int SQRT = int(1E5);
const int INV2 = (MOD + 1)>>1;
inline int add(ll x, ll y) {return (x%MOD + y%MOD)%MOD;}
inline int sub(ll x, ll y) {return add(x%MOD, MOD - y%MOD);}
inline int mul(ll x, ll y) {return (x%MOD)*(y%MOD)%MOD;}
ll n, sq;
int prm[SQRT + 5], pcnt;
bool nprm[SQRT + 5];
void sieve() {
for(int i=2;i<=SQRT;i++) {
if( !nprm[i] ) prm[++pcnt] = i;
for(int j=1;i*prm[j]<=SQRT;j++) {
nprm[i*prm[j]] = true;
if( i % prm[j] == 0 ) break;
}
}
}
ll arr[2*SQRT + 5]; int num1[SQRT + 5], num2[SQRT + 5], cnt;
inline int id(ll x) {return x > sq ? num1[n/x] : num2[x];}
int g[2*SQRT + 5], h[2*SQRT + 5], G[2*SQRT + 5];
void getG() {
for(ll x=1;x<=n;x=n/(n/x)+1) {
arr[++cnt] = n/x;
(n/x > sq ? num1[n/(n/x)] : num2[n/x]) = cnt;
}
for(int i=1;i<=cnt;i++) {
g[i] = sub(mul(add(mul(arr[i], arr[i]), arr[i]), INV2), 1);
h[i] = sub(arr[i], 1);
}
for(int i=1;i<=pcnt;i++)
for(int j=1;j<=cnt&&prm[i]<=arr[j]/prm[i];j++) {
g[j] = sub(g[j], mul(prm[i], sub(g[id(arr[j]/prm[i])], g[id(prm[i-1])])));
h[j] = sub(h[j], sub(h[id(arr[j]/prm[i])], h[id(prm[i-1])]));
}
for(int i=1;i<=cnt;i++) {
G[i] = sub(g[i], h[i]);
if( arr[i] >= 2 ) G[i] = add(G[i], 2);
}
}
int min_25(int n, int k) {
if( arr[n] < prm[k] ) return 0;
int ret = sub(G[n], G[id(prm[k-1])]);
for(int i=k;i<=pcnt&&prm[i]<=arr[n]/prm[i];i++)
for(ll c=1,j=prm[i];j<=arr[n]/prm[i];c++,j*=prm[i])
ret = add(ret, add(mul(prm[i]^c, min_25(id(arr[n]/j), i+1)), prm[i]^(c+1)));
return ret;
}
int main() {
scanf("%lld", &n), sq = sqrt(n);
sieve(), getG();
printf("%d\n", min_25(1, 1) + 1);
}
@4 - powerful number@
WC2019 介绍的，基本不会，就咕了。update in 2020/03/11：会了一点，来更新了。
首先 powerful number 的定义是：唯一分解中每个质因子的幂次 >= 2 的数。
可以发现它可以表示成 \(x = a^2b^3\)（表示方法好像并不唯一），因此估计出 n 以内的 powerful number 数量为 \(O(\sum_{i=1}^{\sqrt{n}} (\frac{n}{i^2})^{\frac{1}{3}})\)，积分得到大概是 \(O(\sqrt{n})\)。
利用 powerful number 的数量很少这一性质可以帮助我们更快地求解积性函数前缀和。
采用拟合的思想。如果我们要求积性函数 F(x) 前缀和，且存在易求前缀和的积性函数 G(x) 满足 G(p) = F(p)（p 为质数）。此时作 H 满足 G*H = F。
由于 F, G 都是积性函数，可以推得 H 也是积性函数，那么 H(1) = 1。
同时由于 H(1)*G(p) + H(p)*G(1) = F(p)，可以得到 H(p) = 0。也就是说 H(x) ≠ 0 当且仅当 x 是 powerful number。
接着推式子：\(\sum_{i=1}^{n}F(i) = \sum_{i=1}^{n}H(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}G(j)\)。
剩下的只需要搜索出所有 powerful number 然后算出 H(i) 的贡献就好了。
举个例子：PE639：注意到 \(f_k(p) = p^k\)，构造 \(g_k(x) = x^k\) 即可。g 求前缀和就是一个自然数幂和的事。剩下的按上面所说来就可以 \(O(k^2\sqrt{n})\) 算出答案来。
再举个例子：PE484：分析一波发现 \(f(k) = \gcd(k, k')\) 是一个积性函数，且 \(f(p^c) = \gcd(p^c, c\times p^{c-1}) = p^{c - [c\mod p \not = 0]}\)。那么有 \(f(p) = 1\)，构造 \(g(x) = 1\) 即可。
嘛，这道题还有一点特殊性：\(h(x) = f(x)\times \mu(x)\)。那么可以根据 \(\mu(p^c) = 0 (c > 1)\) 更快地算出 \(h\)。
最后一个例子：loj6053：用 min_25 做的时候已经分析过，除了 p = 2 以外 \(f(p^c) = p - 1\)。构造 \(g(x) = \phi(x)\) 即可，\(g(x)\) 的前缀和可以min_25筛杜教筛。对于 p = 2 的特例，此时造成的影响只有 h(2) ≠ 0，特判一下 2 的时候即可。
//为什么我的程序跑得比 min_25 还慢（困惑.jpg）
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MOD = 1000000007;
const int INV2 = (MOD + 1) / 2;
const int MAXN = 4641600;
const int PMAX = 325161;
inline int add(int x, int y) {return (x + y >= MOD ? x + y - MOD : x + y);}
inline int sub(int x, int y) {return (x - y < 0 ? x - y + MOD : x - y);}
inline int mul(int x, int y) {return 1LL * x * y % MOD;}
bool nprm[MAXN + 5];
int phi[MAXN + 5], prm[MAXN + 5], pcnt;
void init() {
phi[1] = 1;
for(int i=2;i<=MAXN;i++) {
if( !nprm[i] ) prm[++pcnt] = i, phi[i] = i - 1;
for(int j=1;i*prm[j]<=MAXN;j++) {
nprm[i*prm[j]] = true;
if( i % prm[j] == 0 ) {
phi[i*prm[j]] = phi[i]*prm[j];
break;
}
else phi[i*prm[j]] = phi[i]*phi[prm[j]];
}
}
for(int i=1;i<=MAXN;i++) phi[i] = add(phi[i - 1], phi[i]);
}
ll n;
int f[MAXN + 5], id1[MAXN + 5], id2[MAXN + 5], cnt;
int id(ll x) {return (x <= MAXN ? id1[x] : id2[n/x]);}
void get() {
cnt = 0;
for(ll i=1;i<=n;i=(n/(n/i))+1) {
ll p = n / i;
if( p <= MAXN ) id1[p] = (++cnt);
else id2[n/p] = (++cnt);
f[cnt] = -1;
}
}
int sum(ll m) {
if( m <= MAXN ) return phi[m];
if( f[id(m)] != -1 ) return f[id(m)];
int &x = f[id(m)]; x = mul(mul(m % MOD, (m + 1) % MOD), INV2);
for(ll i=2;i<=m;) {
ll p = m / i, j = m / p;
x = sub(x, mul((j - i + 1) % MOD, sum(p)));
i = j + 1;
}
return x;
}
int g[PMAX + 5][35], h[PMAX + 5][35], ans;
void dfs(ll m, int d, int k) {
ans = add(ans, mul(k, sum(n/m)));
for(int i=d;i<=pcnt;i++) {
if( prm[i] != 2 && (n < 1LL*prm[i]*prm[i]*m) ) break;
ll p = m * prm[i];
for(int j=1;p<=n;j++,p*=prm[i])
if( h[i][j] ) dfs(p, i + 1, mul(k, h[i][j]));
}
}
int main() {
init(), scanf("%lld", &n), get();
for(int i=1;i<=pcnt;i++) {
ll p = 1LL * prm[i] * prm[i];
if( p > n ) break;
g[i][0] = h[i][0] = 1;
g[i][1] = prm[i] - 1, h[i][1] = sub(prm[i] ^ 1, g[i][1]);
for(int j=2;p<=n;p*=prm[i],j++) {
g[i][j] = mul(g[i][j-1], prm[i]), h[i][j] = (prm[i] ^ j);
for(int k=0;k<j;k++)
h[i][j] = sub(h[i][j], mul(h[i][k], g[i][j-k]));
}
}
dfs(1, 1, 1), printf("%d\n", ans);
}
不过这个方法虽然好，但是局限性不是一般地大。。。不过 zzt 好像在 WC2019 的时候讲过怎么拓展这个方法。
但是我什么也没学会.jpg。
