[51nod1238]最小公倍数之和V3

来自FallDream的博客，未经允许，请勿转载，谢谢。

------------------------------------------------------------------------------------------

题意:求$$\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j)   \\\    n\leqslant 10^{10}$$

题解：题目即求$$\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{i*j}{gcd(i,j)}$$

$$=\sum_{d=1}^{n}d*\sum_{i=1}^{\lfloor n/d\rfloor}\sum_{j=1}^{\lfloor n/d\rfloor}i*j*[gcd(i,j)=1]$$

已知$$\sum_{i=1}^{n}i*[gcd(i,n)=1]=\frac{n*\varphi(n)}{2}$$ 

所以所求即为$$\sum_{d=1}^{n}d*\sum_{i=1}^{\lfloor n/d\rfloor}i*i*\varphi(i)$$

 

$\lfloor\frac{n}{d}\rfloor$只有$\sqrt(n)$种取值，那么我们考虑快速求出$g(i)=i^{2}*\varphi(i)$的前缀和$S(i)$。

 

$$\sum_{n|d}\varphi(d)=n$$ 

$$\sum_{i=1}^{n}\sum_{d|n}\varphi(d)=\frac{n(n+1)}{2}$$

$$\sum_{i=1}^{n}\sum_{d|n}\varphi(d)*i^{2}=\frac{n^{2}*(n+1)^{2}}{4}$$

$$\sum_{i=1}^{n}\sum_{d|n}g(d)*(\frac{i}{d})^{2}=\frac{n^{2}*(n+1)^{2}}{4}$$

$$\sum_{i=1}^{n}\sum_{d=1}^{\lfloor n/i\rfloor}g(d)*i^{2}=\frac{n^{2}*(n+1)^{2}}{4}$$

$$S(i)=\frac{n^{2}*(n+1)^{2}}{4}-\sum_{i=2}^{n}i^{2}*S(\lfloor n/i\rfloor)$$

这个可以在$O\left(n^{\frac{2}{3}}\right)$时间内做完。此题得解。

-------------------------------------------------------------------

我好菜啊，推了好久.....

-----

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<map>
#define MAXN 5000000
#define ll long long
#define mod 1000000007
#define ditoly 6666666
#define inv2 500000004
#define inv4 250000002
#define inv6 166666668
using namespace std;
inline ll read()
{
    ll x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0'; ch=getchar();}
    return x*f;
}

struct mymap{
    ll x,ans;int next;
}e[8000000];
ll phi[MAXN+5],n,ans=0;
int s[MAXN/5],num=0,head[ditoly+5];
bool b[MAXN+5];

inline ll getcube(ll x){x%=mod;return x*(x+1)%mod*(x<<1|1)%mod*inv6%mod;}

inline void ins(ll x,ll sum)
{
    int j=x%ditoly;
    e[++num]=(mymap){x,sum,head[j]};head[j]=num;
}

inline ll getsq(ll x){x%=mod;x=x*(x+1)%mod;return x*x%mod*inv4%mod;}

ll calc(ll x)
{
    if(x<=MAXN)return phi[x];
    for(int i=head[x%ditoly];i;i=e[i].next)
        if(e[i].x==x)return e[i].ans;
    ll last,sum=getsq(x);
    for(ll i=2;i<=x;i=last+1)
    {
        last=x/(x/i);
        sum-=(getcube(last)-getcube(i-1)+mod)%mod*calc(x/i)%mod;
        while(sum<0)sum+=mod;
    }
    ins(x,sum);
    return sum;
}

int main()
{
    n=read();phi[1]=1;
    for(int i=2;i<=MAXN;i++)
    {
        if(!b[i]) phi[i]=i-1,s[++num]=i;
        for(int j=1;s[j]*i<=MAXN;j++)
        {
            b[s[j]*i]=1;
            if(i%s[j]==0){phi[s[j]*i]=phi[i]*s[j];break;}
            phi[s[j]*i]=phi[i]*(s[j]-1);
        }
    }num=0;
    for(int i=2;i<=MAXN;i++)
        phi[i]=(phi[i-1]+1LL*i*i%mod*phi[i]%mod)%mod;
    for(ll i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);ll x=(n/i)%mod;
        ans+=x*(x+1)%mod*inv2%mod*((calc(last)-calc(i-1)+mod)%mod)%mod;
        while(ans>=mod)ans-=mod;
    }
    cout<<(ans+mod)%mod;
    return 0;
}