Min_25筛目害学
带着理解抄了一遍模板题解来加深印象
求解积性函数前缀和一般可以使用杜教筛。但当这种积性函数不常见且难于卷出已于求前缀和的函数时,杜教筛是乏力的。这时可以考虑使用Min_25筛。
以下内容中约定 \(p\in\mathbb{P}\) ,\(p_i\) 为第 \(i\) 大的质数。
使用Min_25筛的条件:该函数(不妨令其为 \(F\) ) 是积性函数,且 \(F(p)\) 为关于 \(p\) 的低阶多项式, \(F(p^k)\) 可快速求值。
低次多项式的限制是因为Min_25筛主要用来解决求 \(f(p)=p^k\) 的前缀和,把多项式拆成若干单项式,分别求出前缀和后相加就能得到答案。
Min_25筛求解主要分两个步骤。
Step1 分类
我们按质数与合数将 \(\sum_{i=1}^nf(i)\) 分类。
对后一部分枚举最小质因子来满足限制。设 \(lpf_i\) 为 \(i\) 的最小质因子,
于是对两部分分别求解。
Step2 求解
第一部分
线性筛出所有 \(f(p)\) 肯定是不现实的。我们定义 \(g_{n,i}\) 为 \(n\) 以内所有质数与最小质因子大于 \(p_i\) 的数的 \(k\) 次方和,我们要求的东西就是 \(g(x,\sqrt x)\),考虑DP求这个东西。
\(g(n,0)\) 可以 \(O(1)\) 求,考虑如何从 \(g(n,j-1)\) 转移到 \(g(n,j)\) 。这步操作去掉的部分为最小质因子为 \(p_j\) 且不大于 \(n\) 的合数 \(k\) 次方值。
后一个 \(g\) 去掉了所有小于 \(p_j\) 的质数。
实际上我们需要的只是 \(g(x,p_y)\) ,其中 \(p_y\) 为最后一个小于 \(x\) 的质数。称它为 \(pw_x\) 。
求解答案
继续设 \(S_{n,x}\) 为 \([1,n]\) 内最小质因子大于 \(p_x\) 的数函数值之和。继续分类考虑转移。第一部分是大于 \(p_x\) 的质数,可写作 \(pw_n-pw_x\) 。另一部分是剩余合数,可以递归求解。
然后发现我们实际上只需要 \(\frac{n}{x}\) 处的 \(pw\) 值,因为合数的最小质因子不超过根号,所以 \(pw_x\) 都可以直接预处理出来。于是求解 \(g\) 时只需求 \(O(\sqrt n)\) 个值。
细节
离散化时可以开两个数组,分别记录 \(\frac{n}{x}\) 和 \(\frac{n}{\frac{n}{x}}\) 离散化后的位置,存的时候哪个小就存哪个,可以保证数组大小都是 \(O(\sqrt n)\) 的。
1的贡献单算,最后加一。
\(code:\)
P5325[模板]Min_25筛
#include<bits/stdc++.h>
using namespace std;
namespace IO{
#define int LL
typedef long long LL;
int read(){
int x=0,f=0; char ch=getchar();
while(ch<'0'||ch>'9'){ f|=ch=='-'; ch=getchar(); }
while(ch>='0'&&ch<='9'){ x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); }
return f?-x:x;
} char outp[50];
void write(int x,char sp){
int len=0;
if(x<0) x=-x, putchar('-');
do{ outp[len++]=x%10+'0'; x/=10; }while(x);
for(int i=len-1;~i;i--) putchar(outp[i]); putchar(sp);
}
} using namespace IO;
const int NN=300010,mod=1e9+7,inv3=333333336;
int n,sqr;
void pls(int& x,int y){ (x+=y)>=mod?(x-=mod):(x<0?(x+=mod):0); }
namespace Min_25{
int cnt,pri[NN],p1[NN],p2[NN];
int tot,g1[NN],g2[NN],has[NN],pos1[NN],pos2[NN];
bool vis[NN];
void sieve(int n){
for(int i=2;i<=n;i++){
if(!vis[i]){ pri[++cnt]=i; pls(p1[cnt]=p1[cnt-1],i); pls(p2[cnt]=p2[cnt-1],i*i%mod); }
for(int j=1;j<=cnt&&i*pri[j]<=n;j++){
vis[i*pri[j]]=1;
if(i%pri[j]==0) break;
}
}
}
void divblock(){
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
has[++tot]=n/l;
g1[tot]=has[tot]%mod;
pls(g2[tot]=g1[tot]*(g1[tot]+1)/2%mod*(g1[tot]*2+1)%mod*inv3%mod,-1);
pls(g1[tot]=g1[tot]*(g1[tot]+1)/2%mod,-1);
n/l<=sqr?pos1[n/l]=tot:pos2[r]=tot;
}
}
void getg(){
for(int i=1;i<=cnt;i++)
for(int j=1;j<=tot&&pri[i]*pri[i]<=has[j];j++){
int k=has[j]/pri[i]<=sqr?pos1[has[j]/pri[i]]:pos2[n/(has[j]/pri[i])];
pls(g1[j],-pri[i]*(g1[k]-p1[i-1])%mod);
pls(g2[j],-pri[i]*pri[i]%mod*(g2[k]-p2[i-1])%mod);
}
}
int calc(int x,int y,int res=0){
if(pri[y]>=x) return 0;
int k=x<=sqr?pos1[x]:pos2[n/x];
pls(res,(g2[k]-g1[k]+p1[y]-p2[y])%mod);
for(int i=y+1;i<=cnt&&pri[i]*pri[i]<=x;i++)
for(int j=1,p=pri[i];p<=x;j++,p=p*pri[i]){
int pp=p%mod;
pls(res,pp*(pp-1)%mod*(calc(x/p,i)+(j!=1))%mod);
}
return res;
}
} using namespace Min_25;
signed main(){
n=read(); sqr=sqrt(n);
sieve(sqr); divblock(); getg();
write((calc(n,0)+1)%mod,'\n');
return 0;
}