51nod1227 平均最小公倍数

$
Ans(l,r)=ans(r)-ans(l-1)
\\
ans(n)=\sum\limits_{i=1}^n
\sum\limits_{j=1}^i
\frac{j}{gcd(i,j)}
\\=
\sum\limits_{g=1}^n
\sum\limits_{i=1}^{n/g}
\sum\limits_{j=1}^n
j\cdot [gcd(i,j)=1]
\\=
\sum\limits_{g=1}^n
\sum\limits_{i=1}^{n/g}
\sum\limits_{j=1}^n
j\sum\limits_{d|i\wedge d|j}\mu(d)
\\=
\sum\limits_{g=1}^n
\sum\limits_{d=1}^{n/g}
\mu(d)d\sum\limits_{i=1}^{n/gd}
\sum\limits_{j=1}^n
j
\\=
\sum\limits_{t=1}^n
\sum\limits_{d|t}\mu(d)d
\sum\limits_{i=1}^{n/t}
\sum\limits_{j=1}^n
j
\\=
\sum\limits_{t=1}^n
a(t)S(n/t)
,\;
a(n)=\sum\limits_{d|n}\mu(d)d
,\;
S(n)=n(n+1)(n+2)/6=n(n+1)(2n+1)/12+n(n+1)/4
\\
A(n)=\sum\limits_{i=1}^na(i)
\\
A(n)=n-\sum\limits_{i=2}^nA(n/i)i
$

时间复杂度$O(n^{2/3})$

#include<bits/stdc++.h>
typedef unsigned long long i64;
const int P=1e9+7,I2=(P+1)/2,I6=(P+1)/6;
int B;
const int M=3e6+7;
int ps[M/10],pp=0;
bool np[M];
int va[M],vA[5007],n0;
int S3(int n){
    return i64(n)*(n+1)%P*(n+2)%P;
}
int A(int n){
    if(n<=B)return va[n];
    int&w=vA[n0/n];
    if(w)return w;
    i64 s=0;
    for(int l=1,r,c;l<n;l=r){
        r=n/(c=n/(l+1));
        if((s+=i64(r+l+1)*(r-l)%P*A(c))>i64(1.5e19))s%=P;
    }
    s=s%P*I2%P;
    return w=(n-s+P)%P;
}
int F(int n){
    i64 s=0;
    for(int l=0,r,c,s0=0,s1;l<n;l=r){
        r=n/(c=n/(l+1));
        s1=A(r);
        if((s+=i64(s1-s0+P)*S3(c))>i64(1.5e19))s%=P;
        s0=s1;
    }
    return s%P*I6%P;
}
void pre(){
    va[1]=1;
    for(int i=2;i<=B;++i){
        if(!np[i]){
            ps[pp++]=i;
            va[i]=P+1-i;
        }
        for(int j=0,k;j<pp&&(k=i*ps[j])<=B;++j){
            np[k]=1;
            if(i%ps[j]){
                va[k]=va[i]*i64(P+1-ps[j])%P;
            }else{
                va[k]=va[i];
                break;
            }
        }
    }
    for(int i=1;i<=B;++i)if((va[i]+=va[i-1])>=P)va[i]-=P;
}
int cal(int n){
    n0=n;
    memset(vA,0,sizeof(vA));
    return F(n);
}
int main(){
    int l,r;
    scanf("%d%d",&l,&r);
    B=pow(r,2./3.)*2.5;
    pre();
    printf("%d\n",(cal(r)-cal(l-1)+P)%P);
    return 0;
}