51nod 1238 最小公倍数之和 V3 题解

题意简述

$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n} \operatorname{lcm}(i,j)$

$n\le 10^{10}$ （简短题意，恶臭规模）

思路

总的来说，就是先把式子变一变，然后巧妙的用欧拉函数 $\varphi$ 加上杜教筛和整除分块来算。

变形

首先我们知道 $\operatorname{lcm}(a,b)=\operatorname{lcm}(b,a)$

所以，当 $i\not=j$ 时候，我们计算 $j<i$ 的情况，然后 $\times 2$ 即可

对于 $i=j$ 的情况特判，这时候产生的贡献和为 $\dfrac{n(n+1)}{2}$

所以，答案化为

$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i-1} \operatorname{lcm}(i,j) + \dfrac{n(n+1)}{2}$

我们都知道 $\operatorname{lcm}(i,j)=\dfrac{ij}{\gcd(i,j)}$。这个再换一下即可。

开始硬♂上

开局一公式，结论全靠推。让我们开心的推倒式子吧~

首先枚举 $gcd$，然后枚举互质的倍数 $i,j$。

原式

$=\sum\limits_{d=1}^{n}\sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{i} [\gcd(i,j)=1] ijd$
$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i} [\gcd(i,j)=1] ij \sum\limits_{d=1}^{n/i} d$ （枚举一对互质的 $i,j$ ，计算有多少个 $d$ 算到了这对 $i,j$，加起来）
设 $S_k(n)=\sum\limits_{i=1}^{n} i^k$ 。显然， $k=1,2,3$ 时这个式子都能 $O(1)$ 算。
接下来变成
$=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{i}[\gcd(i,j)=1] ij S_1(n/i)$ （恭喜你少了一个 $d$）
考虑如何求 $\sum\limits_{j=1}^{i} [\gcd(i,j)=1]j$，也就是求 $i$ 以内和 $i$ 互质的数的和。

如何求 i 以内和 i 互质的数的和（会的跳过）

首先数量一共有 $\varphi(i)$ 个。
然后我们打表发现，这东西和等差数列有一个很类似的性质：我们把 $[1,i]$ 中和 $i$ 互质的 $j$ 从小到大排一排，换一行，再从大到小排一排，对应位置和相等（都是 $i$）！
根据高斯当年想出来的办法，和就是（每个对应位置的和）乘以（一共有多少项），然后除二。
换到这题来，这个和就是 $\dfrac{1}{2}i\times \varphi(i)$

它等于 $\dfrac{1}{2}i\times \varphi(i)$。带入原式：

$=\sum\limits_{i=1}^{n} i\times S_1(n/i) \times \dfrac{i\times \varphi(i)}{2}$ （恭喜你又少了一个 $j$，它现在是一维了）
$=\dfrac{1}{2}\sum\limits_{i=1}^{n} i^2\varphi(i)\times S_1(n/i)$

$S_1(n/i)$ 显然可以整除分块，但是整除分块我们要考虑一个东西： 如何求 $i^2\varphi(i)$ 的前缀和？显然，用杜教筛求一下和即可。

具体怎么杜教筛 （老样子，会的跳过）

设 $f(n)=n^2\varphi(n)$
关键就是要配一个好算的 $g$ ，使得 $f\times g$ 也是一个好算的函数 $h$。

由狄利克雷卷积的定义，
$h(n)=\sum\limits_{d|n} f(d)g(\dfrac{n}{d})$
$=\sum\limits_{d|n} d^2\varphi(i)g(\dfrac{n}{d})$
我们发现，令 $g(n)=n^2$，就能消掉一个 $d^2$，原式继续变成：
$=\sum\limits_{d|n} n^2\varphi(i)$ （注意 $g(\dfrac{n}{d})$ 头顶上还有一个 $n$，乘出来变成 $n^2$）
=$n^2\sum\limits_{d|n} varphi(d)=n^3$

所以说， $h(n)=f\times (g(n)=n^2)=(n^3)$

然后用杜教筛的传统艺能就行了。

总复杂度非常恐怖，但是的确能通过 $n=10^10$ 的数据。一定注意处处取膜

代码

复制

#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
    #define N   20000007ll
    #define mod 1000000007ll
    #define i6  166666668ll
    #define i2  500000004ll
    #define int long long
    #define F(i,l,r) for(int i=l;i<=r;++i)
    #define D(i,r,l) for(int i=r;i>=l;--i)
    #define Fs(i,l,r,c) for(int i=l;i<=r;c)
    #define Ds(i,r,l,c) for(int i=r;i>=l;c)
    #define MEM(x,a) memset(x,a,sizeof(x))
    #define FK(x) MEM(x,0)
    #define Tra(i,u) for(int i=G.Start(u),v=G.To(i);~i;i=G.Next(i),v=G.To(i))
    #define p_b push_back
    #define sz(a) ((int)a.size())
    #define iter(a,p) (a.begin()+p)
    int I()
    {
        int x=0;char c=getchar();int f=1;
        while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar();
        while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
        return (x=(f==1)?x:-x);
    }
    void Rd(int cnt,...)
    {
        va_list args; va_start(args,cnt);
        F(i,1,cnt) {int* x=va_arg(args,int*);(*x)=I();}
        va_end(args);
    }
    bitset<N> notp; int primes[N/5];
    int phi[N];
    void Init()
    {
        int n=2e7;
        notp[1]=1; phi[1]=1;
        int &cnt=primes[0];
        F(i,2,n)
        {
            if (!notp[i]) {primes[++cnt]=i; phi[i]=i-1;}
            for(int j=1;j<=cnt and i*primes[j]<=n;++j)
            {
                int u=primes[j];
                notp[i*u]=1;
                if (i%u==0){phi[i*u]=phi[i]*u; break;}
                else phi[i*u]=phi[i]*phi[u];
            }
        }
        F(i,1,n) phi[i]=(phi[i-1]+1ll*phi[i]*i%mod*i%mod)%mod;
    }

    int n;
    void Input()
    {
        n=I();
    }

    map<int,int> rec; 
    int S1(int n) {n%=mod; return 1ll*n%mod*(n+1)%mod*i2%mod;}
    int S2(int n) {n%=mod; return 1ll*n%mod*(n+1)%mod*(2ll*n+1)%mod*i6%mod;}
    int S3(int n) {n%=mod; return 1ll*S1(n)*S1(n)%mod;} // 我就是这里没取膜，调了两小时
    int f(int n) //杜教筛求 f(n)=n^2*phi(n) 的前缀和
    {
        if (n<=2e7) return phi[n];
        if (rec[n]) return rec[n];
        int ans=0;
        for(int l=2,r;l<=n;l=r+1)
        {
            r=n/(n/l);
            ans=(ans+f(n/l)*(S2(r)-S2(l-1)))%mod;
        } 
        n%=mod;
        return rec[n]=(S1(n)*S1(n)-ans)%mod;
    }
    void Soviet()
    {
        int ans=0;
        for(int l=2,r;l<=n;l=r+1)
        {
            r=n/(n/l);
            ans=(ans+S1(n/l)*(f(r)-f(l-1))%mod*i2%mod)%mod;
        }// 整除分块+杜教筛
        printf("%lld\n",((2ll*ans+S1(n))%mod+mod)%mod);
    }

    #define Flan void
    Flan IsMyWife()
    {
        Init();
        Input();
        Soviet();
    }
    #undef int //long long
}
int main()
{
    Flandre_Scarlet::IsMyWife();
    getchar();getchar();
    return 0;
}