[杜教筛] BZOJ3512. DZY Loves Math IV

来源:互联网 发布:ipa版本淘宝的微淘 编辑:程序博客网 时间:2024/05/16 11:20

%%%Vectorxj
http://blog.csdn.net/vectorxj/article/details/78857079

#include <cstdio>#include <iostream>#include <algorithm>#include <map>using namespace std;typedef pair<int,int> par;const int N=1e6+10,P=1e9+7;int mu[N],p[N],phi[N],fac[N],pre[N];inline void Pre(const int &n){  phi[1]=1; mu[1]=1; fac[1]=1;  for(int i=2;i<=n;i++){    if(!p[i])      p[++*p]=i,phi[i]=i-1,mu[i]=-1,fac[i]=i;    for(int j=1;j<=*p && 1LL*i*p[j]<=n;j++){      p[i*p[j]]=1;      if(i%p[j]){    fac[i*p[j]]=fac[i]*p[j];    phi[i*p[j]]=phi[i]*phi[p[j]];    mu[i*p[j]]=-mu[i];      }      else{    fac[i*p[j]]=fac[i];    phi[i*p[j]]=phi[i]*p[j];    break;      }    }  }  for(int i=1;i<=n;i++) pre[i]=(phi[i]+pre[i-1])%P;}map<par,int> M;map<int,int> S;int Phi(int n){  if(n<=1000000) return pre[n];  if(S.count(n)) return S[n];  int ret=1LL*n*(n+1)/2%P;  for(int i=2,j;i<=n;i=j+1){    j=n/(n/i);    ret=(ret-1LL*(j-i+1)*Phi(n/i))%P;  }  return S[n]=ret;}int Sum(int n,int m){  if(!n || !m) return 0;  if(M.count(par(n,m))) return M[par(n,m)];  if(n==1) return Phi(m);  if(m==1) return phi[n];  if(!mu[n])    return 1LL*n/fac[n]*Sum(fac[n],m)%P;  int ret=0;  for(int i=1;1LL*i*i<=n;i++){    if(n%i) continue;    ret=(ret+1LL*phi[n/i]*Sum(i,m/i))%P;    if(i*i!=n)      ret=(ret+1LL*phi[i]*Sum(n/i,m/(n/i)))%P;  }  return M[par(n,m)]=ret;}int main(){  freopen("1.in","r",stdin);  freopen("1.out","w",stdout);  int n,m,ans=0; scanf("%d%d",&n,&m); Pre(1000000);  for(int i=1;i<=n;i++)    ans=(ans+Sum(i,m))%P;  printf("%d\n",(ans+P)%P);  return 0;}
原创粉丝点击