[BZOJ2301][HAOI2011]Problem b

[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
HINT
100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000

Solution
首先把区间问题转化为前缀问题，记答案为
$f(k)=\sum_{1<=x<=n, 1<=y<=m} [gcd(x, y) == k]$
我们记
$g(k)=\sum_{1<=x<=n, 1<=y<=m} [k|gcd(x,y)]$
$考虑一对(x,y)，如果被计数，则x=x'*k, y=y'*k$
所以我们可以发现
$g(k)$
$=\sum_{1<=x'<= \lfloor \frac{n}{k}\rfloor 1<=y'<= \lfloor \frac{m}{k}\rfloor}$ 1
$=\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor$
由两个函数的定义知道
$g(k) = \sum_{1 <= d <= \lfloor \frac{n}{k} \rfloor} f(d * k)=\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor$
由莫比乌斯反演的形式二知道
$f(k) = \sum_{1 <= d <= \lfloor \frac{n}{k} \rfloor} g(d * k) * \mu(d)$
$=\sum_{1 <= d <= \lfloor \frac{n}{k} \rfloor}\lfloor \frac{n}{k*d}\rfloor \lfloor \frac{m}{k*d}\rfloor * \mu(d)$
$令n' =\lfloor \frac{n}{k} \rfloor, m' = \lfloor \frac{m}{k} \rfloor$
$则有f(k) = \sum_{1 <= d <= n'}\lfloor \frac{n'}{d}\rfloor \lfloor \frac{m'}{d}\rfloor * \mu(d)$
利用分块可以根号时间解决

Code

#include <bits/stdc++.h>
using namespace std;
#define rep(i, l, r) for (int i = (l); i <= (r); i++)
#define per(i, r, l) for (int i = (r); i >= (l); i--)
#define X first
#define Y second
#define MS(_) memset(_, 0, sizeof(_))
#define MP make_pair
#define PB push_back
#define debug(...) fprintf(stderr, __VA_ARGS__)
template<typename T> inline void read(T &x){
    x = 0; T f = 1; char ch = getchar();
    while (!isdigit(ch)) {if (ch == '-') f = -1; ch = getchar();}
    while (isdigit(ch))  {x = x * 10 + ch - '0'; ch = getchar();}
    x *= f;
}

const int N = 55555;
int mu[N], check[N], prime[N], sum[N], K;

inline void getmu(){ int tot = 0;
    MS(check); mu[1] = 1;
    rep(i, 2, N){
        if (!check[i]){mu[i] = -1; prime[++tot] = i;}
        rep(j, 1, tot){
            if (i * prime[j] > N) break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0){ mu[i * prime[j]] = 0; break;
            }else mu[i * prime[j]] = -mu[i];
        }
    }
    rep(i, 1, N) sum[i] = sum[i-1] + mu[i];
}
inline int cal(int n, int m){ int ans = 0, pos  = 0;
    if (n > m) swap(n, m);
    n /= K; m /= K;
    for (int i = 1; i <= n; i = pos+1){
        pos = min(n/(n/i), m/(m/i));
        ans += (sum[pos] - sum[i-1]) * (n/i) * (m/i);
    }
    return ans;
}

int main(){
    getmu();    
    int case_cnt; read(case_cnt);
    while (case_cnt--){ int a, b, c, d;
        read(a); read(b); read(c); read(d); read(K);
        printf("%d\n", cal(b, d) - cal(b, c-1) - cal(a-1, d) + cal(a-1, c-1));
    }
    return 0;   
}