最近突然做到一些求积性函数前缀和的题,用到了各种筛,有一题用到 min_25 筛法,于是好好学习了一波,运用极不熟练。后来又遇到一道杜教筛的题,结果发现自己连 phi(x) 前缀和都不会推了?吓得我赶紧复习 + 写博客。
前置技能
常见(完全)积性函数;整除分块; dirichlet 卷积;埃氏筛
这里还是简单介绍一下 dirichlet 卷积:
我们用 * 表示两个函数的 dirichlet 卷积。
定义函数 f 和 g 的 dirichlet 卷积形式为 $(f*g)(n)=sum_{d|n} f(d)gleft (frac{n}{d} ight )$ 。
接下来,设 $f$ 为积性函数,我们要求 $sum_{i=1}^{n} f(i)$ , n 是 $10^9sim 10^{12}$ 左右的级别。本文忽略复杂度分析(需用到积分相关知识,本人太弱不会)。
杜教筛
(先扔链接:https://www.cnblogs.com/zzqsblog/p/5461392.html)
适用范围:存在另两个积性函数 $g,h$ ,且满足 $f*g=h$ , $g,h$ 函数的前缀和比较容易求。
先来举几个例子:
- 比如说要求 $mu$ 的前缀和,我们知道 $mu * 1=e$ ( $e(i)=[i=1]$ ), $e$ 的前缀和显然恒为 $[nge 1]$ , $1$ 的前缀和也很好求,所以可以使用杜教筛;
- 又比如说要求 $varphi$ 的前缀和,我们知道 $varphi * 1=id$ , $1,id$ 的前缀和都很好求 ,所以也可以使用杜教筛。
那么具体怎么求呢?
设 $S(n)=sum_{i=1}^{n} f(i)$ ,那么有
$$egin{aligned} sum_{i=1}^{n} (f*g)(i) &= sum_{i=1}^{n} sum_{d|i} f(d)gleft (frac{i}{d} ight ) \ S(n)=sum_{i=1}^n f(i) &= sum_{i=1}^{n} (f*g)(i) - sum_{i=1}^{n} sum_{d|i,d<i} f(d)gleft (frac{i}{d} ight ) \ &=sum_{i=1}^n (f*g)(i) - sum_{i=2}^n g(i) sum_{d=1}^{left lfloor frac{n}{i} ight floor} f(d)\ &=sum_{i=1}^n (f*g)(i) - sum_{i=2}^n g(i)Sleft (left lfloor frac{n}{i} ight floor ight ) end{aligned}$$
可以发现我们得到了一个 $S(n)$ 的子问题,可以递归下去求解。
至于边界情况,考虑如果 $nleq B$ 的时候我们直接用线性筛线性求解 $S(n)$ ,当 B 取 $n^{frac 23}$ 的时候复杂度最优,大概为 $mathcal{O}(n^{frac23})$ 。足以解决这个问题。
贴代码 (一百年前写的)
1 #include<bits/stdc++.h> 2 #define rep(i,x,y) for (int i=(x);i<=(y);i++) 3 #define ll long long 4 using namespace std; 5 const int N=1e6+5; 6 int p[N/10],cnt,vis[N],mu[N]; map<ll,ll> mp; 7 void init(int n){ 8 mu[1]=1; 9 rep (i,2,n){ 10 if (!vis[i]) p[++cnt]=i,mu[i]=-1; 11 for (int j=1;j<=cnt&&i*p[j]<=n;j++){ 12 vis[i*p[j]]=1; 13 if (i%p[j]==0) break; 14 mu[i*p[j]]=-mu[i]; 15 } 16 } 17 rep (i,2,n) mu[i]+=mu[i-1]; 18 } 19 ll S(ll n){ 20 if (n<N) return mu[n]; 21 if (mp.count(n)) return mp[n]; 22 ll ret=1; 23 for (ll i=2,last;i<=n;i=last+1){ 24 last=n/(n/i); 25 ret-=(last-i+1)*S(n/i); 26 } 27 return mp[n]=ret; 28 } 29 int main(){ 30 init(N-1); 31 ll a,b; cin>>a>>b; cout<<S(b)-S(a-1); 32 return 0; 33 }
min_25 筛
(先扔链接:https://www.cnblogs.com/cjyyb/p/9185093.html)
一些吐槽:这个杜教筛真的是筛??好吧其实我也不知道它为什么叫杜教筛,这确实只是个递归求和的方法。( yhx : 应该叫杜氏求和法)
不管这些,下面介绍一种真正的筛—— min_25 筛。它本质和洲阁筛相同,但小一倍常数。众所周知数论题非常容易被卡常数,所以学好 min_25 筛还是很重要的 ~
适用范围:函数 f 需要满足 $f(p)$ 是一个低阶多项式, $f(p^c)$ 有较为容易的求法,最好是 $mathcal{O}(1)$ 。这里 p 代表任意素数。
首先 min_25 筛和洲阁筛都是基于埃氏筛的。
要求 $sum_{i=1}^{n} f(i)$ ,我们把 $1sim n$ 中的素数和合数分开考虑。
假设素数集合为 $P$ , $P_i$ 为第 i 个素数, $minp(i)$ 表示 i 的最小素因子。
素数
以下假设 $f(p)$ 为完全积性函数(但 $f(i)$ 不一定。原因下面会说)。
要求 $sum_{i=1}^n [i in P]f(i)$ 。
为了方便计算,这里所有数的函数值都定义为 $f(p)$ ,即都和 $f(p)$ 使用一样的表达式。
设 $g(n,j)=sum_{i=1}^n [iin P or minp(i)>P_j]f(i)$ ,即筛到第 j 个素数为止剩下的数和所有素数的函数值之和。
考虑 $g(n,j)$ 的转移:
- $P_{j}^2>n$ :显然 $g(n,j)=g(n,j-1)$ 。
- $P_{j}^2leq n$ :此时需要减去一些贡献。这个贡献怎么求呢?考虑筛 $P_j$ 的过程:找到剩下的数中所有 $P_j$ 的倍数删去,减去它们的函数值。由于 $f$ 为完全积性函数(上面说的原因在这里,如果不是完全积性就无法提出 $f(P_j)$ ,因而无法进行下面的变换),可以提出 $f(P_j)$ ,即 $f(P_j) gleft (left lfloor frac{n}{P_j} ight floor,j-1 ight )$ 。注意到若直接减去这个值,我们会把质数部分减去,而事实上我们需要把质数保留,因此这个贡献实际上是 $$f(P_j) [gleft (left lfloor frac{n}{P_j} ight floor,j-1 ight ) - sum_{i=1}^{j-1} f(P_i)]\ =f(P_j) [gleft (left lfloor frac{n}{P_j} ight floor,j-1 ight ) - g(P_j-1,j-1)]$$
- 故 $$g(n,j)=egin{cases} g(n,j-1) & P_j^2>n \ g(n,j-1)-f(P_j) [gleft (lfloor frac{n}{P_j} floor,j-1 ight) - g(P_j-1,j-1)] & P_j^2leq n end{cases}$$
那么 $sum_{i=1}^n [i in P]f(i)=g(n,|P|)$ ,这里 $P$ 是前 $sqrt{n}$ 个素数的集合。
【稍微解释一下这个过程】
由于我们最终用到的函数值只有 $g(n,|P|)$ ,即我们这步求的是素数的函数值,按照我们的筛法,实际上是先将所有数都按照 $f(p)$ 的式子计算函数值,然后把合数筛掉,最后只剩下素数,那么就能得到正确的函数值之和。但是中间筛的过程中,如果 $f(x)$ 的表达式和 $f(p)$ 不一样的话, $g$ 的值并不是真实的函数值。所以注意这里只求出了素数函数值之和。
合数
这里没有素数那栏里的若干限制。
设 $S(n,j)=sum_{i=1}^n [minp(i)geq P_j]f(i)$ 。和 $g$ 类似的定义,但注意这里取到了等号,且不包含素数。
$S$ 的值分 2 部分计算:质数,合数。 1 的函数值放到最后计算。
质数的函数值即 $g(n,|P|)-sum_{i=1}^{j-1} f(P_i)$ 。至于合数,由于 $f(x)$ 是积性函数,那么可以把最小质因子全部提出来(它和剩下部分互质),所以枚举它的最小质因子 $P_k$ 及幂次 $e$ ,得 $$S(n,j)=g(n,|P|)-sum_{i=1}^{j-1}f(P_i)+sum_{kgeq j}sum_{e>0}^{P_k^{e+1}leq n}(f(P_k^e)Sleft (left lfloor frac{n}{P_k^e} ight floor,k+1 ight )+f(P_k^{e+1}))$$
后面的 $f(P_k^{e+1})$ 是单独处理的特殊情况,由于前面无法统计到。
这样总复杂度为 $mathcal{O}(frac{n^{frac34}}{log n})$ 。
【为什么要求 $f(p)$ 为低阶多项式呢】
我们在处理素数部分的时候说过假设 $f(p)$ 为完全积性函数。那么低阶多项式一定能拆成 $sum_{i} a_ix^i$ ,这里每一项都是完全积性的,所以只需拆项使用上述方法求解即可。复杂度 $ imes 阶数$ 。
栗子
这题的函数比较奇怪,但好在它是积性的; $f(p)$ 除 $p=2$ 时均为 $p-1$ ,是一个低阶多项式; $f(p^c)$ 可以 $mathcal{O}(1)$ 求解。所以可以用 min_25 筛求解。
把 $f(p)$ 拆掉, $p=2$ 的时候特判一下。别的直接套用上面的公式即可。
code:
1 #include<bits/stdc++.h> 2 #define rep(i,x,y) for (int i=(x);i<=(y);i++) 3 #define ll long long 4 5 using namespace std; 6 7 const int N=2e5+10,mod=1e9+7; 8 ll n,w[N]; int m,sq,id1[N],id2[N],h[N],g[N],p[N],cnt,pre[N],vis[N]; 9 10 inline int get_id(ll x){return x<=sq?id1[x]:id2[n/x];} 11 12 void sieve(int n){ 13 rep (i,2,n){ 14 if (!vis[i]) p[++cnt]=i,pre[cnt]=(pre[cnt-1]+i)%mod; 15 for (int j=1;j<=cnt&&i*p[j]<=n;j++){ 16 vis[i*p[j]]=1; 17 if (i%p[j]==0) break; 18 } 19 } 20 } 21 22 int S(ll n,int i){ 23 if (n<=1||p[i]>n) return 0; 24 int k=get_id(n); 25 int res=((ll)g[k]-h[k]-pre[i-1]+i-1)%mod; res=(res+mod)%mod; 26 if (i==1) res+=2; 27 for (int j=i;j<=cnt&&(ll)p[j]*p[j]<=n;j++){ 28 ll pw=p[j]; 29 for (int e=1;pw<=n/p[j];e++) 30 res=(res+((ll)(p[j]^e)*S(n/pw,j+1)%mod+(p[j]^(e+1)))%mod)%mod,pw*=p[j]; 31 } 32 return res; 33 } 34 35 int main(){ 36 scanf("%lld",&n); sq=sqrtl(n); sieve(sq); 37 for (ll i=1,j;i<=n;i=j+1){ 38 ll x=n/i; j=n/x; w[++m]=x; 39 if (x<=sq) id1[x]=m; else id2[j]=m; 40 h[m]=(x-1)%mod,g[m]=(x+2)%mod*((x-1)%mod)%mod; 41 (g[m]&1)?g[m]+=mod:0; g[m]>>=1; 42 } 43 rep (j,1,cnt) 44 for (int i=1;i<=m&&(ll)p[j]*p[j]<=w[i];i++){ 45 int k=get_id(w[i]/p[j]); 46 h[i]=((ll)h[i]-h[k]+j-1)%mod,h[i]=(h[i]+mod)%mod; 47 g[i]=((ll)g[i]+mod-(ll)p[j]*(g[k]+mod-pre[j-1])%mod)%mod; 48 } 49 printf("%d ",S(n,1)+1); 50 return 0; 51 }