（轉載）關於gcd的8題

發現其實有關gcd的題目仍是挺多的，這裏根據作題順序寫出8題。

[bzoj2818: Gcd] gcd(x,y)=質數, 1<=x,y<=n的對數

作這題的時候，懂得了一個很是重要的轉化：求(x, y) = k, 1 <= x, y <= n的對數等於求(x, y) = 1, 1 <= x, y <= n/k的對數！因此，枚舉每一個質數p（線性篩素數的方法見：線性時間內篩素數和歐拉函數），而後求(x, y) = 1, 1 <= x, y <= n/p的個數。

那(x, y) = 1的個數如何求呢？其實就是求互質的數的個數。在[1, y]和y互質的數有phi(y)個，若是咱們令x < y，那麼答案就是sigma(phi(y))。由於x, y是等價的，因此答案*2，又由於(1, 1)只有一對，因此-1。最終答案爲sigma(sigma(phi(n/prime[i])) * 2 - 1)。

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

#include <cstdio> #include <algorithm> const int MAXN = 10000001; int n, primes; long long prime[MAXN], phi[MAXN]; bool com[MAXN]; int main() { scanf("%d", &n); for (int i = 2; i <= n; ++i) { if (!com[i]) { prime[primes++] = i; phi[i] = i-1; } for (int j = 0; j < primes && i*prime[j] <= n; ++j) { com[i*prime[j]] = true; if (i % prime[j]) phi[i*prime[j]] = phi[i] * (prime[j]-1); else { phi[i*prime[j]] = phi[i] * prime[j]; break; } } } phi[1] = 1; for (int i = 2; i <= n; ++i) phi[i] += phi[i-1]; long long ans = 0; for (int i = 0; i < primes; ++i) ans += phi[n/prime[i]]*2-1; printf("%lld", ans); }

[bzoj2190: [SDOI2008]儀仗隊] 求gcd(x,y)=1, 0<=x,y<=n的對數

嗯……彷佛這題在上一題的時候解決了。不過要注意的是，這題範圍是從0開始的，因此，會多出(1, 0)和(0, 1)這兩組，答案就是sigma(phi(i)) - 1 + 2。

[bzoj2005: [Noi2010]能量採集] 求sigma(gcd(x,y)), 1<=x<=n, 1<=y<=m

令f[d]爲(x, y) = d的對數，那麼答案就是sigma(f[i]*((i-1)*2+1))

f[i]怎麼求呢？在1 <= x <= n, 1 <= y <= m中，gcd(x, y) | d的有[n/d] * [m/d]個。不過咱們要扣掉全部的倍數，f[i] = [n/d] * [m/d] - f[2i] - f[3i] - f[4i] - ...。逆序作便可。

後來我在想：

• 爲何上面幾題不用這題的方法呢？嗯，由於x, y的取值不同，這樣的話就不得不把莫比烏斯反演搬出來了。
• 爲何這題不用上面幾題的方法呢？嗯，由於這個的複雜度是O(nlogn)（O(n/1 + n/2 + ... + n/n) = O(nlogn)），而上面的複雜度只有O(n + n/logn)（質數個數是n/logn級別的）。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

#include <cstdio> #include <algorithm> const int MAXN = 100001; int n, m; long long f[MAXN]; int main() { scanf("%d%d", &n, &m); if (n > m) std::swap(n, m); long long ans = 0; for (int i = n; i >= 1; --i) { f[i] = static_cast<long long>(n/i) * (m/i); for (int j = i+i; j <= n; j += i) f[i] -= f[j]; ans += f[i]*(2*i-1); } printf("%lld", ans); }

[SPOJ VLATTICE] gcd(i,j,k)=1, 0<=i,j,k<=n 的個數

莫比烏斯反演（Möbius inversion formula）終於出現了！莫比烏斯反演的基本形式就是：

g(n) = sigma(d|n, f(d))
f(n) = sigma(d|n, μ(n/d) * g(d))

還有另一個形式是：

g(n) = sigma(d|n, f(d))
f(n) = sigma(n|d, μ(n) * g(n/d))

通俗的來講，g(n)和f(n)的關係，就是包含和僅包含的關係。

令g(x)爲知足x | (i, j, k)的個數，f(x)爲知足(i, j, k) = x的個數。顯然g(n) = sigma(d|n, f(d)) = [n/x] * [n/x] * [n/x]，因此答案f(1)就能夠用下面那個公式算出來了，也就是sigma(μ(d) * [n/d] * [n/d] * [n/d])。

不過須要注意的是，這題的範圍能夠取0，因此咱們還須要加上某一維爲0的，某兩維爲0的方案（基本同樣啦）。

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

#include <cstdio> #include <algorithm> const int MAXN = 1000000; inline long long sqr(const long long x) { return x * x; } inline long long cube(const long long x) { return x * x * x; } int Testcase, n, primes, prime[MAXN+1], mu[MAXN+1]; void GetPrimes() { static bool com[MAXN+1]; mu[1] = 1; for (int i = 2; i <= MAXN; ++i) { if (!com[i]) { prime[primes++] = i; mu[i] = -1; } for (int j = 0; j < primes && i*prime[j] <= MAXN; ++j) { com[i*prime[j]] = true; if (i % prime[j]) mu[i*prime[j]] = -mu[i]; else { mu[i*prime[j]] = 0; break; } } } } int main() { GetPrimes(); for (scanf("%d", &Testcase); Testcase--; ) { scanf("%d", &n); long long ans = 3; for (int i = 1; i <= n; ++i) ans += mu[i] * cube(n/i); for (int i = 1; i <= n; ++i) ans += mu[i] * sqr(n/i) * 3; printf("%lld\n", ans); } }

[bzoj1101: [POI2007]Zap] 求有多少對數知足 gcd(x,y)=d, 1<=x<=a, 1<=y<=b

首先確定想到轉化成求gcd(x, y) = 1, 1 <= x <= a/d, 1 <= y <= b/d的對數，和上面一題同樣。雖然上面那題複雜度很優，只有O(n)，可是對於這邊的多組數據徹底就沒法承受了。

其實，只有O(sqrt(n))個不一樣的n/i的取值，所以咱們能夠求mu[i]的前綴和而後分塊作！分塊能夠參考[CQOI2007]餘數之和sum bzoj1257。

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

#include <cstdio> #include <algorithm> const int MAXN = 50000; inline long long sqr(const long long x) { return x * x; } int Testcase, n, m, d, primes, prime[MAXN+1], mu[MAXN+1], sum[MAXN+1]; void GetPrimes() { static bool com[MAXN+1]; mu[1] = 1; for (int i = 2; i <= MAXN; ++i) { if (!com[i]) { prime[primes++] = i; mu[i] = -1; } for (int j = 0; j < primes && i*prime[j] <= MAXN; ++j) { com[i*prime[j]] = true; if (i % prime[j]) mu[i*prime[j]] = -mu[i]; else { mu[i*prime[j]] = 0; break; } } } for (int i = 1; i <= MAXN; ++i) sum[i] = sum[i-1] + mu[i]; } int main() { GetPrimes(); for (scanf("%d", &Testcase); Testcase--; ) { scanf("%d%d%d", &n, &m, &d); if (n > m) std::swap(n, m); long long ans = 0; n /= d, m /= d; for (int i = 1, last; i <= n; i = last+1) { last = std::min(n/(n/i), m/(m/i)); ans += (sum[last]-sum[i-1]) * static_cast<long long>(n/i) * (m/i); } printf("%lld\n", ans); } }

[bzoj2301: [HAOI2011]Problem b] 求有多少對數知足 gcd(x,y)=k, a<=x<=b, c<=y<=d

和上面那題幾乎同樣，只是多了x, y的下界。其實這個能夠在算方案的時候yy一下，可是總以爲可能會哪裏邊界算錯，因而就是用了一種偷懶的方法：容斥原理。事實上，這種方法並不會比只算一次的方法來的慢多少。

記f[a, b]爲知足(x, y) = 1, 1 <= x <= a, 1 <= y <= b的對數，那麼答案就是f[B/K, D/K] - f[(A-1)/K, D/K] - f[B/K, (C-1)/K] + f[(A-1)/K, (C-1)/K]

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

#include <cstdio> #include <algorithm> const int MAXN = 50000; inline long long sqr(const long long x) { return x * x; } int Testcase, A, B, C, D, K, primes, prime[MAXN+1], mu[MAXN+1], sum[MAXN+1]; void GetPrimes() { static bool com[MAXN+1]; mu[1] = 1; for (int i = 2; i <= MAXN; ++i) { if (!com[i]) { prime[primes++] = i; mu[i] = -1; } for (int j = 0; j < primes && i*prime[j] <= MAXN; ++j) { com[i*prime[j]] = true; if (i % prime[j]) mu[i*prime[j]] = -mu[i]; else { mu[i*prime[j]] = 0; break; } } } for (int i = 1; i <= MAXN; ++i) sum[i] = sum[i-1] + mu[i]; } long long Process(int n, int m) { long long res = 0; if (n > m) std::swap(n, m); for (int i = 1, last; i <= n; i = last+1) { last = std::min(n/(n/i), m/(m/i)); res += (sum[last]-sum[i-1]) * static_cast<long long>(n/i) * (m/i); } return res; } int main() { GetPrimes(); for (scanf("%d", &Testcase); Testcase--; ) { scanf("%d%d%d%d%d", &A, &B, &C, &D, &K); long long ans = Process(B/K, D/K); ans -= Process((A-1)/K, D/K); ans -= Process(B/K, (C-1)/K); ans += Process((A-1)/K, (C-1)/K); printf("%lld\n", ans); } }

[bzoj2820: YY的GCD] 求1<=x<=N, 1<=y<=M且gcd(x,y)爲質數有多少對

若是枚舉質數的話就要給這個多組數據跪成傻逼了。

ans = sigma(p, sigma(d, μ(d) * (n/pd) * (m/pd)))

Let s = pd, then

ans = sigma(s, sigma(p, μ(s/p) * (n/s) * (m/s)))
    = sigma(s, (n/s) * (m/s) * sigma(p, μ(s/p)))

Let g(x) = sigma(p, μ(x/p)), then

ans = sigma(s, (n/s) * (m/s) * g(s))

若是咱們能預處理g(x)的話就能和前面同樣分塊搞了。這個時候咱們多麼但願g(x)和μ(x)同樣是積性函數。看完題解後，發現有一個不是積性函數，勝似積性函數的性質。因爲題解沒有給證實，因此就意淫了一個證實。

考慮質數p'，g(p'x) = sigma(p | p'x, μ(p'x/p))。

• 當x % p' == 0，那麼考慮sigma中的變量p的全部取值，它和g(x)中p是相同的。而μ(x)這個函數，若是x有平方因子的話就等於0，所以當p != p'時μ(p'x/p) = 0，當p == p'是，μ(p'x/p) = μ(x)。因此g(p'x) = μ(x)。
• 當x % p' != 0，一樣考慮p，會發現它的取值只比g(x)中的p多出一個p'。同理按照p是否等於p'討論，能夠獲得g(p'x) = -g(x) + μ(x)。

所以g(x)能夠在線性篩素數的時候算出。剩下的就是前綴和、分塊了。

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

#include <cstdio> #include <algorithm> const int MAXN = 10000000; int Testcase, n, m, primes, prime[MAXN+1], mu[MAXN+1], g[MAXN+1], sum[MAXN+1]; void GetPrimes() { static bool com[MAXN+1]; mu[1] = 1; for (int i = 2; i <= MAXN; ++i) { if (!com[i]) prime[primes++] = i, mu[i] = -1, g[i] = 1; for (int j = 0; j < primes && i*prime[j] <= MAXN; ++j) { com[i*prime[j]] = true; if (i % prime[j]) mu[i*prime[j]] = -mu[i], g[i*prime[j]] = -g[i] + mu[i]; else { mu[i*prime[j]] = 0, g[i*prime[j]] = mu[i]; break; } } } for (int i = 1; i <= MAXN; ++i) sum[i] = sum[i-1] + g[i]; } int main() { GetPrimes(); for (scanf("%d", &Testcase); Testcase--; ) { scanf("%d%d", &n, &m); if (n > m) std::swap(n, m); long long ans = 0; for (int i = 1, last; i <= n; i = last+1) { last = std::min(n/(n/i), m/(m/i)); ans += (n/i) * static_cast<long long>(m/i) * (sum[last]-sum[i-1]); } printf("%lld\n", ans); } }

[bzoj2705: [SDOI2012]Longge的問題] 求 sigma(gcd(i,n), 1<=i<=n<2^32)

又是令gcd(i, n) = d，答案就是sigma(phi(n/d))，可是咱們不能預處理出phi[]數組，由於開不了數組……

注意到因數個數是O(2sqrt(n))級別的，咱們枚舉全部的n/d，一邊dfs一邊算phi。

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

#include <cmath> #include <cstdio> int n, cnt, p[30], c[30]; long long ans = 0; void dfs(const int step, int pdt, int phi) { if (step == cnt) { ans += phi; return; } dfs(step+1, pdt, phi); phi = phi / p[step] * (p[step]-1); for (int i = 1; i <= c[step]; ++i) dfs(step+1, pdt *= p[step], phi); } int main() { scanf("%d", &n); int x = n; for (int i = 2; i*i <= x; ++i) if (x % i == 0) { for (; x % i == 0; x /= i) ++c[cnt]; p[cnt++] = i; } if (x > 1) c[cnt] = 1, p[cnt++] = x; dfs(0, 1, n); printf("%lld\n", ans); }