「HAOI 2011」Problem B

Description

题目链接：BZOJ 2301

求如下式子的值：

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

本题 $T$ 组数据。

数据范围：$1\le T,x,y,n,m,k\le 5\times 10^4$

---

Solution

根据容斥原理，原式可以分成 $4$ 块来处理，每一块的式子形式都为

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

考虑化简该式子

$$\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)=1]$$

因为 $\gcd(i,j)=1$ 时对答案才用贡献，于是我们可以将其替换为 $\varepsilon(\gcd(i,j))$（$\varepsilon(n)$ 当且仅当 $n=1$ 时值为 $1$ 否则为 $0$ ），故原式化为

$$\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\varepsilon(\gcd(i,j))$$

将 $\varepsilon$ 函数展开得到

$$\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d\mid \gcd(i,j)}\mu(d)$$

变换求和顺序，先枚举 $d\mid gcd(i,j)$ 可得（其中 $[d\mid i]$ 表示 $i$ 是 $d$ 的倍数时对答案有 $1$ 的贡献）

$$\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}[d\mid i]\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[d\mid j]$$

易知 $1\sim\lfloor\frac{n}{k}\rfloor$ 中 $d$ 的倍数有 $\lfloor\frac{n}{kd}\rfloor$ 个，故原式化为

$$\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d) \lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor$$

很显然，式子可以数论分块求解（过程中默认 $n\le m$）。

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

---

Code

#include <cstdio>
#include <algorithm>
const int N=50000;
int mu[N+5],p[N+5];
bool flg[N+5];
void init() {
    int tot=0;
    mu[1]=1;
    for(int i=2;i<=N;++i) {
        if(!flg[i]) {
            p[++tot]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=tot&&i*p[j]<=N;++j) {
            flg[i*p[j]]=1;
            if(i%p[j]==0) {
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;++i) mu[i]+=mu[i-1];
}
int solve(int n,int m) {
    int res=0;
    for(int i=1,j;i<=std::min(n,m);i=j+1) {
        j=std::min(n/(n/i),m/(m/i));
        res+=(mu[j]-mu[i-1])*(n/i)*(m/i);
    }
    return res;
}
int main() {
    int T,a,b,c,d,k;
    init();
    for(scanf("%d",&T);T;--T) {
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        printf("%d\n",solve(b/k,d/k)-solve(b/k,(c-1)/k)-solve((a-1)/k,d/k)+solve((a-1)/k,(c-1)/k));
    }
    return 0;
}