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

作者: wjyyy 分类: 数学,狄利克雷卷积,莫比乌斯反演,解题报告 发布时间: 2019-02-19 08:46

点击量:52

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

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

说点什么

avatar
  Subscribe  
提醒
/* */