51nod 1227 平均最小公倍数

51nod 1227 平均最小公倍数

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

$$A(n)=\frac{1}{n}\sum_{i=1}^nlcm(i,n)$$

$$\sum_{k=1}^nA(k)\\=\sum_{k=1}^n\frac{1}{k}\sum_{i=1}^klcm(i,k)\\=\sum_{k=1}^n\frac{1}{k}\sum_{i=1}^k\frac{ki}{gcd(i,k)}\\=\sum_{k=1}^n\sum_{i=1}^k\frac{i}{gcd(i,k)}\\=\sum_{k=1}^n\sum_{d|k}\sum_{i=1}^k\frac{i}{d}[gcd(i,k)=d]\\=\sum_{k=1}^n\sum_{d|k}\sum_{i=1}^{\frac{k}{d}}i[i\perp \frac{k}{d}]$$

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

因为 $a\perp n$ 则：$(n-a)\perp n$

所以： $$\sum_{i=1}^n[i\perp n]i=\frac{\varphi(n)*n}{2}$$

特别的，$n=1$时： $$\sum_{i=1}^n[i\perp n]i=\frac{\varphi(n)*n+1}{2}$$

所以： $$\sum_{k=1}^n\sum_{d|k}\sum_{i=1}^{\frac{k}{d}}i[i\perp \frac{k}{d}]\\=\sum_{k=1}^n\Big(\frac{1}{2}+\sum_{d|k}\frac{\varphi(d)d}{2}\Big)\\=\frac{n}{2}+\frac{1}{2}\sum_{d=1}^n\varphi(d)d\Big\lfloor\frac{n}{d}\Big\rfloor$$

令： $$g(n)=\varphi(n)n\\id_k(n)=n^k$$

显而下面的狄利克雷积成立： $$g*id_1=id_2$$

$f$的前缀和为：$S_f$

则： $$S_g(n)=\frac{n(n+1)(2n+1)}{6}-\sum_{i=2}^niS_g(\Big\lfloor\frac{n}{i}\Big\rfloor)$$

那么对于： $$\sum_{d=1}^n\varphi(d)d\Big\lfloor\frac{n}{d}\Big\rfloor\\=\sum_{L=1}^{m-1}L\Big(S_g(\Big\lfloor\frac{n}{L}\Big\rfloor)-S_g(\Big\lfloor\frac{n}{L+1}\Big\rfloor)\Big)+\sum_{d=1}^{\lfloor\frac{n}{m}\rfloor}g(d)\Big\lfloor\frac{n}{d}\Big\rfloor$$

代码：

#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;
}