[杜教筛] 51Nod 1238 最小公倍数之和 V3

最新推荐文章于 2019-02-20 20:08:09 发布

里阿奴摩西

最新推荐文章于 2019-02-20 20:08:09 发布

阅读量1.2k

点赞数

分类专栏：莫比乌斯反演&杜教筛

本文链接：https://blog.csdn.net/u014609452/article/details/54747619

版权

莫比乌斯反演&杜教筛专栏收录该内容

38 篇文章 0 订阅

订阅专栏

比最大公约数之和要难搞一些

吉丽博客传送门：http://jiruyi910387714.is-programmer.com/posts/195270.html

这道题两个要点

首先

∑1<=i<=n [(n,i)==1]*i == 1/2*([n==1]+n*phi(n))

还有就是对于phi(i)*i*i也就是phi·id·id的前缀和

我们单独卷一下phi

(phi*1)·id·id=id·id·id

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<tr1/unordered_map>
using namespace std;
using namespace std::tr1;
typedef long long ll;

const int P=1000000007;
const int inv2=(P+1)/2;
const int inv3=(P+1)/3;

const int maxn=5e6;
int prime[(int)1e6],num;
int vst[maxn+5],phi[maxn+5],sum[maxn+5];

inline void Pre(){
  phi[1]=1;
  for (int i=2;i<=maxn;i++){
    if (!vst[i]) phi[i]=i-1,prime[++num]=i;
    for (int j=1;j<=num && (ll)i*prime[j]<=maxn;j++){
      vst[i*prime[j]]=1;
      if (i%prime[j]==0){
	phi[i*prime[j]]=phi[i]*prime[j];
	break;
      }else{
	phi[i*prime[j]]=phi[i]*(prime[j]-1);
      }
    }
  }
  for (int i=1;i<=maxn;i++) sum[i]=((ll)phi[i]*i%P*i%P+sum[i-1])%P;
}

inline ll sum1(ll n){ return (n+1)%P*(n%P)%P*inv2%P;}
inline ll sum2(ll n){ return (n%P)*((n+1)%P)%P*((2*n+1)%P)%P*inv2%P*inv3%P; }
inline ll sum3(ll n){ return sum1(n)*sum1(n)%P; }
inline ll sum1(ll l,ll r){ return (sum1(r)+P-sum1(l-1))%P; }
inline ll sum2(ll l,ll r){ return (sum2(r)+P-sum2(l-1))%P; }
inline ll sum3(ll l,ll r){ return (sum3(r)+P-sum3(l-1))%P; }

unordered_map<ll,int> S;

inline int Sum(ll n){
  if (n<=maxn) return sum[n];
  if (S.find(n)!=S.end()) return S[n];
  int tem=sum3(n);
  ll l,r;
  for (l=2;l*l<=n;l++) (tem+=P-l%P*(l%P)%P*Sum(n/l)%P)%=P;
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,(tem+=P-sum2(l,r)*Sum(t)%P)%=P;
  return S[n]=tem;
}

int main(){
  ll n; int Ans=0,tem=0;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  Pre();
  scanf("%lld",&n);
  Ans=sum1(n)*inv2%P;
  ll l,r;
  tem=0;
  //for (l=1;l<=n;l++) {
  //  (tem+=(ll)l*l%P*phi[l]%P*sum1(n/l))%=P;
  //printf("%d %d\n",tem,l);
  //}
  for (l=1;l*l<=n;l++) (tem+=(ll)l*l%P*phi[l]%P*sum1(n/l)%P)%=P;
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,(tem+=(Sum(r)+P-Sum(l-1))%P*sum1(t)%P)%=P;
  (Ans+=(ll)tem*inv2%P)%=P;
  Ans=(Ans*2%P+P-sum1(n))%P;
  printf("%d\n",Ans);
  return 0;
}