@总结

目录

• @0 - 参考资料@
• @1 - 线性筛@
• @2 - 杜教筛@
• @3 - Min-25 筛@
• @4 - powerful number@

---

@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) 来求解，加快我们的速度：

[sum_{i=1}^{n}id(i) = sum_{i=1}^{n}sum_{d|i}phi(frac{i}{d})*I(d) ]

[frac{n*(n-1)}{2} = sum_{a*ble n}phi(a) = sum_{b=1}^{n}sum_{i=1}^{lfloorfrac{n}{b} floor}phi(i) ]

不妨记 (S(n) = sum_{i=1}^{n}phi(i))，则上式可以进一步转化：

[S(n) = frac{n*(n-1)}{2} - sum_{i=2}^{n}S(lfloorfrac{n}{i} floor) ]

这样就建立了 S(n) 与其他 S 之间的关系。
假如使用记忆化搜索，则用到的 (lfloorfrac{n}{i} floor) 只有 (2*sqrt{n}) 种，且其中有 (sqrt{n}) 种 (le sqrt{n})。

我们分析一下时间复杂度，设 T(n) 表示求解 S(n) 的时间复杂度。因为用的是记忆化搜索，所以有：

[T(n) = sum_{i=1}^{sqrt{n}}(sqrt{i} + sqrt{frac{n}{i}}) ]

用积分去拟合一下，得到 (T(n) = O(n^{frac{3}{4}}))。
如果预处理前 k 个正整数的 S，且 (k ge sqrt{n})，则可以得到：

[T(n) = sum_{i=1}^{frac{n}{k}}sqrt{frac{n}{i}} = O(frac{n}{sqrt{k}}) ]

当 (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) 的前缀和是很好求的。
则仿照上面的欧拉函数前缀和，可以推导出如下结果：

[g(1)*S(n) = sum_{i=1}^{n}h(i) - sum_{i=2}^{n}g(i)*S(lfloorfrac{n}{i} floor) ]

预处理出 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
", (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)：

[h(i, n) = sum_{k 为质数||k 的最小质因子 > prime[i]}^{1 < kle n}f(k) ]

注：这个以及下面的 f(k) 并不是它的真实积性函数值，而是根据 f(p) 的多项式表达式算出来。
实现时，必须 f(k) 拆成若干个单项式分别求和，这样才能同时满足积性函数与多项式两条性质。

其实就是一个埃氏筛法的过程。h 的转移分为以下两类：
（1）当 (prime[i]^2 le n) 时：

[h(i, n) = h(i-1,n) - f(prime[i])*(h(i-1,lfloorfrac{n}{prime[i]} floor) - sum_{k=1}^{i-1}f(prime[k])) ]

因为 h(i, n) 包含两个部分：质数、最小值质因子 (>prime[i]) 的，所以要减去质数的情况。
我们通过这个操作就可以筛去 (prime[i]) 的倍数中还未筛掉的数。

（2）当 (prime[i]^2 > n) 时：

[h(i, n) = h(i-1, n) ]

这个很显然。根据埃筛的过程，这个情况根本不会筛去任何数。
可以发现第（2）种情况可以直接滚动数组滚掉。

设 (le sqrt{n}) 的质数个数为 pcnt，则最后的 (h(pcnt, n)) 即为 (g(n))
同时，一样的套路，(lfloorfrac{n}{prime[i]} floor) 只有 (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))，类似地，定义为：

[S(i, n) = sum_{j 的最小质因子 ge prime[i]}^{1 < j le n}f(j) ]

（注：这里的 f 就是它的真实积性函数值了）

首先处理质数部分，即：

[A = g(n, pcnt) - sum_{j=1}^{i-1}f(prime[j]) ]

这个比较容易理解。

然后合数部分，我们去枚举合数的最小质因子及其幂：

[B = sum_{j=i}^{prime[j]^2 le n}sum_{e=1}^{prime[j]^{e+1} le n}(f(prime[j]^e)*S(lfloorfrac{n}{prime[j]^e} floor, j+1) + f(prime[j]^{e+1})) ]

看起来有点长，其实因为合数也分为两种：一种是形如 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
", min_25(1, 1) + 1);
}

@4 - powerful number@

WC2019 介绍的，基本不会，就咕了。update in 2020/03/11：会了一点，来更新了。

参考的zzq的博客(orz)。

首先 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}^{lfloorfrac{n}{i} floor}G(j))。
剩下的只需要搜索出所有 powerful number 然后算出 H(i) 的贡献就好了。

举个例子：PE639：注意到 (f_k(p) = p^k)，构造 (g_k(x) = x^k) 即可。g 求前缀和就是一个自然数幂和的事。剩下的按上面所说来就可以 (O(k^2sqrt{n})) 算出答案来。

再举个例子：PE484：分析一波发现 (f(k) = gcd(k, k')) 是一个积性函数，且 (f(p^c) = gcd(p^c, c imes p^{c-1}) = p^{c - [cmod p ot = 0]})。那么有 (f(p) = 1)，构造 (g(x) = 1) 即可。
嘛，这道题还有一点特殊性：(h(x) = f(x) imes 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
", ans);
}

不过这个方法虽然好，但是局限性不是一般地大。。。不过 zzt 好像在 WC2019 的时候讲过怎么拓展这个方法。
但是我什么也没学会.jpg。