Min25筛法学习小记
问题引入
已知积性函数\(f(x)\)
要求其前\(n\)项和(\(n\leq 10^{10}\))
其中对于所有质数\(p\)
满足\(f(p)=a_0+a_1p+a_2p^2……\)
也就是可以表示成一个低阶多项式
且\(f(p^k)\)可以快速计算
I
定义\(lst(x)\)表示\(x\)的最小质因子，\(p_k\)表示第\(k\)个质数，特别的,\(p_0=1\)
我们把答案拆成质数,合数,与1
也就是
\(·\) 1的贡献是1
\(.\) 质数的贡献：
\(·\) 合数的贡献
我们可以枚举合数的最小质因子和质因子的指数
内层换成枚举\(p^c\)的倍数\(i\)，则\(lst(d)=p\)等价于\(lst(i)>p\)
注意因为枚举倍数，所以都会有\(f(p^c)\)，因为是积性函数，所以可以提出来
但是还有问题
因为\(p^c\)本身也有贡献，如果\(c!=1\)的话
因为\(c==1\)也就是质数，我们已经贡献了
最终的式子是
中间一坨看着很不爽
我们设\(S(n,k)=\sum_{i=1}^n[lst(i)>p_k]f(i)\)
那么答案就是
考虑用这个式子代替上式
考虑如何推\(S(n,k)\)
和上面一样，考虑质数和合数的贡献
质数的贡献
合数的贡献
总结一下就是
因为我们最初给出了限制:\(f(p_k)\)可以快速计算，所以这个可以忽略
注意到以\(p_k\)为\(lst(x)\)的最小合数\(x\)是\(p_k^2\)
所以后边的\(i\)只需枚举到\(\sqrt n\),总体的这个\(k\)也是只到\(\sqrt n\)
因此后半部分只需要筛出来小于等于\(\sqrt(n)\)的质数就可以递归进行了
问题在于前半部分
II
我们现在要求
根据上边可以知道这里的\(k\)是小于等于\(\sqrt n\)的
我们可以预处理出\(sum(k)\)表示前\(k\)个质数的\(f(p)\)之和
那么上式等价于
也就是我们只要能知道所有质数的函数就行了
又因为我们一开始限制，函数在质数处的取值是低阶多项式
所以我们拆成一些单项式，分别求和就行了
也就是说我们只需对形如\(g(p)=p^k\)在质数处求和就行了
注意到这是一个完全积性函数哦
与\(I\)类似，我们设
通俗来说就是用前\(k\)个质数做线性筛剩下的数的函数之和
那么所有质数的\(g(p)\)之和，就是\(G(n,r)\)
其中\(r\)是满足\(p_r^2\leq n\)的最大\(r\)，也就是我们筛出来的最后一个质数
考虑怎么递推\(G(n,k)\)
它肯定是在\(G(n,k-1)\)的基础上多筛出了一下数
那么筛出了那些数呢，就是以\(lst(x)= p_k\)的\(x\)
我们依旧可以提出一个\(p_k\)
直接枚举倍数
注意，此时后面可能还有\(p_k\)，但是因为\(g(p_k)\)是完全积性函数，所以没事
\(\geq p_k\to>p_{k-1}\)
似乎后面的东西和我们的\(G\)有点像，观察其实就是少算了那些\(\leq p_{k-1}\)的质数
那么直接用\(G\)减掉就行了
以\(k\)为阶段，递推式就出来了
不过我们似乎存不下这么多值
观察到因为\(k\leq \sqrt n\)
所以最后一项可以预处理
问题就在于\(G(\lfloor \frac{n}{p_k}\rfloor,k-1)\)
但是我们又发现
也就是说，我们访问到的只会是形如\(\lfloor \frac{n}{d}\rfloor\)
这样子的位置的值
根据整除分块的理论，这种取值只会有\(\sqrt n\)个
就开的下了
不过因为每个数都很大
所以要么可以离散化，要么可以用这种方法
设\(x\)为某一个可能的取值的位置
若\(x\leq \sqrt n\)，直接记录\(x\)
否则，记录\(\lfloor \frac{n}{x}\rfloor\)
根据某些证明，该算法的时间复杂度是\(O(\frac{n^{\frac{3}{4}}}{\log n})\)
但是实际的效率是很高的
III
应用
P5325 【模板】Min_25筛
直接把在质数处的表达式告诉我们了
我们直接拆成\(g1(p)=p^2\),和\(g2(p)=p\)
用\(g1(p)-g2(p)\)就行了
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N = 1e6+7;
const LL mod = 1e9+7;
const LL inv2=5e8+4,inv3=333333336;
LL prime[N],tot=0;
LL sum[N],Sum[N];
LL Id1[N],Id2[N];
LL v[N];
void init(LL n)
{
v[1]=0;
for(LL i=2;i<=n;i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
sum[tot]=(sum[tot-1]+i)%mod;
Sum[tot]=(Sum[tot-1]+i*i%mod)%mod;
}
for(LL j=1;j<=tot;j++)
{
if(i*prime[j]>n||prime[j]>v[i]) continue;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0) break;
}
}
}
LL B;
LL num[N],cnt=0;
LL G1[N],G2[N];
LL W=0;
LL S(LL n,LL m)
{
if(prime[m]>=n) return 0;
LL id=(n<=B?Id1[n]:Id2[W/n]);
LL ans=(G2[id]-G1[id]+mod)%mod-(Sum[m]-sum[m]+mod)%mod;
ans=(ans+mod)%mod;
// cout<<n<<' '<<m<<' '<<id<<' '<<G2[id]-G1[id]<<endl;
for(LL k=m+1;k<=tot&&prime[k]*prime[k]<=n;k++)
{
LL cur=prime[k];
for(LL c=1;cur<=n;c++,cur=cur*prime[k])
{
LL x=cur%mod;
LL v=x*(x-1)%mod;
ans=(ans+v*(S(n/cur,k)+(c!=1)%mod)%mod)%mod;
}
}
return ans;
}
int main()
{
LL n;
cin>>n;
B=sqrt(n);
init(B);
W=n;
LL l=1,r;
// cout<<3*inv3%mod;
for(;l<=n;l=r+1)
{
r=(n/(n/l));
num[++cnt]=(n/l);
LL tmp=num[cnt]%mod;
G1[cnt]=tmp*(tmp+1)%mod*inv2%mod-1;
G2[cnt]=tmp*(tmp+1)%mod*inv2%mod*(2*tmp%mod+1)%mod*inv3%mod-1;
if(n/l<=B) Id1[n/l]=cnt;
else Id2[n/(n/l)]=cnt;
}
for(LL i=1;i<=tot;i++)
{
for(LL j=1;j<=cnt&&prime[i]*prime[i]<=num[j];j++)
{
LL v=num[j]/prime[i],id;
if(v<=B) id=Id1[v];
else id=Id2[n/v];
// cout<<j<<' '<<id<<' '<<G2[id]<<endl;
G1[j]-=prime[i]*(G1[id]-sum[i-1]+mod)%mod;
G2[j]-=prime[i]*prime[i]%mod*(G2[id]-Sum[i-1]+mod)%mod;
G1[j]%=mod;
G2[j]%=mod;
if(G1[j]<0) G1[j]+=mod;
if(G2[j]<0) G2[j]+=mod;
}
}
cout<<(S(n,0)+1)%mod;
return 0;
}
同样也可解决许多杜教筛可以解决的问题
筛\(\mu\)
\(\mu(p)=-1\)
直接算就行了
筛\(\varphi\)
\(\varphi(p)=p-1\)
直接算就行了
筛除数函数\(\sigma_k\)
\(\sigma_k(p)=p^k\)
应该也可以筛一些奇奇怪怪的东西
比如
浙公网安备 33010602011771号