洛谷 [SDOI2015]约数个数和 解题报告

[SDOI2015]约数个数和

题目描述

设\(d(x)\)为\(x\)的约数个数，给定\(N,M\)，求$ \sum\limits^{N_{i=1}\sum\limits}M_{j=1}d(ij)$

输入输出格式

输入格式：

输入文件包含多组测试数据。第一行，一个整数\(T\)，表示测试数据的组数。接下来的\(T\)行，每行两个整数\(N,M\)。

输出格式：

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

说明

\(1 \le N, M \le 50000\)

\(1 \le T \le 50000\)

---

Solution

引理\(1\)：

\[\sum_{d|gcd(a,b)} \mu(d)=[gcd(a,b)=1] \]

引理\(2\):

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

可以通过\(d\)唯一分解后的计算式感性理解一下

剩下的暴力推个式子

\[\sum_{i=1}^a\sum_{j=1}^bd(ij) \]

\[=\sum_{i=1}^a\sum_{j=1}^b\sum_{x|i}\sum_{y|j}\sum_{d|gcd(x,y)}\mu(d) \]

暴力更换不太好枚举的一些东西（比如谁整除谁）

\[=\sum_{i=1}^a\sum_{j=1}^b\sum_{x|i}\sum_{y|j}\sum_{d=1}^{min(a,b)}\mu(d)[d|gcd(x,y)] \]

\[=\sum_{d=1}^{min(a,b)}\mu(d)\sum_{i=1}^a\sum_{j=1}^b\sum_{x|i}\sum_{y|j}[d|gcd(a,b)] \]

调整求和顺序

\[=\sum_{d=1}^{min(a,b)}\mu(d)\sum_{x=1}^a\sum_{y=1}^b[d|gcd(a,b)]\sum_{x|i}^a\sum_{y|j}^b 1 \]

\[=\sum_{d=1}^{min(a,b)}\mu(d)\sum_{x=1}^a\sum_{y=1}^b[d|gcd(a,b)]\lfloor\frac{a}{x}\rfloor\lfloor\frac{b}{y}\rfloor \]

某一项太不好弄了，通过更改枚举项拿掉

\[=\sum_{d=1}^{min(a,b)}\mu(d)\sum_{x=1}^{\lfloor\frac{a}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{b}{d}\rfloor}\lfloor\frac{a}{dx}\rfloor\lfloor\frac{b}{dy}\rfloor \]

发现求和项也带有下取整，预处理前缀和以后直接整除分块就可以了。

---

Code:

#include <cstdio>
#define ll long long
const int N=5e4;
int pri[N+10],mu[N+10],ispri[N+10],f[N+10],cnt,T,a,b;
void init()
{
    mu[1]=1;
    for(int i=2;i<=N;i++)
    {
        if(!ispri[i])
        {
            pri[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&pri[j]*i<=N;j++)
        {
            ispri[i*pri[j]]=1;
            if(i%pri[j]==0) break;
            else mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++)
    {
        mu[i]+=mu[i-1];
        for(int l=1,r;l<=i;l=r+1)
        {
            r=i/(i/l);
            f[i]+=i/l*(r-l+1);
        }
    }
}
int min(int x,int y){return x<y?x:y;}
int main()
{
    init();
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d",&a,&b);
        ll ans=0;
        for(int l=1,r;l<=min(a,b);l=r+1)
        {
            r=min(a/(a/l),b/(b/l));
            ans+=1ll*(mu[r]-mu[l-1])*f[a/l]*f[b/l];
        }
        printf("%lld\n",ans);
    }
    return 0;
}

---

2018.10.20