Powerful Number 筛学习笔记
用途
\(Powerful\ number\) 筛可以用来求出一类积性函数的前缀和，最快可以达到根号复杂度。
实现
\(Powerful\ number\) 的定义是每个质因子次数都 \(\ge 2\) 的数。
有如下的性质：
\(1\)、一个 \(Powerful\ number\) 一定可以表示为 \(a^2b^3\) 的形式。
\(2\)、\(n\) 以内的 \(Powerful\ number\) 个数是 \(O(\sqrt n)\) 级别的。
所以找 \(Powerful\ number\) 可以直接暴力 \(dfs\)。
如果要求的函数是 \(f\)，那么我们需要找到一个积性函数 \(g\)，使得 \(f\) 和 \(g\) 在质数处的取值相同。
同时还要找到一个积性函数 \(h\)，使得 \(f=g*h\)。
根据狄利克雷卷积的定义
\(f(p)=g(p)h(1)+g(1)h(p)=g(p)+h(p)=f(p)+h(p)\)。
所以 \(h(p)=0\)，因为 \(h\) 是一个积性函数，所以所有非 \(Powerful\ number\) 在 \(h\) 函数中的取值都是 \(0\)。
\(\sum_{i=1}^nf(i)=\sum_{i=1}^n\sum_{d|i}h(d)g(\frac{i}{d})=\sum_{d=1}^nh(d)\sum_{j=1}^{\frac{n}{d}}g(j)\)。
因为 \(h\) 只在 \(Powerful\ number\) 处有值，所以我们只需要求出 \(\sqrt{n}\) 个 \(g\) 函数的前缀和即可。
例题
题目描述
给定一个积性函数 \(f\)，满足 \(f(1)=1\)，并且对于任意质数 \(p\) 和正整数 \(e\)，都有 \(f(p^e)=p^k\)，\(k\) 为给定的数，\(n \leq 10^{13},k \leq 20\)。
分析
构造积性函数 \(g\)，满足对于任意 \(x\)，都有 \(g(x)=x^k\)。
\(f(p^2)=g(p^2)h(1)+g(p)h(p)+g(1)h(p^2)\)。
那么 \(p^{k}=p^{2k}+h(p^2)\)，\(h(p^2)=p^{k}-p^{2k}\)。
\(f(p^3)=g(p^3)h(1)+g(p^2)h(p)+g(p)h(p^2)+g(1)h(p^3)\)。
\(p^{k}=p^{3k}+p^k(p^{k}-p^{2k})+h(p^3)\)，\(h(p^3)=p^{k}-p^{2k}\)。
多算几项就会发现 \(h(p^e)=p^{k}-p^{2k},e \ge 2\),
直接做就行了。
代码
#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cmath>
#include<map>
#define rg register
template<typename T>void read(rg T& x){
x=0;rg int fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
x*=fh;
}
const int maxn=1e7+5,maxm=35,mod=1e9+7;
inline int addmod(rg int now1,rg int now2){
return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
return now1*=now2,now1>=mod?now1%mod:now1;
}
inline int ksm(rg int ds,rg int zs){
rg int nans=1;
while(zs){
if(zs&1) nans=mulmod(nans,ds);
ds=mulmod(ds,ds);
zs>>=1;
}
return nans;
}
int pri[maxn],k,sqr,b[maxm],c[maxm][maxm],ny[maxm],tot,val[maxn],mi[maxn],ans;
bool not_pri[maxn];
long long n,w[maxn];
void pre(){
for(rg int i=0;i<maxm;i++) c[i][0]=1;
for(rg int i=1;i<maxm;i++){
for(rg int j=1;j<=i;j++){
c[i][j]=addmod(c[i-1][j-1],c[i-1][j]);
}
}
ny[1]=1;
for(rg int i=2;i<maxm;i++) ny[i]=mulmod(mod-mod/i,ny[mod%i]);
b[0]=1;
for(rg int i=1;i<=20;i++){
for(rg int j=0;j<=i-1;j++){
b[i]=addmod(b[i],mulmod(c[i+1][j],b[j]));
}
b[i]=delmod(0,mulmod(b[i],ny[i+1]));
}
}
std::map<int,int> mp;
int getsum(rg int val){
if(mp.find(val)!=mp.end()) return mp[val];
val++;
rg int nans=0,tmp=val;
for(rg int i=k;i>=0;i--){
nans=addmod(nans,mulmod(c[k+1][i],mulmod(b[i],tmp)));
tmp=mulmod(tmp,val);
}
nans=mulmod(nans,ny[k+1]);
return mp[val-1]=nans;
}
void xxs(rg int mmax){
not_pri[0]=not_pri[1]=1;
for(rg int i=2;i<=mmax;i++){
if(!not_pri[i]){
pri[++pri[0]]=i;
mi[pri[0]]=delmod(ksm(i,k),ksm(i,k<<1));
}
for(rg int j=1;j<=pri[0] && i*pri[j]<=mmax;j++){
not_pri[i*pri[j]]=1;
if(i%pri[j]==0) break;
}
}
}
void dfs(rg int now,rg long long nw,rg int nval){
w[++tot]=nw,val[tot]=nval;
for(rg int i=now;i<=pri[0] && nw<=n/pri[i]/pri[i];i++){
rg long long tmp=1LL*nw*pri[i];
for(;tmp<=n/pri[i];tmp*=pri[i]) dfs(i+1,tmp*pri[i],mulmod(nval,mi[i]));
}
}
int main(){
pre();
read(n),read(k);
sqr=sqrt(n);
xxs(sqr+5);
dfs(1,1,1);
for(rg int i=1;i<=tot;i++){
ans=addmod(ans,mulmod(val[i],getsum((n/w[i])%mod)));
}
printf("%d\n",ans);
return 0;
}