min_25筛学习笔记
min_25筛是一个能快速求解积性函数前缀和的东西。
要保证 $f(p)(p\text{是质数})$ 是个关于 $p$ 的多项式(次数也不要太高),并且 $f(p^k)$ 能快速计算。
以下以洛谷的模板为例:($f(p^k)=(p^k)^2-p^k(p\text{是质数})$,求前 $N(N\le 10^{10})$ 项的和 $\bmod 10^9+7$ 的值)
也就是 $\sum\limits^N_{i=1}f(i)$。
我们考虑分成质数和合数计算,同时枚举合数的最小质因子及其次数(因为是积性函数,所以直接提出来)。
答案为 $\sum\limits^N_{p=1}f(p)[p\in prime]+\sum f(p^e)[p\in prime,p^e\le N](\sum\limits_{i=1}^{\lfloor\frac{N}{p^e}\rfloor}f(i)[i\text{的最小质因子}>p])$
首先考虑质数部分:
我们对每个次数不同的项分别计算。假设我们正在计算 $k$ 次方项。
我们不知道为什么但就是设 $g(n,j)$ 为 $\sum\limits^n_{i=1}i^k[i\text{是质数或者}i\text{的最小质因子}>p_j]$。
其中 $p_j$ 是第 $j$ 小的质数。
那么有初始状态 $g(n,0)=\sum\limits^n_{i=2}i^k$(不把 $1$ 算进去)。
转移,发现不满足 $g(n,j)$ 的限制且满足 $g(n,j-1)$ 的限制的数就是最小质因子恰为 $p_j$ 的合数。
首先发现 $p_j^2>n$ 时这样的数不存在,此时 $g(n,j)=g(n,j-1)$。
否则由于 $k$ 次方是完全积性函数,所以可以把 $p_j^k$ 提出来,那么这些数的 $k$ 次方和就是 $p_j^k(g(\lfloor\dfrac{n}{p_j}\rfloor,j)-g(p_{j-1},j-1))$。
为什么前面的 $g$ 第二维是 $j$ 呢?因为我们只提出了一个 $p_j$,剩下可能还有。
为什么要减掉后面的 $g$ 呢?因为后面是 $<p_j$ 的质数的 $k$ 次方和,如果把这一部分算上,那么就相当于整个答案多算了一些最小质因子小于 $p_j$ 的数。所以要减掉。(其实这一块就是质数 $k$ 次方的前缀和,可以预处理。以后记作 $s_k(j-1)$)
最后用 $g(n,j-1)$ 减掉不满足 $g(n,j)$ 的限制且满足 $g(n,j-1)$ 的限制的数即可。
方程出来了:
$$g(n,j)=\begin{cases}g(n,j-1)&p_j^2>n\\g(n,j-1)-p_j^k(g(\lfloor\dfrac{n}{p_j}\rfloor,j)-s_k(j-1))&p_j^2\le n\end{cases}$$
这个可以上滚动数组滚掉 $j$,同时对于 $p_j^2>n$ 的不用动就行了,时间复杂度也降到了……等等,这玩意复杂度明明比暴力还高好吗!
但是我们可以发现,一个 $g(n,j)$ 可以由若干个 $g(\lfloor\dfrac{n}{x}\rfloor,k)$ 表示。
这启示我们把所有 $\lfloor\dfrac{N}{n}\rfloor$ 相同的 $n$ 归为一类一起处理。具体实现,把 $g$ 的第一维下标变为 $\lfloor\dfrac{N}{n}\rfloor$。然而这样还是可能太大,于是可以对 $\le\sqrt{n}$ 和 $>\sqrt{n}$ 的开两种编号。
这一部分时间复杂度约为 $O(\frac{N^{3/4}}{\log N})$。
怎么算的?我们发现类似埃氏筛,一个数 $x$ 被更新过 $\frac{\sqrt{x}}{\log\sqrt{x}}$ 次。
那么复杂度 $O(\sum\limits_{i=1}^{\sqrt{N}}\frac{\sqrt{i}}{\log\sqrt{i}}+\sum\limits_{i=1}^{\sqrt{N}}\frac{\sqrt{\frac{N}{i}}}{\log\sqrt{i}})$。经过积分爆算,发现略小于 $O(\frac{N^{3/4}}{\log N})$。
先放一下模板的代码实现:
inline int id(ll x){return x<=sq?id1[x]:id2[n/x];} void init(){ sq=sqrt(n); FOR(i,2,sq){ if(!vis[i]){ pri[++pl]=i; s1[pl]=(s1[pl-1]+i)%mod; s2[pl]=(s2[pl-1]+(ll)i*i)%mod; //质数一次项和二次项前缀和 } FOR(j,1,pl){ if(i*pri[j]>sq) break; vis[i*pri[j]]=true; if(i%pri[j]==0) break; } } for(ll l=1,r;l<=n;l=r+1){ r=n/(n/l); w[++tot]=n/l; if(n/l<=sq) id1[n/l]=tot; else id2[n/(n/l)]=tot; //设置编号 int x=w[tot]%mod; g1[tot]=(1ll*x*(x+1)%mod*inv2%mod-1+mod)%mod; g2[tot]=(1ll*x*(x+1)%mod*(2*x+1)%mod*inv6%mod-1+mod)%mod; //2到x的和 } } void calc_g(){ FOR(i,1,pl){ //遍历质数(上文中的j) FOR(j,1,tot){ //遍历上文中的n if((ll)pri[i]*pri[i]>w[j]) break; //只管p_j^2<=n的 g1[j]=(g1[j]-1ll*pri[i]*(g1[id(w[j]/pri[i])]-s1[i-1]+mod)%mod+mod)%mod; g2[j]=(g2[j]-1ll*pri[i]*pri[i]%mod*(g2[id(w[j]/pri[i])]-s2[i-1]+mod)%mod+mod)%mod; //上面这两行就死磕吧……记得g1,g2和w下标的含义就好 } } }
再看整个式子:
我们吸取了经验知道要绞尽脑汁地设 $S(n,j)$ 表示 $\sum\limits_{i=2}^n f(i)[i\text{的最小质因子}>p_j]$。注意这里是 $f$ 而不是 $k$ 次方,下界是 $2$。
那么要求即为 $S(N,0)+f(1)$。边界 $S(n,j)=0(p_j>n)$。
转移,对答案有贡献的质数必定大于 $p_j$,那么贡献就是 $g(n,pl)-\sum s_k(j)$。($pl$ 是 $\le\sqrt{N}$ 的质数个数)
按照套路再枚举合数的最小质因子和次数,就是 $\sum f(p^e)[p\in prime,p>p_j,p^e\le n](S(\lfloor\dfrac{n}{p^e}\rfloor)+[e\ne 1])$。
(为何后面要有 $[e\ne 1]$ 呢?因为 $e=1$ 时,这个值在前面算质数时已经算过,而如果 $e\ne 1$,那么前面算质数时没算过,而因为 $S$ 下界为 $2$ 后面也不会算到,所以要加上 $f(p^e)$)。
即使不记忆化,这一部分复杂度也是 $O(\frac{N^{3/4}}{\log N})$。(这个复杂度我是真不会分析了)
贴上整题代码:
#include<bits/stdc++.h> using namespace std; typedef long long ll; const int maxn=1000100,mod=1000000007,inv2=500000004,inv6=166666668; #define FOR(i,a,b) for(int i=(a);i<=(b);i++) #define ROF(i,a,b) for(int i=(a);i>=(b);i--) #define MEM(x,v) memset(x,v,sizeof(x)) int sq,pri[maxn],pl,tot,g1[maxn],g2[maxn],id1[maxn],id2[maxn],s1[maxn],s2[maxn]; bool vis[maxn]; ll n,w[maxn]; inline int add(int a,int b){return a+b<mod?a+b:a+b-mod;} inline int sub(int a,int b){return a<b?a-b+mod:a-b;} inline int mul(int a,int b){return (ll)a*b%mod;} inline int id(ll x){return x<=sq?id1[x]:id2[n/x];} void init(){ sq=sqrt(n); FOR(i,2,sq){ if(!vis[i]) pri[++pl]=i,s1[pl]=add(s1[pl-1],i),s2[pl]=add(s2[pl-1],mul(i,i)); FOR(j,1,pl){ if(i*pri[j]>sq) break; vis[i*pri[j]]=true; if(i%pri[j]==0) break; } } for(ll l=1,r;l<=n;l=r+1){ r=n/(n/l); w[++tot]=n/l; if(n/l<=sq) id1[n/l]=tot; else id2[n/(n/l)]=tot; g1[tot]=sub(mul(mul(w[tot]%mod,(w[tot]+1)%mod),inv2),1); g2[tot]=sub(mul(mul(mul(w[tot]%mod,(w[tot]+1)%mod),(2*w[tot]+1)%mod),inv6),1); } } void calc_g(){ FOR(i,1,pl){ FOR(j,1,tot){ if((ll)pri[i]*pri[i]>w[j]) break; g1[j]=sub(g1[j],mul(pri[i],sub(g1[id(w[j]/pri[i])],s1[i-1]))); g2[j]=sub(g2[j],mul(mul(pri[i],pri[i]),sub(g2[id(w[j]/pri[i])],s2[i-1]))); } } } int solve(ll nn,int x){ //S(nn,x) if(pri[x]>=nn) return 0; int ans=sub(sub(g2[id(nn)],s2[x]),sub(g1[id(nn)],s1[x])); //题目中f为二次项减一次项 FOR(i,x+1,pl){ //p if((ll)pri[i]*pri[i]>nn) break; ll pro=1; FOR(j,1,50){ //e pro*=pri[i]; //p^e if(pro>nn) break; ans=add(ans,mul(mul(pro%mod,(pro-1)%mod),add(solve(nn/pro,i),j!=1))); } } return ans; } int main(){ scanf("%lld",&n); init(); calc_g(); printf("%d\n",add(solve(n,0),1)); }
一些题目:
LOJ6235 区间素数个数 题解
LOJ6053 简单的函数 题解