洛谷P3327 [SDOI2015]约数个数和 题解

题目链接：P3327 [SDOI2015]约数个数和

题意：

设 $d(x)$ 为 $x$ 的约数个数，给定 $n,m$，求

$$\sum_{i=1}^n\sum_{j=1}^md(ij)$$

输入格式：

输入文件包含多组测试数据。

第一行，一个整数 $T$，表示测试数据的组数。

接下来的 $T$ 行，每行两个整数 $n,m$。

输出格式：

$T$ 行，每行一个整数，表示你所求的答案。

数据范围：

$1\le T,n,m \le 50000$。

感觉我前几篇题解的无解释推导有点不可读

所以这篇还是写的清楚一点叭

这题需要一个引理（感觉没必要记住证明的方法）

引理：设 $d(x)$ 为 $x$ 的约数个数，则有

$$d(ij)=\sum_{x \mid i}\sum_{y \mid j}[\gcd(x,y)=1]$$

证明：摘自 https://siyuan.blog.luogu.org/solution-p3327

我们考虑把每个因子一一映射。

如果 $ij$ 的因子 $k$ 中有一个因子 $p^c$，$i$ 中有因子 $p^a$，$j$ 中有因子 $p^b$。我们规定：

• 如果 $c\le a$，那么在 $i$ 中选择。
• 如果 $c>a$，那么我们把 $c$ 减去 $a$，在 $j$ 中选择 $p^{c-a}$（在 jj 中选择 $p^e$ 表示的是 $p^{a+e}$ ）

对于 $ij$ 的因子 $k$ 的其他因子同理。于是对于任何一个 $k$ 有一个唯一的映射，且每一个选择对应着唯一的 $k$。

通过如上过程，我们发现：对于 $ij$的因子 $k=\prod {p_i}^{c_i}$，我们不可能同时在 $i$ 和 $j$ 中选择 $p_i$（优先在 $i$ 中选择，如果不够就只在 $j$ 中选择不够的指数），故 $x$ 和 $y$ 必须互质。

等式得证。

有了这个引理，就可以开始搞了

不妨假设 $n \le m$

$$\sum_{i=1}^{n}\sum_{j=1}^{m}d(ij)$$

代入

$$\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y \mid j}[\gcd(x,y)=1]$$

莫比乌斯反演

$$\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y \mid j}\sum_{d\mid \gcd(x,y)}\mu(d)$$

这个 $d \mid \gcd(x,y)$ 看着很难受，把它换掉

$$\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y \mid j}\sum_{d=1}^{n}\mu(d)[d\mid \gcd(x,y)]$$

啊它怎么又回来了，不急，继续推

移一下

$$\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y \mid j}[d\mid \gcd(x,y)]$$

可以发现这个正整数 $x$ 的取值范围在 $[1,n]$

而它的出现次数其实就是 $[1,n]$ 中有多少个数有因数 $x$

也就是 $\sum_{x=1}^{n}[x \mid n] = \sum_{x=1}^{n}\left\lfloor\frac{n}{x}\right\rfloor$

$y$ 同理，则可以得到

$$\sum_{d=1}^{n}\mu(d)\sum_{x=1}^{n}\sum_{y=1}^{m}\left\lfloor{\dfrac{n}{x}}\right\rfloor\left\lfloor{\dfrac{m}{y}}\right\rfloor[d\mid \gcd(x,y)]$$

可以发现，满足 $d\mid \gcd(x,y)$ 的 $d$ 都是 $x,y$ 的因数，反过来也可以说 $x,y$ 都是 $d$ 的倍数。

于是我们改为枚举 $d$ 的倍数 $dx$ 和 $dy$ 即可（注意这里 $x,y$ 的意义改变了）

$$\sum_{d=1}^{n}\mu(d)\sum_{x=1}^{\left\lfloor{\frac{n}{d}}\right\rfloor}\left\lfloor{\dfrac{n}{dx}}\right\rfloor\sum_{y=1}^{\left\lfloor{\frac{m}{d}}\right\rfloor}\left\lfloor{\dfrac{m}{dy}}\right\rfloor$$

然后就可以数论分块啦！

注意到这里的 $dx,dy$ 其实你把它的 $d$ 放上去就是可以算的 $\left\lfloor\frac{\left\lfloor\frac{n}{d}\right\rfloor}{x}\right\rfloor$ 了

可以预处理一下这个东西为 $g\left(\left\lfloor\frac{n}{d}\right\rfloor\right)$ ，然后每组询问查表就好了。

时间复杂度 $\mathcal{O}(T\sqrt{n} + n\sqrt{n})$

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
typedef long long ll;
#define INF 0x3f3f3f3f3f3f3f3f
namespace FastIO
{
    #define gc() readchar()
    #define pc(a) putchar(a)
    #define SIZ (int)(1e6+15)
    char buf1[SIZ],*p1,*p2;
    char readchar()
    {
        if(p1==p2)p1=buf1,p2=buf1+fread(buf1,1,SIZ,stdin);
        return p1==p2?EOF:*p1++;
    }
    template<typename T>void read(T &k)
    {
        char ch=gc();T x=0,f=1;
        while(!isdigit(ch)){if(ch=='-')f=-1;ch=gc();}
        while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=gc();}
        k=x*f;
    }
    template<typename T>void write(T k)
    {
        if(k<0){k=-k;pc('-');}
        static T stk[66];T top=0;
        do{stk[top++]=k%10,k/=10;}while(k);
        while(top){pc(stk[--top]+'0');}
    }
}using namespace FastIO;
#define N (int)(5e4+15)
int prime[N],sum[N],pcnt,mu[N],g[N];
bool ck[N];
int solve(int n)
{
    int res=0;
    for(int l=1,r; l<=n; l=r+1)
    {
        r=(n/(n/l));
        res+=(r-l+1)*(n/l);
    }
    return res;
}
void Mobius()
{
    mu[1]=1;
    for(int i=2; i<=N-5; i++)
    {
        if(!ck[i])
        {
            prime[++pcnt]=i;
            mu[i]=-1;
        }
        for(int j=1; j<=pcnt&&i*prime[j]<=N-5; j++)
        {
            int pos=i*prime[j];
            ck[pos]=1;
            if(i%prime[j])
            {
                mu[pos]=-mu[i];
            }else
            {
                mu[pos]=0;
                break;
            }
        }
    }
    for(int i=1; i<=N-5; i++)
        sum[i]+=sum[i-1]+mu[i];
    for(int i=1; i<=N-5; i++)
        g[i]=solve(i);
}
int Q,n,m;
signed main()
{
    Mobius();
    read(Q);
    while(Q--)
    {
        read(n);read(m);
        if(n>m)swap(n,m);
        int res=0;
        for(int l=1,r; l<=n; l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            res+=(sum[r]-sum[l-1])*g[n/l]*g[m/l];
        }
        write(res);pc('\n');
    }
    return 0;
}