【51nod 1238】 最小公倍数之和 V3（杜教筛）

基准时间限制：8 秒 空间限制：262144 KB 分值: 640 难度：8级算法题
出一个数N，输出小于等于N的所有数，两两之间的最小公倍数之和。

相当于计算这段程序（程序中的lcm(i,j)表示i与j的最小公倍数）：
由于结果很大，输出Mod 1000000007的结果。

G=0;
for(i=1;i<=N;i++)
for(j=1;j<=N;j++)
{
G = (G + lcm(i,j)) % 1000000007;
}
Input
输入一个数N。(2 <= N <= 10^10)
Output
输出G Mod 1000000007的结果。
Input示例
4
Output示例
72

之前推过一个差不多的，但是那个的数量级是10的5次的，而且不是正方形，是给出的n和m，牵扯上了莫比乌斯函数的，但是这个就不能这样搞了，反正先推一推吧。

$$ans=\sum_{i=1}^n\sum_{j=1}^n\frac{ij}{gcd(i,j)}\\ =\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)==d]\frac{ij}{d}\\ =\sum_{d=1}^{n}d\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{n}{d} \right \rfloor}[gcd(i,j)==1]ij$$
设： $$F(n)=\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)==1]ij$$
观察发现若i==j，且gcd为1，只有i和j均为1时，所以我们不妨考虑这个式子，设： $$F'(x)=\sum_{i=1}^n\sum_{j=1}^i[gcd(i,j)==1]ij$$
易知$F(x)=2F'(x)-1$，所以我们先来考虑$F'(x)$
那很显然啊… $$F'(x)=\frac{1}{2}+\sum_{i=1}^n\frac{i^2\phi(i)}{2}$$

所以：

$$F(x)=\sum_{i=1}^ni^2\phi(i)$$ .
感觉这个式子好简洁，那么就是求： $$ans=\sum_{i=1}^niF(\lfloor\frac{n}{i}\rfloor)$$
似乎没什么实质化简，用杜教搞一搞啦。
想到$\sum_{d|n}\phi(d)=n$,看看能不能利用这个性质了。
设：$f(x)=x^2\phi(x),g(x)=x^2$,我们拿这两个函数做狄利克雷卷积，有：
$$f*g(n)=\sum_{d|n}d^2\phi(d)\frac{n^2}{d^2}=n^3=\sum_{d|n}g(d)f(\frac{n}{d})$$
对卷积和做一下前缀和：
$$\sum_{i=1}^ni^3=\sum_{i=1}^n\sum_{j|i}g(j)f(\frac{i}{j})=\sum_{j=1}^ng(j)\sum_{k=1}^{\lfloor\frac{n}{j}\rfloor}f(k)=\sum_{j=1}^ng(j)F(\lfloor\frac{n}{j}\rfloor)=\frac{n^2(n+1)^2}{4}$$
所以： $$g(1)F(n)=\frac{n^2(n+1)^2}{4}-\sum_{j=2}^ng(j)F(\lfloor\frac{n}{j}\rfloor)$$

推导部分就到这了，现在就是递归分块求了。

代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define maxx 5000000
#define mod 1000000007
#define ll long long
using namespace std;
bool isP[maxx];
int prime[400000];
int cnt;
int phi[maxx];
ll F[maxx];
ll inv2=500000004,inv4=250000002,inv6=166666668;
ll get(ll x)
{
    x%=mod;
    return (x+1)*x%mod*(2*x+1)%mod*inv6%mod;
}
void init()
{
    phi[1]=1;
    for(int i=2;i<maxx;i++)
    {
        if(!isP[i]){prime[cnt++]=i;phi[i]=i-1;}
        for(int j=0;j<cnt&&(ll)i*prime[j]<maxx;j++)
        {
            isP[i*prime[j]]=true;
            if(i%prime[j])
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            else
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
        }
    }
    for(ll i=1;i<maxx;i++)
        F[i]=i*i%mod*phi[i]%mod;
    for(int i=1;i<maxx;i++)
        F[i]=(F[i]+F[i-1])%mod;
}
map<ll,ll>M;
ll work(ll x)
{
    if(x<maxx)return F[x];
    if(M[x])return M[x];
    ll  t=x%mod;
    ll ans=t*t%mod*(t+1)%mod*(t+1)%mod*inv4%mod;
    for(ll i=2,last;i<=x;i=last+1)
    {
        last=x/(x/i);
        ans=(ans-(get(last)-get(i-1))*work(x/i)%mod)%mod;
    }
    if(ans<0)ans=(ans+mod)%mod;
    M[x]=ans;
    return ans;
}
ll _get(ll a,ll b)
{

    return (b-a+1)%mod*((a+b)%mod)%mod*inv2%mod;
}
int main()
{
    init();
    //cout<<cnt<<endl;
    ll n;
    cin>>n;
    ll ans=0;
    for(ll i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans=(ans+_get(i,last)*work(n/i)%mod)%mod;
    }
    cout<<ans<<endl;
    return 0;
}