[HAOI2011]Problem b

化成數學代數式就是這樣：

$$\sum_{x=a}^b\sum_{y=c}^d[\gcd(x,y)=k]$$

根據差分的思想，能夠將該式變成：

$$\sum_{x=1}^b\sum_{y=1}^d[\gcd(x,y)=k]-\sum_{x=1}^{a-1}\sum_{y=1}^d[\gcd(x,y)=k]-\sum_{x=1}^b\sum_{y=1}^{c-1}[\gcd(x,y)=k]+\sum_{x=1}^{a-1}\sum_{y=1}^{b-1}[\gcd(x,y)=k]$$

只要求出：

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

最後的答案就能求解。

咱們來推導這個式子：

$$\begin{array}{ll} &\sum\limits_{i=1}^N\sum\limits_{j=1}^M[\gcd(i,j)=k]\\=&\sum\limits_{i=1}^{\lfloor\frac Nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac Mk\rfloor}[\gcd(i,j)=1]\\ =&\sum\limits_{i=1}^{\lfloor\frac Nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac Mk\rfloor}\epsilon(\gcd(i,j)) & \text{根據單位函數ϵ的定義}\\=&\sum\limits_{i=1}^{\lfloor\frac Nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac Mk\rfloor}\sum\limits_{d|gcd(i,j)}\mu(d)&\text{注意這一步的變換是根據} \epsilon=\mu*1 \text{得來的}\\=&\sum\limits_{i=1}^{\lfloor\frac Nk\rfloor}\sum\limits_{j=1}^{\lfloor\frac Mk\rfloor}\sum\limits_{d|i,d|j}\mu(d)\\=&\sum\limits_{d=1}^{\min(N,M)}\mu(d)\sum\limits_{d|i,i\leq \lfloor\frac Nk\rfloor}\sum\limits_{d|j,j\leq \lfloor\frac Mk\rfloor}&\text{調換求和順序} \\=&\sum\limits_{d=1}^{\min(N,M)}\mu(d) \lfloor \frac{\lfloor \frac Nk \rfloor}d \rfloor \lfloor \frac{\lfloor \frac Mk \rfloor}d \rfloor \end{array}$$

如今能夠用$O(\min(N,M))$的複雜度解決，但還不夠。咱們要用一種科技：整除分塊。

有一個性質：一個數$N$整除任意數，商最多有$2\sqrt{N}$種不一樣的數；

證實：設除數$d$，當$d\leq \sqrt{N}$，最多有$\sqrt{N}$種除數；當$d\leq \sqrt{N}$，$0\leq\left\lfloor \dfrac Nd \right\rfloor\leq \sqrt{N}$最多隻有$\sqrt{N}$種，綜上最多有$2\sqrt{N}$種。

再一個性質：一個數$N$與另外一個正整數$l$整除爲$d$，在$d$不變的狀況下，$l$的最大值$l'$爲$\left\lfloor \dfrac N{\lfloor \frac Nl \rfloor} \right\rfloor$；

證實：回想下小學的東西：被除數必定，商越小，除數越大。由於是整除，因此商是整數，若是不是整除，商$d'$必定知足$d\leq d' <d+1$，因此除數最大時最小的商$d'=d=\left\lfloor\dfrac Nl\right\rfloor$。除數的最大值$l'$就是$\left\lfloor \dfrac N{\lfloor \frac Nl \rfloor} \right\rfloor$了（加下取整是由於$l'$是整數）。

根據上面兩個性質，算法就應運而生了。

代碼實現（這裏用$r$代替了上面的$l'$）：

for (int l = 1, r; l <= N; l = r+1) {
    r = N / (N / l);
    // do something
}

 複雜度爲$O(\sqrt{N})$。

而後咱們能夠經過分塊求出當商相等的$d$的區間再乘上$\mu$的前綴差分便可。

關於$\mu$的求解，由於是積性函數，能夠用線性篩篩出，複雜度$O(\text{size})$，而後直接前綴處理，最後分塊統計便可。

總體複雜度：$O(T\sqrt N)$，這裏$T$是數據組數，$N$爲$a,b,c,d$的範圍。詳見底下的代碼。

 1 #include <bits/stdc++.h>
 2 
 3 using namespace std;
 4 
 5 #define re register
 6 #define rep(i, a, b) for (re int i = a; i <= b; ++i)
 7 #define repd(i, a, b) for (re int i = a; i >= b; --i)
 8 #define srep(i, a, b, c) for (re int i = a; i <= b; c)
 9 #define maxx(a, b) a = max(a, b);
10 #define minn(a, b) a = min(a, b);
11 #define LL long long
12 #define INF (1 << 30)
13 
14 inline int read() {
15     int w = 0, f = 1; char c = getchar();
16     while (!isdigit(c)) f = c == '-' ? -1 : f, c = getchar();
17     while (isdigit(c)) w = (w << 3) + (w << 1) + (c ^ '0'), c = getchar();
18     return w * f;
19 }
20 
21 const int maxn = 5e4 + 5;
22 
23 int mu[maxn], vis[maxn], p[maxn], prime_tot;
24 
25 void get_mu() {
26     prime_tot = 0;
27     memset(vis, 0, sizeof(vis));
28     mu[1] = 1; int size = 50000;
29     rep(i, 2, size) {
30         if (!vis[i]) {
31             p[prime_tot++] = i;
32             mu[i] = -1;
33         }
34         rep(j, 0, prime_tot-1) {
35             if (i * p[j] > size) break;
36             vis[i * p[j]] = 1;
37             if (!(i % p[j])) {
38                 mu[i * p[j]] = 0;
39                 break;
40             }
41             mu[i * p[j]] = -mu[i];
42         }
43     }
44     rep(i, 2, size) mu[i] += mu[i-1];
45 }
46 
47 LL calc(int m, int n, int k) {
48     m /= k, n /= k;
49     int loop = min(m, n);
50     LL ans = 0;
51     for (int l = 1, r; l <= loop; l = r+1) {
52         r = min(n / (n / l), m / (m / l));
53         ans += 1ll * (mu[r]-mu[l-1]) * (1ll * (n/l)*(m/l));
54     }
55     return ans;
56 }
57 
58 int a, b, c, d, e, n;
59 
60 int main() {
61     n = read();
62     get_mu();
63     while (n--) {
64         a = read(), b = read(), c = read(), d = read(), e = read();
65         printf("%lld\n", calc(b, d, e) - calc(a-1, d, e) - calc(b, c-1, e) + calc(a-1, c-1, e));
66     }
67     return 0;
68 }