2693: jzptab/2154: Crash的数字表格 莫比乌斯反演

神题……不会做啊……抄题解吧……

以下题解来自网上（某爬虫网站的我也不知道原作者是谁）

---

下边均设$n<=m$

$$\sum_{i=1}^{n} \sum_{j=1}^{m} lcm(i, j) = \sum_{i=1}^{n} \sum_{j=1}^{m} \frac{ij}{(i,j)}$$

然后我们要想枚举$d=(i,j)$，那么就要确定$ij$怎么取，显然我们只需要先除去$i$和$j$的$d$，也就是$(i/d,j/d)=1$就行了，那么设

$$F(x, y) = \sum_{i=1}^{x} \sum_{j=1}^{y} ij[(i,j)=1]$$

那么原式变成

$$\sum_{d=1}^{n} \frac{d^2 F(\lfloor \frac{n}{d} \rfloor, \lfloor \frac{m}{d} \rfloor)}{d} = \sum_{d=1}^{n} d F(\lfloor \frac{n}{d} \rfloor, \lfloor \frac{m}{d} \rfloor)$$

考虑求$F(x,y)$

$$\begin{align} F(x, y) & = \sum_{i=1}^{x} \sum_{j=1}^{y} ij[(i,j)=1] \\ & = \sum_{i=1}^{x} \sum_{j=1}^{y} ij \sum_{d|(i,j)} \mu (d) \\ & = \sum_{d=1}^{x} \mu (d) \sum_{d|i}^{x} i \sum_{d|j}^{y} j \\ & = \sum_{d=1}^{x} \mu (d) d^2 \sum_{i=1}^{\lfloor \frac{x}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{y}{d} \rfloor} 1 \\ & = \sum_{d=1}^{x} \mu (d) d^2 \frac{\lfloor \frac{x}{d} \rfloor (\lfloor \frac{x}{d} \rfloor +1)}{2} \frac{\lfloor \frac{y}{d} \rfloor (\lfloor \frac{y}{d} \rfloor +1)}{2} \\ \end{align}$$

带回原式得

$$\begin{align} & \sum_{d=1}^{n} d F(\lfloor \frac{n}{d} \rfloor, \lfloor \frac{m}{d} \rfloor) \\ = & \sum_{d=1}^{n} d \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \mu (i) i^2 \frac{\lfloor \frac{\lfloor \frac{n}{d} \rfloor}{i} \rfloor (\lfloor \frac{\lfloor \frac{n}{d} \rfloor}{i} \rfloor +1)}{2} \frac{\lfloor \frac{\lfloor \frac{m}{d} \rfloor}{i} \rfloor (\lfloor \frac{\lfloor \frac{m}{d} \rfloor}{i} \rfloor +1)}{2} \\ = & \sum_{d=1}^{n} d \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \mu (i) i^2 \frac{\lfloor \frac{n}{di} \rfloor (\lfloor \frac{n}{di} \rfloor +1)}{2} \frac{\lfloor \frac{m}{di} \rfloor (\lfloor \frac{m}{di} \rfloor +1)}{2} \\ \end{align}$$

现在已经可以$O(\sqrt n \sqrt n) = O(n)$单次查询了，但是不够理想，我们继续化简

令$T=di$，则$i|T, d=T/i$，换掉指标，得

$$\begin{align} & \sum_{d=1}^{n} d \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \mu (i) i^2 \frac{\lfloor \frac{n}{di} \rfloor (\lfloor \frac{n}{di} \rfloor +1)}{2} \frac{\lfloor \frac{m}{di} \rfloor (\lfloor \frac{m}{di} \rfloor +1)}{2} \\ = & \sum_{T=1}^{n} \frac{\lfloor \frac{n}{T} \rfloor (\lfloor \frac{n}{T} \rfloor +1)}{2} \frac{\lfloor \frac{m}{T} \rfloor (\lfloor \frac{m}{T} \rfloor +1)}{2} \sum_{i|T} \frac{T}{i} \mu (i) i^2 \\ = & \sum_{T=1}^{n} \frac{\lfloor \frac{n}{T} \rfloor (\lfloor \frac{n}{T} \rfloor +1)}{2} \frac{\lfloor \frac{m}{T} \rfloor (\lfloor \frac{m}{T} \rfloor +1)}{2} T \sum_{i|T} \mu (i) i \\ \end{align}$$

设$g(T)=T \sum_{i|T} \mu (i) i$，我们再设$f(T)=\sum_{i|T} \mu (i) i$那么$g(T)=Tf(T)$，考虑求$f(T)$

在线性筛中，外层为$k$，内层为$p_y$，所以求$f(kp_y)=\sum_{i|T} \mu(i) i$
当$p_y|k$时
当$i$取的数的因子中不包含新加入的$p_y$时，答案就是$f(k)$
当$i$取包含新加入的因子$p_y$时，由于此时$p_y$指数已经$>=2$，所以$\mu (i)=0$，因此贡献为0
综上，当$p_y|k$时，答案为$f(k)$

当$p_y \nmid k$时
当$i$取的数的因子中不包含新加入的$p_y$时，同上，答案是$f(k)$
当$i$取的数的因子包含新加入的$p_y$时，由于指数为$1$，所以我们考虑$i=ap_y$，原式变为

$$\begin{align} & \sum_{i|T} \mu(i) i \\ = & \sum_{ap_y|kp_y} \mu (ap_y) ap_y \\ = & p_y \sum_{a|k} \mu (a) \mu(p_y) a \\ = & -p_y \sum_{a|k} \mu(a) a \\ = & -p_y f(k) \\ \end{align}$$

综上，当$p_y \nmid k$时，答案为$(1-p_y)f(k)$

然后线性筛随便搞搞即可，最后答案就是

$$\sum_{T=1}^{n} \frac{\lfloor \frac{n}{T} \rfloor (\lfloor \frac{n}{T} \rfloor +1)}{2} \frac{\lfloor \frac{m}{T} \rfloor (\lfloor \frac{m}{T} \rfloor +1)}{2} g(T)$$

---

$$F(x,y)= \sum_{i=1}^{x} \sum_{j=1}^{y} ij[(i,j)=1] = \sum_{i=1}^{x} \sum_{j=1}^{y} ij \sum_{d|(i,j)} \mu (d) \\$$
感觉这一步转化好神啊，完全想不到。

#include<iostream>
#include<cstdio>
#define ll long long 
#define M 100000009
#define N 10000005
using namespace std;
int prime[N],mobius[N];
bool flag[N];
ll f[N];
inline ll read()
{
    ll a=0,f=1; char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1; c=getchar();}
    while (c>='0'&&c<='9') {a=a*10+c-'0'; c=getchar();}
    return a*f;
}
inline void prepare()
{
    mobius[1]=1; f[1]=1;
    for (int i=2;i<N;i++)
    {
        if (!flag[i]) prime[++prime[0]]=i,mobius[i]=-1,f[i]=1-i;
        for (int j=1;j<=prime[0]&&i*prime[j]<N;j++)
        {
            flag[i*prime[j]]=1;
            if (i%prime[j]==0) {mobius[i*prime[j]]=0; f[i*prime[j]]=f[i]; break;}
            mobius[i*prime[j]]=-mobius[i]; f[i*prime[j]]=(f[i]*f[prime[j]])%M;
        }
    }
    for (int i=1;i<N;i++)
        f[i]=(f[i-1]+f[i]*i%M)%M;
}   
inline ll sum(ll n,ll m)
{
    return (n*(n+1)/2%M)*(m*(m+1)/2%M)%M;
}
inline ll work(ll n,ll m)
{
    ll ans=0;
    for (ll i=1,pos;i<=n;i=pos+1)
    {
        pos=min(n/(n/i),m/(m/i));
        ans=(ans+(f[pos]-f[i-1]+M)%M*sum(n/i,m/i)%M)%M;
    }
    return ans;
}
int main()
{
    int testcase=read();
    prepare();
    while (testcase--)
    {
        ll n=read(),m=read();
        if (n>m) swap(n,m);
        printf("%lld\n",(work(n,m)+M)%M);
    }
    return 0;
}