1. 程式人生 > >[杜教篩 約數和字首和] 51Nod 1220 約數之和

[杜教篩 約數和字首和] 51Nod 1220 約數之和

套用公式後反演 然後杜教篩求和


比較有意思的是其間我算了兩個值

後來發現這兩個值竟然是相等的 都可以由約數和字首和推導過來


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

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

int prime[2000000],num;
int miu[maxn+5];
ll sumd[maxn+5],t1[maxn+5],t2[maxn+5]; int d[maxn+5];

inline void Pre(){
  miu[1]=1; sumd[1]=1;
  for (int i=2;i<=maxn;i++){
    if (!d[i]) d[i]=prime[++num]=i,sumd[i]=1+i,t1[i]=1+i,t2[i]=i,miu[i]=-1;
    for (int j=1;j<=num && (ll)i*prime[j]<=maxn;j++){
      d[i*prime[j]]=prime[j]; int k=i*prime[j];
      if (i%prime[j]==0){
	miu[k]=0;
	t2[k]=t2[i]*prime[j],t1[k]=t1[i]+t2[k],sumd[k]=sumd[i]/t1[i]*t1[k];
	break;
      }
      miu[i*prime[j]]=miu[i]*miu[prime[j]];
      sumd[k]=sumd[i]*sumd[prime[j]],t1[k]=1+prime[j],t2[k]=prime[j];
    }
  }
  for (int i=1;i<=maxn;i++)
    miu[i]=((ll)i*((P+miu[i])%P)%P+miu[i-1])%P,(sumd[i]+=sumd[i-1])%=P;
}

unordered_map<ll,int> S;

inline void add(int &x,int y){
  x+=y; if (x>=P) x-=P;
}

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

inline int Sum(ll n){
  if (n<=maxn) return miu[n];
  if (S.find(n)!=S.end()) return S[n];
  int tem=1; ll l,r;
  for (l=2;l*l<=n;l++) add(tem,P-l*Sum(n/l)%P);
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,add(tem,P-(ll)sum1(l,r)*Sum(t)%P);
  return S[n]=tem;
}
/*
inline int F(ll n){
  int t1=0,t2=0; ll l,r;
  for (l=1;l*l<=n;l++) add(t1,l%P*(n/l)%P),add(t2,sum1(n/l)%P);
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,add(t1,(ll)sum1(l,r)*(n/l)%P),add(t2,(r-l+1)%P*sum1(n/l)%P);
  return (ll)t1*t2%P;
}
*/

inline int F(ll n){
  if (n<=maxn) return (ll)sumd[n]*sumd[n]%P;
  int t1=0; ll l,r;
  for (l=1;l*l<=n;l++) add(t1,l%P*(n/l)%P);
  for (ll t=n/l;l<=n;l=r+1,t--)
    r=n/t,add(t1,(ll)sum1(l,r)*(n/l)%P);
  return (ll)t1*t1%P;
}

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