zoukankan      html  css  js  c++  java
  • powerful number求积性函数前缀和

    ( ext{powerful number})即每个质因子出现次数都(geq 2)的数;一个非常重要的性质是(leq n)( ext{powerful number})(O(sqrt{n}))级别的。如果想要求出所有(leq n)( ext{powerful number})的话,只需要写一个类似于( ext{min_25})筛第二部分的爆搜即可。大概就是记录一下当前要找的( ext{powerful number})的最小质数是什么,之后枚举质数次幂即可。

    这个东西比较厉害的一点是可以用来优化求积性函数前缀和的过程。

    (F(p^q)=p^k),且(F(x))是一个积性函数,求(sum_{i=1}^nF(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) eq 0),那么(x)一定是一个( ext{powerful number})。而(displaystyle sum_{i=1}^nF(i)=sum_{i imes jleq n}H(i)G(j)=sum_{i=1}^nH(i)sum_{j=1}^{lfloor frac{n}{i} floor}G(j)),于是我们只需要在(O(sqrt{n}))( ext{powerful number})处计算即可。

    现在的两个问题是求出(G(i))的前缀和,以及(H(i))的值。

    对于这道题,我们令(G(x)=x^k),那么显然有(G(p)=F(p)=p^k),于是(G(i))的前缀和就是一个自然数幂次方和,可以利用插值等多种多样的方法求出。

    (H(i)),显然有(F(p^q)=sum_{i=0}^qG(p^i)H(p^{q-i})),注意到这实际上是一个多项式卷积的形式,已知(F,G)的情况下求(H)就是一个多项式求逆问题,一个(O(q^2))的暴力即可。我们在搜( ext{powerful number})的时候实际上在搜质因子组成,于是把各个质因子的函数值乘在一起就好了。

    时间复杂度是(O(sqrt{n}))乘上计算(G(i))前缀和的复杂度,对于上面的问题,复杂度就是(O(ksqrt{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
    ",ans);return 0;
    }
    

    难点在于构造一个保证(G(p)=F(p))积性函数(G(i)),同时这个函数还能快速求前缀和,所以灵活程度上可能并不如( ext{min_25})筛。

  • 相关阅读:
    SpringBoot Mybatis的驼峰命名
    设计模式之单例模式
    Java调用TSC打印机进行打印
    DOM与Jquery方法对照表(versions:Itcast)
    jQuery选择器过滤器
    算法的力量转自李开复
    源码中的相对路径和绝对路径
    C语言I博客作业03
    firstprogram
    C语言I博客作业02
  • 原文地址:https://www.cnblogs.com/asuldb/p/13067903.html
Copyright © 2011-2022 走看看