P3455

发表于 更新于 分类于
[POI2007]ZAP-Queries

[POI2007]ZAP-Queries

首先我们知道求的东西是：

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

可以化为：

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

这个东西显然可以莫比乌斯反演：

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

然后交换求和顺序：

$$\sum_{d=1}\mu(d)\sum_{i=1}^{\lfloor\frac{a}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{b}{k}\rfloor}[d\mid \gcd(i,j)]$$

显然这个式子就等价于：

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

又因为有：

$$\sum_{i=1}^{\lfloor\frac{a}{k}\rfloor}[d\mid i]=\lfloor\dfrac{a}{kd}\rfloor$$

上式可化为：

$$\sum_{d=1}\mu(d)\lfloor\dfrac{a}{kd}\rfloor\lfloor\dfrac{b}{kd}\rfloor$$

然后就是个二维数论分块。

注意不要过度使用除法，好像小卡常数。。。

时间复杂度 $O(n+t\sqrt{n})$。

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include<iostream>
#include<cstdio>
#define ll long long
using namespace std;

const ll N=5e4;

ll n,a,b,k,ans,cnt;

ll mu[N+5],s[N+5],prime[N+5];
bool f[N+5];

inline void Init() {
  mu[1]=1;f[1]=1;
  for(ll i=2;i<=N;i++) {
    if(!f[i]) {prime[++cnt]=i;mu[i]=-1;}
    for(ll j=1;j<=cnt&&prime[j]*i<=N;j++) {
      f[i*prime[j]]=1;
      if(i%prime[j]==0) {mu[i*prime[j]]=0;break;}
      mu[i*prime[j]]=-mu[i];
    }
  }
  for(ll i=1;i<=N;i++) {s[i]=s[i-1]+mu[i];}
}

inline ll read() {
  ll ret=0,f=1;char ch=getchar();
  while(ch<'0'||ch>'9') {if(ch=='-') f=-f;ch=getchar();}
  while(ch>='0'&&ch<='9') {ret=(ret<<3)+(ret<<1)+ch-'0';ch=getchar();}
  return ret*f;
}

inline void write(ll x) {
  static char buf[22];static ll len=-1;
  if(x>=0) {do{buf[++len]=x%10+48;x/=10;}while(x);}
  else {putchar('-');do{buf[++len]=-(x%10)+48;x/=10;}while(x);}
  while(len>=0) putchar(buf[len--]);
}

inline void writeln(ll x) {write(x);putchar('\n');}

int main() {

  Init();

  n=read();

  while(n--) {
    a=read();b=read();k=read();ans=0;
    ll tmp1=a/k,tmp2=b/k;
    for(ll i=1,j;i<=tmp1&&i<=tmp2;i=j+1) {
      j=min(tmp1/(tmp1/i),tmp2/(tmp2/i));
      ans+=(s[j]-s[i-1])*(tmp1/i)*(tmp2/i);
    }
    writeln(ans);
  }

  return 0;
}