【BZOJ】T1101 Zap_zap bzoj-CSDN博客

本文详细解析了求解给定范围内所有i和j的组合中，使得gcd(i,j)==d的个数的问题。通过数学推导，利用莫比乌斯函数的性质简化了问题，并给出了具体的实现代码。

题意：

k次询问，每次询问 $\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)==d$ ，即gcd(i,j)结果为d的个数。
k,n,m≤5e5。

题解：

原式显然可以化简为 $\sum_{i=1}^{N}\sum_{j=1}^{M}gcd(i,j)==1$ ，其中 $N=\lfloor \frac n d \rfloor,M=\lfloor \frac m d \rfloor$ 。
由于莫比乌斯函数有一个性质：

\sum d ∣ n μ (d) = {1, n = 1. 0, n \neq 1.

$\sum_{d \mid n} \mu(d) = \left \lbrace _ {0,n\not=1.} ^ {1,n=1.} \right .$

n=1时，显然原式值为1；
n≠1时，设 $d=\prod_{1}^{k} {p_i}^{a_i}$ ，p为质数，a≠0，讨论每个d对答案的贡献。
- $\exists a_i \ge 2$ ，d对答案的贡献=0；
- $\forall a_i=1$ ，d对答案的贡献= $\binom k 0-\binom k 1+\binom k 2-\binom k 3+\ldots$
根据二项式定理，当d不为1时对答案的贡献均为0，所以性质得证。

所以原式即可化为

\sum i = 1 N \sum j = 1 M \sum d ∣ g c d (i, j) μ (d)

$\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{d\mid gcd(i,j)}\mu(d)$
优先枚举d的取值，考虑

d∣gcd(i,j)d∣gcd(i,j) $d\mid gcd(i,j)$ ，可得

(d∣i)∧(d∣j)(d∣i)∧(d∣j) $(d\mid i) \land (d\mid j)$ ；
即

μ(d)μ(d) $\mu(d)$ 被统计进答案，当且仅当

i,ji,j $i,j$ 均为

dd $d$ 的整数倍
所以答案为

\sum_{d = 1}^{m i n (N, M)} μ (d) \times ⌊ \frac{N}{d} ⌋ \times ⌊ \frac{M}{d} ⌋

$\sum_{d=1}^{min(N,M)}\mu(d)\times\lfloor \frac N d \rfloor\times\lfloor \frac M d \rfloor$
由于

⌊Nd⌋⌊Nd⌋ $\lfloor \frac N d \rfloor$ ,

⌊Md⌋⌊Md⌋ $\lfloor \frac M d \rfloor$ 取值的阶梯性，我们可以在

O(n−−√)O(n) $O(\sqrt n)$ 的时间复杂度下枚举

⌊Nd⌋⌊Nd⌋ $\lfloor \frac N d \rfloor$ ,

⌊Md⌋⌊Md⌋ $\lfloor \frac M d \rfloor$ 的不同取值。
总时间复杂度

O(nn−−√)O(nn) $O(n\sqrt n)$ 。

#include<cstdio>
#include<algorithm>
#define N 55000
#define mx 50000
using namespace std;
int mobius[N],prime[N],flag[N],tot=0,sum[N];
void sol(){
    mobius[1]=1;
    for(int i=2;i<=mx;i++){
        if(!flag[i]){prime[tot++]=i;mobius[i]=-1;}
        for(int j=0;j<tot;j++){
            if(prime[j]*i>50000)break;
            flag[prime[j]*i]=true;
            if(i%prime[j]==0){mobius[prime[j]*i]=0;break;}
            else mobius[prime[j]*i]=-mobius[i];
        }
    }
    for(int i=1;i<=mx;i++)sum[i]=mobius[i]+sum[i-1];
    return;
}
int main(){
    int n;
    scanf("%d",&n);
    sol();
    while(n--){
        int n,m,d;
        scanf("%d%d%d",&n,&m,&d);
        n/=d,m/=d;
        //sigma(mobius(d)*(n/d)*(m/d))
        int ans=0,t;
        for(int i=1;(n/i)&&(m/i);i=t+1){
            t=min(n/(n/i),m/(m/i));
            ans+=(sum[t]-sum[i-1])*(n/i)*(m/i);
        }
        printf("%d\n",ans);
    }
}