51nod 1227 平均最小公倍数

原题链接：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1227

A(n)=1n∑i=1nlcm(i,n)

∑k=1nA(k)=∑k=1n1k∑i=1klcm(i,k)=∑k=1n1k∑i=1kkigcd(i,k)=∑k=1n∑i=1kigcd(i,k)=∑k=1n∑d|k∑i=1kid[gcd(i,k)=d]=∑k=1n∑d|k∑i=1kdi[i⊥kd]

对于计算[1,n]与n互素的数字之和可以类比等差数列的倒叙相加。

因为 a⊥n 则：(n−a)⊥n

所以：

∑i=1n[i⊥n]i=φ(n)∗n2

特别的，n=1时：

∑i=1n[i⊥n]i=φ(n)∗n+12

所以：

∑k=1n∑d|k∑i=1kdi[i⊥kd]=∑k=1n(12+∑d|kφ(d)d2)=n2+12∑d=1nφ(d)d⌊nd⌋

令：

g(n)=φ(n)nidk(n)=nk

显而下面的狄利克雷积成立：

g∗id1=id2

f的前缀和为：Sf

则：

Sg(n)=n(n+1)(2n+1)6−∑i=2niSg(⌊ni⌋)

那么对于：

∑d=1nφ(d)d⌊nd⌋=∑L=1m−1L(Sg(⌊nL⌋)−Sg(⌊nL+1⌋))+∑d=1⌊nm⌋g(d)⌊nd⌋

代码：

#include <stdio.h>#include <algorithm>#include <string.h>#include <cmath>#define MAXN 1000000#define SQR 40000#define S(a) ((LL)(a)*((a)+1)%P*Iv[2]%P)using namespace std;typedef long long LL;const int P=1e9+7;int tmp[SQR];int g[MAXN];int phi[MAXN];int Iv[SQR];void init (){    Iv[1]=1;    for(int i=2;i<SQR;i++)Iv[i]=(P-(LL)Iv[P%i]*(P/i)%P)%P;    for(int i=1;i<MAXN;i++)    {        phi[i]=i-phi[i];        for(int j=i<<1;j<MAXN;j+=i)phi[j]+=phi[i];        phi[i]=(LL)phi[i]*i%P;        g[i]=phi[i]+g[i-1];        if(g[i]>=P)g[i]-=P;    }}void clat(int n,int d){    int m=sqrt(n)+1,ans=0;    for(int L=1;L<m;L++)    {        ans+=(LL)g[L]*((S(n/L)-S(n/(L+1))+P)%P)%P;        if(ans>=P)ans-=P;    }    for(int i=n/m;i>1;i--)    {        int u=n/i;        if(u<MAXN)            ans+=(LL)g[u]*i%P;        else            ans+=(LL)tmp[i*d]*i%P;        if(ans>=P)ans-=P;    }    m=(LL)n*(n+1)%P*(2*n+1)%P*Iv[6]%P;    tmp[d]=(P+m-ans)%P;}void sum(int n){    if(n<MAXN)return ;    for(int i=n/MAXN;i;i--)clat(n/i,i);    return ;}int solve(int n){    sum(n);    int m=sqrt(n)+1,ans=0;    for(int L=1;L<m;L++)    {        int u1=n/L;        int u2=n/(L+1);        if(u1<MAXN)            u1=g[u1];        else            u1=tmp[L];        if(u2<MAXN)            u2=g[u2];        else            u2=tmp[L+1];        ans+=(LL)L*(u1-u2+P)%P;        if(ans>=P)ans-=P;    }    for(int i=n/m;i;i--)    {        ans+=(LL)phi[i]*(n/i)%P;        if(ans>=P)ans-=P;    }    return ans;}int main (){    init();    int a,b;    scanf("%d %d",&a,&b);    a--;    int A=solve(a);    int B=solve(b);    A=(LL)A*Iv[2]%P;    A+=(LL)a*Iv[2]%P;    if(A>=P)A-=P;    B=(LL)B*Iv[2]%P;    B+=(LL)b*Iv[2]%P;    if(B>=P)B-=P;    printf("%d\n",(B-A+P)%P);    return 0;}