莫比乌斯反演，杜教筛

Posted 2020-09-30 z1j1n1

BZOJ4176

#include <cstdio>
#include <map>
#define LL long long 
using namespace std;

  map <LL,LL> mpa;
  map <LL,LL> mpb;

  const LL mo=1e9+7;

  int b[10000001],phi[10000001],ss[10000001];
  int miu[10000001],summiu[10000001],cnt,n;

  void euler(int n=1e7){
    b[1]=1;miu[1]=1;summiu[1]=1;
    for (LL i=2;i<=n;i++){
      if (b[i]==0) {ss[++cnt]=i;miu[i]=-1;}
      summiu[i]=summiu[i-1]+miu[i];
         for (LL j=1;j<=cnt;j++){
            if (i*ss[j]>n) break;
            b[i*ss[j]]=1;
            if (i%ss[j]==0) {miu[i*ss[j]]=0;break;}
            else miu[i*ss[j]]=miu[i]*(-1);
         }
       }
  } 
  
  LL getsummiu(LL num){
      if (num<=1e7) return(summiu[num]);
      if (mpb[num]!=0) return(mpa[num]);//可以不使用map，可能出现的值为n/i的根号n个值
      LL res=1;//筛phi时为n*(n-1)/2
      for (LL i=2,last;i<=num;i=last+1){
        last=min((num/(num/i)),num);
      res-=(last-i+1)*getsummiu(num/i)%mo;
      res%=mo;
    }
    mpb[num]=1;mpa[num]=res;
    return(res);
  }
  
  LL f(LL n){
      LL ret=0;
      for (int i=1,last;i<=n;i=last+1){
        last=min(n/(n/i),n);
        ret+=(last-i+1)*(n/i)%mo;
        ret%=mo;
    }
    return(ret*ret%mo);
  }
  
  int main(){
      scanf("%d",&n);
      euler();
      
      LL ans=0;
      for (int d=1,last;d<=n;d=last+1){
        last=min(n/(n/d),n);
      ans+=(getsummiu(last)-getsummiu(d-1))*f(n/d)%mo;
      ans%=mo;
    }
    ans%=mo;ans+=mo;ans%=mo;
    printf("%lld\n",ans);
  }