powerful number求积性函数前缀和
powerful number即每个质因子出现次数都≥2的数;一个非常重要的性质是≤n的powerful number是O(√n)级别的。如果想要求出所有≤n的powerful number的话,只需要写一个类似于min_25筛第二部分的爆搜即可。大概就是记录一下当前要找的powerful number的最小质数是什么,之后枚举质数次幂即可。
这个东西比较厉害的一点是可以用来优化求积性函数前缀和的过程。
设F(pq)=pk,且F(x)是一个积性函数,求∑ni=1F(i)。
考虑找两个积性函数G(x),H(x),满足F(x)=G(x)H(x),同时对于任意质数p,均满足F(p)=G(p)。不难发现,F(p)=G(1)H(p)+H(1)G(p),很显然H(1)=G(1)=1,于是F(p)=H(p)+G(p),又因为G(p)=F(p),所以显然有H(p)=0。
所以如果x满足H(x)≠0,那么x一定是一个powerful number。而n∑i=1F(i)=∑i×j≤nH(i)G(j)=n∑i=1H(i)⌊ni⌋∑j=1G(j),于是我们只需要在O(√n)个powerful number处计算即可。
现在的两个问题是求出G(i)的前缀和,以及H(i)的值。
对于这道题,我们令G(x)=xk,那么显然有G(p)=F(p)=pk,于是G(i)的前缀和就是一个自然数幂次方和,可以利用插值等多种多样的方法求出。
对H(i),显然有F(pq)=∑qi=0G(pi)H(pq−i),注意到这实际上是一个多项式卷积的形式,已知F,G的情况下求H就是一个多项式求逆问题,一个O(q2)的暴力即可。我们在搜powerful number的时候实际上在搜质因子组成,于是把各个质因子的函数值乘在一起就好了。
时间复杂度是O(√n)乘上计算G(i)前缀和的复杂度,对于上面的问题,复杂度就是O(k√n)。
上面那道题的代码
#include<cmath>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define re register
#define LL long long
const int maxn=3.3e6+5;
const int mod=1e9+7;
std::vector<int> h[300005];
LL n;int k,Sqr,ans,f[maxn],sz[maxn<<1],vis[maxn<<1],p[maxn>>1];
int A[105],B[105],g[105],ifac[105],suf[105],pre[105];
inline int dqm(int x) {return x<0?x+mod:x;}
inline int qm(int x) {return x>=mod?x-mod:x;}
inline int ksm(int a,int b) {
int S=1;for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)S=1ll*S*a%mod;return S;
}
inline int get(LL nw) {
int id=nw<=Sqr?nw:Sqr+n/nw;
if(vis[id]) return sz[id];
nw%=mod;int tot=0;vis[id]=1;
pre[0]=nw;for(re int i=1;i<=k;i++)pre[i]=1ll*pre[i-1]*dqm(nw-i)%mod;
suf[k+1]=1;for(re int i=k;i;--i)suf[i]=1ll*suf[i+1]*dqm(nw-i)%mod;
for(re int i=1;i<=k;i++) {
int v=1ll*pre[i-1]*suf[i+1]%mod*g[i]%mod*ifac[i]%mod*ifac[k-i]%mod;
if((k-i)&1) v=dqm(-v);tot=qm(tot+v);
}
return sz[id]=tot;
}
void dfs(int x,int v,LL res) {//表示当前要找的powerful number的最小质数次幂>=第x个质数
ans=qm(ans+1ll*v*get(res)%mod);
if(x==p[0]+1)return;if(res/p[x]/p[x]<=0) return;LL R=res/p[x]/p[x];
for(re int i=x;i<=p[0]&&res/p[i]/p[i]>0;++i,R=res/(p[i]?p[i]:1)/(p[i]?p[i]:1))
for(re int j=2;R>0;++j,R/=p[i]) dfs(i+1,1ll*v*h[i][j]%mod,R);
}
int main() {
scanf("%lld%d",&n,&k);Sqr=1+sqrt(n);
for(re int i=2;i<=Sqr;i++) {
if(!f[i]) p[++p[0]]=i;
for(re int j=1;j<=p[0]&&p[j]*i<=Sqr;++j) {
f[p[j]*i]=1;if(i%p[j]==0)break;
}
}
for(re int i=1;i<=p[0];++i) {
int t=0;LL res=n;while(res>=p[i])res/=p[i],++t;
A[0]=B[0]=1;A[1]=ksm(p[i],k);
for(re int j=2;j<=t;++j)A[j]=A[j-1];
for(re int j=1;j<=t;++j)B[j]=1ll*B[j-1]*A[1]%mod;
h[i].push_back(1);
for(re int j=1;j<=t;++j) {
int nw=A[j];
for(re int k=0;k<j;++k)nw=dqm(nw-1ll*h[i][k]*B[j-k]%mod);
h[i].push_back(nw);
}
}++k;
for(re int i=1;i<=k;i++)g[i]=qm(g[i-1]+ksm(i,k-1));ifac[0]=ifac[1]=1;
for(re int i=2;i<=k;i++)ifac[i]=1ll*(mod-mod/i)*ifac[mod%i]%mod;
for(re int i=2;i<=k;i++)ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;
dfs(1,1,n);printf("%d\n",ans);return 0;
}
难点在于构造一个保证G(p)=F(p)积性函数G(i),同时这个函数还能快速求前缀和,所以灵活程度上可能并不如min_25筛。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· go语言实现终端里的倒计时
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 分享 3 个 .NET 开源的文件压缩处理库,助力快速实现文件压缩解压功能!
· Ollama——大语言模型本地部署的极速利器
· DeepSeek如何颠覆传统软件测试?测试工程师会被淘汰吗?