@总结 - 13@ 筛

た入场券 2023-06-01 14:56 5阅读 0赞

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

--------------------

## @0 - 参考资料@ ##

[zsy的线性筛blog][zsy_blog]  
[唐老师's杜教筛(orz)][s_orz]  
[jiry's杜教筛(orz)][jiry_s_orz]  
[fjzzq's杜教筛(orz)][fjzzq_s_orz]  
[fjzzq'sMin\_25筛(orz)][fjzzq_sMin_25_orz]  
[一篇关于Min\_25筛拓展功能的blog(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。

这个时候时间复杂度依然控制在线性复杂度。

[一个需要任意筛线性函数的例子。][Link 1]

## @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\*b\\le n\}\\phi(a) = \\sum\_\{b=1\}^\{n\}\\sum\_\{i=1\}^\{\\lfloor\\frac\{n\}\{b\}\\rfloor\}\\phi(i)\\\]

不妨记 \\(S(n) = \\sum\_\{i=1\}^\{n\}\\phi(i)\\)，则上式可以进一步转化：  
\\\[S(n) = \\frac\{n\*(n-1)\}\{2\} - \\sum\_\{i=2\}^\{n\}S(\\lfloor\\frac\{n\}\{i\}\\rfloor)\\\]

这样就建立了 S(n) 与其他 S 之间的关系。  
假如使用记忆化搜索，则用到的 \\(\\lfloor\\frac\{n\}\{i\}\\rfloor\\) 只有 \\(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(\\lfloor\\frac\{n\}\{i\}\\rfloor)\\\]

预处理出 S 的前 \\(O(n^\{\\frac\{2\}\{3\}\})\\) 项即可实行杜教筛。  
当然，你问我这个 \\(g\\) 怎么构造。。。emmm这我也不好说，自己感受一下就可以构造出来了。  
（因为 \\(g\\) 的构造往往是难点所在啊喂）

[一个例题。][Link 2]  
一份参考代码：

#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)：  
\\\[h(i, n) = \\sum\_\{k 为质数||k 的最小质因子 > prime\[i\]\}^\{1 < k\\le 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,\\lfloor\\frac\{n\}\{prime\[i\]\}\\rfloor) - \\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)\\)  
同时，一样的套路，\\(\\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)\\)，类似地，定义为：  
\\\[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 i\}\\sum\_\{e=1\}^\{prime\[j\]^\{e+1\} \\le i\}(f(prime\[j\]^e)\*S(\\lfloor\\frac\{n\}\{prime\[j\]^e\}\\rfloor, j+1) + f(prime\[j\]^\{e+1\}))\\\]

看起来有点长，其实因为合数也分为两种：一种是形如 p^k (k > 1) 式的，另一种则包含两种以上的质因子。  
因为我们 f(1) 是特殊计算的，所以这种分类是必要的。

最后就可以得到 \\(S(i, n) = A + B\\)。计算合数的复杂度就是 O(玄学)。

[一个例题。][Link 3]  
可以发现除了 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 介绍的。。。

基本不会，就咕了。

转载于:https://www.cnblogs.com/Tiw-Air-OAO/p/11604767.html

[zsy_blog]: https://www.cnblogs.com/zhoushuyu/p/8275530.html
[s_orz]: https://blog.csdn.net/skywalkert/article/details/50500009
[jiry_s_orz]: http://jiruyi910387714.is-programmer.com/posts/195270.html
[fjzzq_s_orz]: https://www.cnblogs.com/zzqsblog/p/5461392.html
[fjzzq_sMin_25_orz]: https://www.cnblogs.com/zzqsblog/p/8302815.html
[Min_25_blog_orz]: https://www.cnblogs.com/GuessYCB/p/10061411.html
[Link 1]: https://www.cnblogs.com/Tiw-Air-OAO/p/11304614.html
[Link 2]: https://www.51nod.com/Challenge/Problem.html#problemId=1239
[Link 3]: https://loj.ac/problem/6053