就是实现一下Min-25筛 筛积性函数的操作
首先要得到
$G(M,j)=sum_{t=j}^{cnt} sum_{e=1}^{p_t^{e+1}<=M} [phi(p_t^e)*G([M/(p_t^e)],t+1)+phi(p_t^{(e+1)})]$
$+(F(M)-(F(p_{j-1})))$
先要预处理后面的部分,得到$F(M)$和$F(p_{j-1})$
$F(p_{j-1})$可以直接筛素数的时候前缀和计算一下
$F(M)$就要利用第一步的筛法了
发现,除了2之外的质数都是奇数,所以f(p^1)=p xor 1=p-1
对于2要特判
对于G,直接根据式子大力计算即可。
递归处理。由于值还是比较分散的,所以没有记忆化的必要。(而且状态很多,对空间极为不友好)
剪枝:pri[t]的平方大于n就不用继续算了。
代码:
#include<bits/stdc++.h> #define il inline #define reg register int #define int long long #define numb (ch^'0') using namespace std; typedef long long ll; il void rd(int &x){ char ch;bool fl=false; while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } namespace Miracle{ const int N=5e5+2; const int M=5e5+2; const int mod=1e9+7; int pri[M],tot; int sum[M];//pre of prime bool vis[N]; int sqr; ll f[N],g[N],h[N]; void sieve(int n){ for(reg i=2;i<=n;++i){ if(!vis[i]){ vis[i]=1; pri[++tot]=i; } for(reg j=1;j<=tot;++j){ if(i*pri[j]>n) break; vis[i*pri[j]]=1; if(i%pri[j]==0) break; } } for(reg i=1;i<=tot;++i){ sum[i]=(sum[i-1]+pri[i])%mod; g[i]=(g[i-1]+(pri[i]^1))%mod; } } int id1[N],id2[N]; ll val[N]; ll n; int S(int x,int j){ if(x<=1||x<pri[j]) return 0; cout<<" xx "<<x<<" jj "<<j<<endl; int d=(x<=sqr)?id1[x]:id2[n/x]; int ret=(f[d]-g[j-1]+mod)%mod; for(reg t=j;t<=tot&&pri[t]*pri[t]<=x;++t){ int now=pri[t]; for(reg e=1;now*pri[t]<=x;now=now*pri[t],++e){ ret=(ret+(pri[t]^e)*S(x/now,t+1)%mod+(pri[t]^(e+1))%mod)%mod; } } return ret; } int main(){ scanf("%lld",&n); if(n==1){ puts("1");return 0; } sqr=sqrt(n); // cout<<" sqr "<<sqr<<endl; sieve(sqr); // cout<<" after sieve "<<endl; int m=0; for(ll i=1,x;i<=n;i=x+1){ x=n/(n/i); val[++m]=n/i; if(val[m]<=sqr) id1[val[m]]=m; else id2[n/val[m]]=m; } for(reg i=1;i<=m;++i){ f[i]=val[i]-1;h[i]=(((ll)val[i]%mod*(val[i]%mod+1))/2-1+mod)%mod; } for(reg j=1;j<=tot;++j){ for(reg i=1;i<=m&&(ll)pri[j]*pri[j]<=val[i];++i){ int to=(val[i]/pri[j])<=sqr?id1[val[i]/pri[j]]:id2[n/(val[i]/pri[j])]; f[i]=(f[i]-(f[to]-(j-1))+mod+mod)%mod; h[i]=(h[i]-pri[j]*(h[to]-sum[j-1]+mod)%mod+mod)%mod; } } for(reg i=1;i<=m;++i){ if(val[i]>=2) f[i]=(h[i]-f[i]+2+mod)%mod; else f[i]=0; } //cout<<" after prewrk "<<endl; printf("%lld",(S(n,1)+1)%mod); return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2019/1/13 17:03:03 */