【51nod 1238】最小公倍数之和 V3（杜教筛）-优快云博客

本文链接：https://blog.youkuaiyun.com/Coldfresh/article/details/81273715

本文解析了一道关于计算小于等于N的所有数两两之间最小公倍数之和的数论算法题，通过数学推导简化了问题，并给出了具体的实现代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

基准时间限制：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，牵扯上了莫比乌斯函数的，但是这个就不能这样搞了，反正先推一推吧。

a n s = \sum i = 1 n \sum j = 1 n i j g c d ( i , j ) = \sum d = 1 n \sum i = 1 n \sum j = 1 n [g c d (i, j) = = d] i j d = \sum d = 1 n d \sum i = 1 ⌊ n d ⌋ \sum j = 1 ⌊ n d ⌋ [g c d (i, j) = = 1] i j

$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 [g c d (i, j) = = 1] i j

$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 [g c d (i, j) = = 1] i j

$F'(x)=\sum_{i=1}^n\sum_{j=1}^i[gcd(i,j)==1]ij$
易知

F(x)=2F′(x)−1F(x)=2F′(x)−1 $F(x)=2F'(x)-1$ ，所以我们先来考虑

F′(x)F′(x) $F'(x)$
那很显然啊…

F' (x) = 1 2 + \sum i = 1 n i 2 ϕ ( i ) 2

$F'(x)=\frac{1}{2}+\sum_{i=1}^n\frac{i^2\phi(i)}{2}$

所以：

F (x) = \sum i = 1 n i 2 ϕ (i)

$F(x)=\sum_{i=1}^ni^2\phi(i)$ .
感觉这个式子好简洁，那么就是求：

a n s = \sum i = 1 n i F (⌊ n i ⌋)

$ans=\sum_{i=1}^niF(\lfloor\frac{n}{i}\rfloor)$
似乎没什么实质化简，用杜教搞一搞啦。
想到

∑d|nϕ(d)=n∑d|nϕ(d)=n $\sum_{d|n}\phi(d)=n$ ,看看能不能利用这个性质了。
设：

f(x)=x2ϕ(x),g(x)=x2f(x)=x2ϕ(x),g(x)=x2 $f(x)=x^2\phi(x),g(x)=x^2$ ,我们拿这两个函数做狄利克雷卷积，有：

f * g (n) = \sum d | n d 2 ϕ (d) n 2 d 2 = n 3 = \sum d | n g (d) f (n d)

$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 n i 3 = \sum i = 1 n \sum j | i g (j) f (i j) = \sum j = 1 n g (j) \sum k = 1 ⌊ n j ⌋ f (k) = \sum j = 1 n g (j) F (⌊ n j ⌋) = n 2 ( n + 1 ) 2 4

$\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) = n 2 ( n + 1 ) 2 4 - \sum j = 2 n g (j) F (⌊ n j ⌋)

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