Description：

求： $$ \sum_{n=L}^R\sum_{x=1}^n\sum_{y=1}^n[\text{lcm}(x,y)=n] $$ $1\leqslant L\leqslant R\leqslant 10^{11}$

Solution：

先化式子： $$ \begin{align} &ans=\sum_{n=L}^Rf(n),f(n)=\sum_{x=1}^n\sum_{y=1}^n[\text{lcm}(x,y)=n]\\ &设g(n)=\sum_{x=1}^n\sum_{y=1}^n[\text{lcm}(x,y)|n]=\sum_{x=1}^n\sum_{y=1}^n[x|n][y|n]=\sigma_0^2(n)\\ &那么有：g(n)=\sum_{d|n}f(d)\\ &根据莫比乌斯反演定理：f(n)=\sum_{d|n}\mu(\frac n d)g(d)=\sum_{d|n}\mu(\frac n d)\sigma_0^2(d)\\ \end{align} $$ 然后我们要证明一个结论： $$ \sum_{d|n}\mu(\frac nd)\sigma_0^2(d)=\sigma_0(n^2) $$ 假设上式成立，先莫比乌斯反演： $$ \begin{align} \sigma_0^2(n)&=\sum_{d|n}\sigma_0(d^2)\\ \prod_{i=1}^{k_n}(q_{n,i}+1)(q_{n,i}+1)&=\sum_{d|n}\prod_{i=1}^{k_n}(2q_{d,i}+1)\\ \prod_{i=1}^{k_n}(q_{n,i}+1)(q_{n,i}+1)&=\prod_{i=1}^{k_n}\sum_{j=0}^{q_{n,i}}(2j+1)\\ \prod_{i=1}^{k_n}(q_{n,i}^2+2q_{n,i}+1)&=\prod_{i=1}^{k_n}\frac{(2q_{n,i}+1+1)\times(q_{n,i}+1)}{2}\\ \end{align} $$ 所以我们要求的就是： $$ \sum_{i=1}^n\sigma_0(n^2) $$ 考虑怎么把平方去掉： $$ \begin{align} \sigma_0(n^2)&=\prod_{i=1}^k(2q_i+1)\\ &=\sum_{S\subseteq \{1,2,\dots,k\}}2^{|S|}\prod_{i\in S}q_i\\ &=\sum_{d|n}2^{w(d)} \end{align} $$ $w(n)$ 代表 $n$ 的不同质因子数， $2^{w(n)}$ 也就是 $n$ 的无平方因子数，那么有： $$ 2^{w(n)}=\sum_{d|n}\mu^2(d) $$ 设 $f(n)=\sigma_0(n^2),g(n)=2^{w(n)},h(n)=\mu^2(n)$ 。

那么有： $$ f(n)=1\times g(n),g(n)=1\times h(n),f(n)=h(n)\times 1\times1=h(n)\times \sigma_0(n) $$ 根据这个式子我们可以线性筛，具体来说： $$ \begin{align} f(xy)&=f(x)f(y)(\gcd(x,y)=1)\\ f(p)&=h(p)+\sigma_0(p)=1+2=3\\ f(p^c)&=\sum_{i=0}^ch(p^i)\sigma_0(p^{c-i})=h(1)\sigma_0(p^c)+h(p)\sigma_0(p^{c-1})=c+1+c=2c+1 \end{align} $$ 但是还不够： $$ \sum_{i=1}^nf(i)=\sum_{i=1}^n\sum_{d|i}h(d)\sigma_0(\frac id)=\sum_{d=1}^n\mu^2(d)\sum_{k=1}^{\lfloor\frac n d\rfloor}\sigma_0(k)=\sum_{d=1}^n\mu^2(d)s(\lfloor\frac n d\rfloor) $$ 所以如果我们能计算出 $\sigma_0$ 和 $\mu^2$ 的前缀和就可以除法分块了。 $$ \sum_{i=1}^n\sigma_0(i)=\sum_{i=1}^n\lfloor\frac n i\rfloor $$ 于是 $\sigma_0$ 可以除法分块。

由于 $\sum \mu^2(d)$ 表示的是 $n$ 以内无平方因子数，因此我们可以考虑容斥： $$ \sum_{i=1}^{\sqrt n}\mu(i)\lfloor\frac n{i^2}\rfloor $$ 也就是枚举 $n$ 的平方因子，可以发现容斥系数是莫比乌斯函数，这样所有平方因子都会被计算 $0$ 次，而 $1$ 会被计算一次。

但是这样显然复杂度还是不对，于是我们可以先线性筛预处理 $n^{\frac 2 3}$ 的部分，然后后面的出发分块做就行了。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<map>

#include<cctype>

#include<cstring>

using namespace std;

#define MAXN 5000010

#define N 5000000

typedef long long ll;

int mindivc[MAXN];

int prime[MAXN],tot = 0;

bool isprime[MAXN];

int mu[MAXN],mus[MAXN];

ll sig[MAXN];

#define R register

#define I inline

I void init()

{

    for(R int i = 2;i <= N;++i)isprime[i] = true;

    sig[1] = mu[1] = mus[1] = 1;

    for(R int i = 2;i <= N;++i)

    {

        if(isprime[i])prime[++tot] = i,mindivc[i] = 1,mu[i] = -1,sig[i] = 2;

        for(R int j = 1;j <= tot && i * prime[j] <= N;++j)

        {

            R int k = i * prime[j];

            isprime[k] = false;

            if(i % prime[j] == 0)

            {

                mindivc[k] = mindivc[i] + 1;mu[k] = 0;

                sig[k] = sig[i] / (mindivc[i] + 1) * (mindivc[i] + 2);

                break;

            }

            else

            {

                mindivc[k] = 1;mu[k] = -mu[i];sig[k] = sig[i] * 2;

            }

        }

        mus[i] = mus[i - 1] + abs(mu[i]);

        sig[i] = sig[i - 1] + sig[i];

    }

    return;

}

I ll sigma(ll n)

{

    if(n <= N)return sig[n];

    R ll ans = 0;

    for(R ll l = 1,r;l <= n;l = r + 1)

    {

        r = n / (n / l);

        ans = ans + (r - l + 1) * (n / l);

    }

    return ans;

}

struct hash

{

    #define MOD 1000007

    int head[MOD];

    ll state[MOD],val[MOD];

    int nxt[MOD],tot;

    hash(){memset(head,0,sizeof(head));tot = 0;}

    I ll find(ll x)

    {

        R int mod = x % MOD;

        for(R int i = head[mod];i != 0;i = nxt[i])if(state[i] == x)return val[i];

        return -1;

    }

    I void set(ll p,ll v)

    {

        R int mod = p % MOD;

        ++tot;state[tot] = p;val[tot] = v;nxt[tot] = head[mod];head[mod] = tot;

        return;

    }

}smu;

I ll summu(ll n)

{

    if(n <= N)return mus[n];

    R ll val = smu.find(n);

    if(val != -1)return val;

    R ll ans = 0;

    for(R ll i = 1;i * i <= n;++i)

    {

        ans += mu[i] * (n / (i * i));

    }

    smu.set(n,ans);

    return ans;

}

I ll calc(ll n)

{

    R ll ans = 0;

    for(R ll l = 1,r;l <= n;l = r + 1)

    {

        r = n / (n / l);//if(r / 1000 != (l - 1) / 1000)cout << l << " " << r << endl;

        ans = ans + (summu(r) - summu(l - 1)) * sigma(n / l);

    }

    return ans;

}

ll query(ll n)

{

    if(n == 0)return 0;

    return (calc(n) + n) / 2;

}

int main()

{

    init();

    ll l,r;scanf("%lld%lld",&l,&r);

    cout << query(r) - query(l - 1) << endl;

    return 0;

}