powerful number求积性函数前缀和
类似杜教筛的构造思想。
构造容易求的东西求出f的前缀和
不同之处是,杜教筛构造出两个可以快速求前缀和
但是powerful number(简称pn)构造方法有所不同。
对于这个题
数据范围,min_25筛显然是过不去的。最多跑到1e11
$f(p^e)=p^k$启发
构造$g(x)=x^k$和相似一些
再构造$h(x)$,使得$f=h*g$这里表示狄利克雷卷积
那么$\sum_i f(i)=\sum_i\sum_{d|i}h(i)\times g(i/d)=\sum_d h(d)\sum_{i=1}^{n/d}g(i)$
g的前缀和好算 。拉格朗日差值。
h?
先考虑h是什么
首先,f是积性函数,g也是,不妨考虑能不能构造h也是
要使得对于$f=h*g$
只需要对于$f(p^e)$都是成立的。
首先$h(1)=1$
$f(p)=h(1)\times g(p)+h(p)\times g(1)=p$故$h(p)=0$
每个质因子次数都>=2的才有值,
这些数就叫做powerful number
每个数都可以写成$a^2b^3$的形式。
那么发现,powerful number只有$\sqrt(n)$个
可以暴力枚举!
h有值情况是什么?
$f(p^e)=\sum_{i=0}^{e} h(p^i)\times g(p^{e-i})$
$h(p^e)=f(p^e)-\sum_{i=0}^{e-1} h(p^i)\times g(p^{e-i})$
已经可以前缀和优化求出$h(p^e)$了
根据本题的g的定义,归纳可以得到$h(p^e)=p^k-p^{2k}$
然后可以通过本题。
具体的,预处理$<=sqrt(n)$的质数。
预处理拉格朗日差值需要的东西。
枚举h的话,
法一:枚举$a^2b^3$再进行质因数分解。
法二:枚举每个质因子的次数。从2开始。dfs枚举,或者类似min_25筛的方法枚举也一样
(注意枚举方法:
枚举下一个质因子>=2的位置!不要什么dfs(d,x,h)->dfs(d+1,x,h)表示不选这样。会多余很多不选择的情况。
而枚举下一个,可以O(总个数)枚举所有的。见代码
)
$O(sqrt(n)*k)$
加强一下,
预处理$<=sqrt(n)$的g的前缀和(可以用线性筛)
可以做到:$O(sqrt(n)+sqrt(n)+sqrt(sqrt(n))\times k)$
k=5000,n=1e14也可以3s内跑过!
#include<bits/stdc++.h> #define reg register int #define il inline #define fi first #define se second #define mk(a,b) make_pair(a,b) #define numb (ch^'0') #define pb push_back #define solid const auto & #define enter cout<<endl #define pii pair<int,int> //#define int long long using namespace std; typedef long long ll; template<class T>il void rd(T &x){ char ch;x=0;bool fl=false;while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb);(fl==true)&&(x=-x);} template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');} template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');} template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar('\n');} namespace Modulo{ const int mod=1e9+7; il int ad(int x,int y){return x+y>=mod?x+y-mod:x+y;} il int sub(int x,int y){return ad(x,mod-y);} il int mul(int x,int y){return (ll)x*y%mod;} il void inc(int &x,int y){x=ad(x,y);} il void inc2(int &x,int y){x=mul(x,y);} il int qm(int x,int y=mod-2){int ret=1;while(y){if(y&1) ret=mul(x,ret);x=mul(x,x);y>>=1;}return ret;} template<class ...Args>il int ad(const int a,const int b,const Args &...args) {return ad(ad(a,b),args...);} template<class ...Args>il int mul(const int a,const int b,const Args &...args) {return mul(mul(a,b),args...);} } using namespace Modulo; namespace Miracle{ const int N=1e7+4; const int K=5005; ll n,k; int mi[N]; int pri[N],cnt; int vis[N]; int sum[N]; int sqr; void sieve(ll n){ sum[1]=1; for(reg i=2;i<=n;++i){ if(!vis[i]){ pri[++cnt]=i; sum[i]=mi[cnt]=qm(i,k); } for(reg j=1;j<=cnt;++j){ if(i*pri[j]>n) break; vis[i*pri[j]]=1; sum[i*pri[j]]=mul(sum[i],sum[pri[j]]); if(i%pri[j]==0){ break; } } } for(reg i=1;i<=n;++i){ inc(sum[i],sum[i-1]); } } int Y[K]; int zh[K],fu[K]; void pre(){ zh[0]=1;fu[0]=1; for(reg i=1;i<=k+2;++i) { inc(Y[i],ad(Y[i-1],qm(i,k))); zh[i]=mul(zh[i-1],qm(i)); fu[i]=mul(fu[i-1],qm(mod-i)); } // prt(Y,1,k+2); // prt(zh,1,k+2); // prt(fu,1,k+2); } int pr[K],bc[K]; int num; int calc(ll n){ // cout<<" calc "<<n<<endl; if(n<=sqr) { // cout<<" sum "<<sum[n]<<endl; return sum[n]; } if(n<=k+2) { // cout<<" Y "<<Y[n]<<endl; return Y[n]; } n%=mod; // ++num; pr[0]=1;bc[k+3]=1; for(reg i=1;i<=k+2;++i){ pr[i]=mul(pr[i-1],n-i); } for(reg i=k+2;i>=1;--i){ bc[i]=mul(bc[i+1],n-i); } int ret=0; for(reg i=1;i<=k+2;++i){ inc(ret,mul(Y[i],pr[i-1],bc[i+1],zh[i-1],fu[k+2-i])); } // cout<<" ret "<<ret<<endl; return ret; } int ans; int tot=0; void dfs(int d,ll x,int h){ if(d>cnt) return; for(reg y=d;y<=cnt&&(ll)x<=(n/pri[y]/pri[y]);++y){ int tmp=mul(h,sub(mi[y],mul(mi[y],mi[y]))); ll now=x*pri[y]*pri[y]; for(reg e=2;;++e){ // cout<<" now "<<now<<endl; inc(ans,mul(tmp,calc(n/now))); dfs(y+1,now,tmp); if(now>n/pri[y]) break; now*=pri[y]; } } } int main(){ rd(n);rd(k); sqr=sqrt(n); sieve(sqr); pre(); // cout<<" nd "<<endl; ans=calc(n); dfs(1,1,1); cout<<ans;//<<" time "<<tot<<" cha "<<num<<endl; return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* */
对于LOJ简单的函数——Min_25筛
也可以用powerful number做的
令$G(x)=\phi$,p=2需要特殊处理。
这里h就没有O(1)公式了,还不能前缀和优化,只能暴力卷积
h[p][e]表示h(p^e)那么实际上计算h的总复杂度也不是很大。
好吧其实感觉挺大的,但是看不懂zzq在写什么神仙操作。。。。
提供了一种处理前缀和的思路。
也是构造,但是复杂度降低的原因是,g前缀和好求,h有值的位置不多。且$h(p)=0$
一般地,如果可以找到合适的g去拟合(大概至少使得g(p)=f(p)?且g为积性函数),使得能构造出积性函数h使得$h(p)=0$,那么都是可以这样做的!
如果没有好的办法,$h(p^e)$理论上可以暴力计算。