bzoj 2301[HAOI2011]Problem b

Description

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

Input

第一行一个整数n，接下来n行每行五个整数，分别表示a、b、c、d、k

Output

共n行，每行一个整数表示满足要求的数对(x,y)的个数

Sample Input

2
2 5 1 5 1
1 5 1 5 2

Sample Output

14
3

最近在学习莫比乌斯反演，终于看到一道裸题可以让我A掉了
询问 a≤x≤b，c≤y≤d，且gcd(x,y) = k 的个数
即$F(a, b, c, d) = sum_a^b{x}sum_c^d{y}[gcd(x, y) == k]$

很明显我们可以用容斥转换成
$F(a, b, c, d)= F(1, b, 1, d) - F(1, b, 1, c - 1) - F(1, a - 1, 1, d) + F(1, a - 1, 1, c - 1)$

所以只需要求出:

$f(k, b, d) = sum_1^b{x}sum_1^d{y}[gcd(x, y) == k]$

通过莫比乌斯反演:

$f(k, b, d) = sum_ {k | t} g(t)μ(t / k)$

所以

$f(k, b, d) = sum_ {k | x}μ(frac{x}{k}){[frac{b}{x}]}[frac{d}{x}]$

$ = sum_ {i = 1}^{[frac{n}{k}]}μ(i){[frac{frac{b}{k}}{i}]}[frac{frac{d}{k}}{i}]$

然后预处理μ(n)的前缀和，后面可以$sqrt{n}$ 算出
总复杂度 $nsqrt{n}$

 1 #include <iostream>
 2 #include <cstdio>
 3 #include <cstring>
 4 #include <cmath>
 5 #include <algorithm>
 6 
 7 #define LL long long
 8 
 9 using namespace std;
10 
11 const int MAXN = 5e4 + 10;
12 int N;
13 int t1, t2;
14 int a, b, k, c, d;
15 int cnt = 0;
16 bool f[MAXN];
17 LL sm[MAXN];
18 int prime[MAXN], mul[MAXN];
19 inline LL read()
20 {
21     LL x = 0, w = 1; char ch = 0;
22     while(ch < '0' || ch > '9') {
23         if(ch == '-') {
24             w = -1;
25         }
26         ch = getchar();
27     }
28     while(ch >= '0' && ch <= '9') {
29         x = x * 10 + ch - '0';
30         ch = getchar();
31     }
32     return x * w;
33 }
34 
35 void init()
36 {
37     mul[1] = 1;
38     for(int i = 2; i <= 50000; i++) {
39         if(!f[i]) {
40             prime[++cnt] = i;
41             mul[i] = -1;
42         } 
43         for(int j = 1; j <= cnt && prime[j] * i <= 50000; j++) {
44             f[prime[j] * i] = 1;
45             if(i % prime[j] == 0) {
46                 mul[prime[j] * i] = 0;
47                 break;
48             }
49             mul[prime[j] * i] = mul[i] * (-1);
50         }
51     }
52     for(int i = 1; i <= 50000; i++) {
53         sm[i] = sm[i - 1] + mul[i];
54     }
55 }
56 int main()
57 {
58     init();
59     
60     N = read();
61     while(N--) {
62         long long ans = 0;
63         a = read(), b = read(), c = read(), d = read(), k = read();
64         t1 = b, t2 = d;
65         if(t1 > t2) 
66             swap(t1, t2);
67         t1 = t1 / k, t2 = t2 / k;
68         for(int i = 1, j; i <= t1; i = j + 1) {
69             j = min(t1 / (t1 / i), t2 / (t2 / i));
70             ans += (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i);
71         }
72         t1 = a - 1, t2 = d;
73         if(t1 > t2) 
74             swap(t1, t2);
75         t1 = t1 / k, t2 = t2 / k;
76         for(int i = 1, j; i <= t1; i = j + 1) {
77             j = min(t1 / (t1 / i), t2 / (t2 / i));
78             ans -= (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i);
79         }
80         t1 = b, t2 = c - 1;
81         if(t1 > t2) 
82             swap(t1, t2);
83         t1 = t1 / k, t2 = t2 / k;
84         for(int i = 1, j; i <= t1; i = j + 1) {
85             j = min(t1 / (t1 / i), t2 / (t2 / i));
86             ans -= (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i);
87         }
88         t1 = a - 1, t2 = c - 1;
89         if(t1 > t2) 
90             swap(t1, t2);
91         t1 = t1 / k, t2 = t2 / k;
92         for(int i = 1, j; i <= t1; i = j + 1) {
93             j = min(t1 / (t1 / i), t2 / (t2 / i));
94             ans += (sm[j] - sm[i - 1]) * (t1 / i) * (t2 / i);
95         }
96         printf("%lld
", ans);
97     }
98     return 0;
99 }

View Code