bzoj 2693 jzptab / 洛谷 P1829 Crash的数字表格 题解【莫比乌斯反演】【狄利克雷卷积】

点击量：173

莫比乌斯反演进阶+优化。

目录

Description

求
$$
\sum_{i=1}^n\sum_{j=1}^m\operatorname{lcm}(i,j)
$$
对 $ 100000009$ 取模输出。

Input

一个正整数 $ T$ 表示数据组数，

接下来 $ T$ 行，每行两个正整数，表示 $ N,M​$。

Output

$ T$ 行 每行一个整数 表示第i组数据的结果

Sample Input

1
4 5

Sample Output

122

HINT

$ T\le 10000$

$ N, M\le 10000000​$

Source

版权所有者：倪泽堃

题解：

$ \operatorname{lcm}$ 和 $ \gcd$ 关系很近，因此在反演中有可能有较大联系。要尽可能把 $ \operatorname{lcm}$ 转化为 $ \gcd$。

首先考虑题目，求
$$
\sum_{i=1}^n\sum_{j=1}^m\operatorname{lcm}(i,j)
$$
令 $ a=\min(n,m),b=\max(n,m)$（方便后面的定界），并根据定义，有
$$
\begin{aligned}
\text{原式}&=\sum_{i=1}^a\sum_{j=1}^b\operatorname{lcm}(i,j)\\
&=\sum_{i=1}^a\sum_{j=1}^b\frac{ij}{\gcd(i,j)}
\end{aligned}
$$
此时出现了$ \gcd​$，根据之前做过题的套路，可以把它转化为布尔表达式，原式变为
$$
\sum_{k=1}^{a}\sum_{i=1}^a\sum_{j=1}^b[\gcd(i,j)=k]\frac{ij}k
$$
把 $ k$ 提前，得
$$
\sum_{k=1}^a\frac 1k\sum_{i=1}^a\sum_{j=1}^b[\gcd(i,j)=k]ij
$$
再把 $ k$ 提出来，<这里不提出来也能推下去，但是到最后式子会变麻烦？>
$$
\begin{aligned}
&\sum_{k=1}^a\frac 1k\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\sum_{j=1}^\left\lfloor\frac bk\right\rfloor[\gcd(i,j)=1]ijk^2\\
=&\sum_{k=1}^ak\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\sum_{j=1}^\left\lfloor\frac bk\right\rfloor[\gcd(i,j)=1]ij
\end{aligned}
$$
然后是套路？令
$$
f(n)=\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\sum_{j=1}^\left\lfloor\frac bk\right\rfloor[\gcd(i,j)=n]ij
$$
则 $ f(1)$ 就是上面式子的后半段。继续反演设（$ k$ 这个参数可以视作临时全局变量）
$$
\begin{aligned}
g(n)&=\sum_{n|d}f(d)\\
&=\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\sum_{j=1}^\left\lfloor\frac bk\right\rfloor[\gcd(i,j)=n,2n,3n,\cdots]ij\\
&=\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\sum_{j=1}^\left\lfloor\frac bk\right\rfloor[n|\gcd(i,j)]ij
\end{aligned}
$$
把 $ n$ 提上来
$$
g(n)=\sum_{i=1}^\left\lfloor\frac a{kn}\right\rfloor\sum_{j=1}^\left\lfloor\frac b{kn}\right\rfloor[1|\gcd(i,j)]ij
$$
布尔表达式恒成立，则 $ g(n)$ 为
$$
(1+2+\cdots+\left\lfloor\frac a{kn}\right\rfloor)\times(1+2+\cdots+\left\lfloor\frac b{kn}\right\rfloor)=\frac{\left\lfloor\frac a{kn}\right\rfloor(\left\lfloor\frac a{kn}\right\rfloor+1)\left\lfloor\frac b{kn}\right\rfloor(\left\lfloor\frac b{kn}\right\rfloor+1)}{4}
$$
反演回去，得
$$
\begin{aligned}
f(n)&=\sum_{n|d}\mu\left(\frac dn\right)g(d)\\
&=\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\mu(i)g(ni)
\end{aligned}
$$
因为我们要求 $ f(1)$ 来得到原式，所以代入 $ n=1$
$$
\begin{aligned}
f(1)&=\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\mu(i)g(i)\\
&=\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\mu(i)i^2\frac{\left\lfloor\frac a{ki}\right\rfloor\left(\left\lfloor\frac a{ki}\right\rfloor+1\right)\left\lfloor\frac b{ki}\right\rfloor\left(\left\lfloor\frac b{ki}\right\rfloor+1\right)}{4}
\end{aligned}
$$
带回原式得
$$
\sum_{k=1}^ak\sum_{i=1}^\left\lfloor\frac ak\right\rfloor\mu(i)i^2\frac{\left\lfloor\frac a{ki}\right\rfloor\left(\left\lfloor\frac a{ki}\right\rfloor+1\right)\left\lfloor\frac b{ki}\right\rfloor\left(\left\lfloor\frac b{ki}\right\rfloor+1\right)}{4}
$$

此时预处理 $ \mu(i)i^2$ 的前缀和，上面的式子就能用两个数论分块很快地求出来。

首先 $ \left\lfloor\frac ak \right\rfloor$ 可以数论分块，然后把 $ \left\lfloor\frac ak \right\rfloor$ 看成一个整体，即把 $ \left\lfloor\frac a{ki} \right\rfloor$ 看成 $ \left\lfloor\frac {\left\lfloor\frac ak\right\rfloor} i\right\rfloor$ ，再数论分块。

此时时间复杂度为 $ O(\sqrt n\sqrt n)=O(n)$ ，可以过掉 Crash 的数字表格 这道题了。

但是jzptab这道题有 $ T\le 10000$ 个询问，因此需要在 $ O(\sqrt n)$ 的时间内回答每个询问。

我们继续对原式进行变形。令
$$
S(n)=\frac{n(n+1)}2
$$
则原式为
$$
\sum_{k=1}^a\sum_{i=1}^\left\lfloor\frac ak\right\rfloor k\mu(i)i^2S\left(\left\lfloor\frac a{ki}\right\rfloor\right)S\left(\left\lfloor\frac b{ki}\right\rfloor\right)
$$
考虑令 $ ki=T$，则前面的枚举就是在枚举 $ T$ 和它的因数 $ i$。可以转化为
$$
\sum_{T=1}^a\sum_{d|T}Td\mu(d)S\left(\left\lfloor\frac aT\right\rfloor\right)S\left(\left\lfloor\frac bT\right\rfloor\right)
$$
此时后面的 $ S$ 中只包含 $ T$，那么我们把这个式子提到外面来，有
$$
\sum_{T=1}^aS\left(\left\lfloor\frac aT\right\rfloor\right)S\left(\left\lfloor\frac bT\right\rfloor\right)\sum_{d|T}Td\mu(d)
$$
此时外面的 $ \sum$ 可以数论分块了，那么后面的式子能不能 $ O(1)$ 完成呢？

看到后面的 $ \sum$ 是枚举因数形式，联想到卷积，但是卷积一般是 $ g\left(\frac Td\right)h(d)$ 的形式，因此把它转换过去。

令（新的 $ f$）
$$
\begin{aligned}
f(T)&=\sum_{d|T}Td\mu(d)\\
&=\sum_{d|T}\frac Td d^2\mu(d)
\end{aligned}
$$
可以看出 $ g(T)=T,h(T)=T^2\mu(T)$，则 $ f=g*h$。而 $ g(T)$ 和 $ h(T)$ 都是积性函数，因此 $ f(T)$ 也是积性函数，就可以类似莫比乌斯函数一样筛出来了。

通过 $ f$ 在质数幂处的值可以求出其所有函数值，$ f\left(p^k\right)=p^k-p^{k+1}$。

此时原式为
$$
\sum_{T=1}^aS\left(\left\lfloor\frac aT\right\rfloor\right)S\left(\left\lfloor\frac bT\right\rfloor\right)f(T)
$$

其中 $ S$ 和 $ f$ 都是 $ O(1)$ 的。

Code：

#include<cstdio>
#include<cstring>
#define p 100000009
int Min(int x,int y){return x<y?x:y;}
int Plus(int x,int y){return x+y>=p?x+y-p:x+y;}
int Mul(int x,int y){return 1ll*x*y%p;}
int pri[10001000],f[10001000],cnt=0;
bool is[10001000];
int main()
{
    is[0]=is[1]=1;
    f[1]=1;
    for(int i=2;i<=10000000;++i)
    {
        if(!is[i])
        {
            pri[++cnt]=i;
            f[i]=Plus(i,p-Mul(i,i));
        }
        for(int j=1;j<=cnt&&i*pri[j]<=10000000;++j)
        {
            is[i*pri[j]]=1;
            if(i%pri[j])
                f[i*pri[j]]=Mul(f[i],f[pri[j]]);
            else
            {
                f[i*pri[j]]=Mul(f[i],pri[j]);
                break;
            }
        }
        //f[i]=Plus(f[i-1],f[i]);注意不能顺便加，因为后面可能还要用
    }
    for(int i=1;i<=10000000;++i)
        f[i]=Plus(f[i-1],f[i]);
    int T,n,m;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d",&n,&m);
        if(n>m)
        {
            int t=n;
            n=m;
            m=t;
        }
        int ans=0;
        for(int i=1;i<=n;++i)
        {
            int N=n/i,M=m/i;
            int nxt=Min(n/N,m/M);
            ans=Plus(ans,Mul(Mul(Mul(N,N+1),Mul(M,M+1)),Plus(f[nxt],p-f[i-1])));
            i=nxt;
        }
        printf("%d\n",Mul(ans,75000007));
    }
    return 0;
}