【 HAOI2011】 Problem b （数论，分块，莫比乌斯函数）

· · 题解

题意

对于给出的n个询问，每次求有多少个数对(x,y)，满足a≤x≤b，c≤y≤d，且gcd(x,y) = k，gcd(x,y)函数为x和y的最大公约数。

分析

所求 ：\sum_{i=a}^b \sum_{j=c}^d [gcd(i,j)==k]

二维区间和的问题通常会转化为二维前缀和来容斥

所求： solve(n,m) = \sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)==k]

为了同一推导，默认n\leq m ，如果 n >m可以交换n,m,因为solve(n,m) == solve(m,n)

我们常常研究互质的情况和性质，如果互质就容易联想到一系列积性函数

变形：\sum_{i=1}^{\frac n k} \sum_{j=1}^{\frac m k} [gcd(i,j)==1]

联想到\epsilon (n) 也是在gcd == 1 的时候取值为1,其他情况取值为0

变形：\sum_{i=1}^{\frac n k} \sum_{j=1}^{\frac m k} \epsilon (gcd(i,j))

因为 1* \mu = \epsilon （这里*指狄利克雷卷积）

变形：\sum_{i=1}^{\frac n k} \sum_{j=1}^{\frac m k} \sum_{d|gcd(i,j)} \mu(d)

交换求和顺序：通常可以从意义出发来推，实在不会可以展开找规律,在我们想要把\mu(d)提到最前面，便于计数 参考另一篇博客的和式知识部分

因为gcd(i,j) <= min(i,j)

变形：\sum_{d=1}^{\frac n k} \mu(d) \sum_{i=1}^{\frac n k}[d|i] \sum_{i=1}^{\frac m k}[d|j]

考虑：已知d, 在[1,n] 中，共有 \lfloor \frac n d\rfloor 个数满足能整除于d，

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

这里我们变成了可以O(n)求解的式子

因为向下取整其实是阶梯函数，会有很多连续相等的部分，我们可以考虑聚合连续相等的部分一起算（也就是数论分块）

考虑含有\lfloor \frac n i\rfloor 的式子（ｎ为常数）

对于任意一个i 我们需要找到一个最大的j 使得\lfloor \frac n i \rfloor = \lfloor \frac n j \rfloor

则j = \lfloor \frac n {\lfloor \frac n i \rfloor} \rfloor

利用上述结论，每次以[i,j]为一块

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5+5;
bool not_prime[N];
int read() {
    int x=0;char ch = getchar();
    while(!isdigit(ch)) ch=getchar();
    while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
    return x;
}
int prime[N],cnt,mu[N];
int sum[N];
void mus(int n) {
    mu[1] = 1;
    cnt=0;
    for(int i=2;i<=n;i++) {
        if(!not_prime[i]) {
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for(int j=1;j<=cnt && i*prime[j] <= n;j++) {
            not_prime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                mu[i*prime[j]] = 0;
                break;
            } 
            else{
                mu[i*prime[j]] = -mu[i];
            }
        }
    }
    sum[0]=0;
    for(int i=1;i<=n;i++) {
        sum[i] = sum[i-1] + mu[i];
    }
}
int solve(int n,int m,int k) {
    if(n < k || m < k ) return 0;
    if(n > m) swap(n,m);
    int nn = n/k;
    int mm = m/k;
    int res=0;
    for(int d=1;d<=nn;) {
        int enn = nn/(nn/d);
        int enm = mm/(mm/d);
        int en = min(enn,enm);
        res += (sum[en]-sum[d-1])*(nn/d)*(mm/d);
        d = en+1;
    }

    return res;
}
int main() {
    mus(100000);
    int _=read();
    while(_--) {
        int a=read(),b=read(),c=read(),d=read(),k=read();
        int ans = solve(b,d,k) - solve(a-1,d,k) - solve(b,c-1,k) + solve(a-1,c-1,k);
        printf("%d\n",ans);
    }
}