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}^nsum_{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
",ans);
return 0;
}