题意:
这道题和POJ 3090很相似,求|x|≤a,|y|≤b 中站在原点可见的整点的个数K,所有的整点个数为N(除去原点),求K/N
分析:
坐标轴上有四个可见的点,因为每个象限可见的点数都是一样的,所以我们只要求出第一象限可见的点数然后×4+4,即是K。
可见的点满足gcd(x, y) = 1,于是将问题转化为x∈[1, a], y∈[1, b],求gcd(x, y) = 1的个数。
类比HDU 1695可以用莫比乌斯反演来做,我还写了普通的和分块加速的两份代码,交上去发现运行时间相差并不是太多。
1 #include <cstdio> 2 #include <algorithm> 3 typedef long long LL; 4 5 const int maxn = 2000; 6 int mu[maxn+10], vis[maxn+10], prime[maxn]; 7 8 void Mobius() 9 { 10 mu[1] = 1; 11 int cnt = 0; 12 for(int i = 2; i <= maxn; ++i) 13 { 14 if(!vis[i]) 15 { 16 mu[i] = -1; 17 prime[cnt++] = i; 18 } 19 for(int j = 0; j < cnt && (LL)i*prime[j] <= maxn; ++j) 20 { 21 vis[i*prime[j]] = 1; 22 if(i % prime[j] != 0) mu[i*prime[j]] = -mu[i]; 23 else 24 { 25 mu[i*prime[j]] = 0; 26 break; 27 } 28 } 29 } 30 //计算前缀和,用于分块加速 31 for(int i = 2; i <= 2015; ++i) mu[i] += mu[i-1]; 32 33 } 34 35 int main() 36 { 37 Mobius(); 38 int a, b; 39 while(scanf("%d%d", &a, &b) == 2) 40 { 41 if(a == 0 && b == 0) break; 42 LL K = 0, N = (LL)(a*2+1) * (b*2+1) - 1; 43 if(a > b) std::swap(a, b); 44 for(int i = 1, j; i <= a; i = j + 1) 45 { 46 j = std::min(a/(a/i), b/(b/i)); 47 K += (LL)(mu[j] - mu[i-1]) * (a/i) * (b/i); 48 } 49 //for(int i = 1; i <= a; ++i) K += (LL) mu[i] * (a/i) * (b/i); 50 K = (K+1)*4; 51 52 printf("%.7f\n", (double) K / N); 53 } 54 55 return 0; 56 }
代码君
也可以按照紫书上的思路求欧拉函数,因为a的范围比较小,所以可以逐列统计。不过这个方法要比上面的莫比乌斯反演慢得多
1 #include <cstdio> 2 #include <cmath> 3 typedef long long LL; 4 5 int phi(int n) 6 { 7 int m = sqrt(n + 0.5); 8 int ans = n; 9 for(int i = 2; i <= m; ++i) if(n % i == 0) 10 { 11 ans = ans / i * (i-1); 12 while(n % i == 0) n /= i; 13 } 14 if(n > 1) ans = ans / n * (n-1); 15 return ans; 16 } 17 18 int gcd(int a, int b) 19 { 20 return b == 0 ? a : gcd(b, a%b); 21 } 22 23 int main() 24 { 25 int a, b; 26 while(scanf("%d%d", &a, &b) == 2 && a) 27 { 28 LL N = (LL)(a*2+1) * (b*2+1) - 1; 29 LL K = 0; 30 for(int x = 1; x <= a; ++x) 31 { 32 int k = b / x; 33 K += phi(x) * k; 34 for(int y = k*x+1; y <= b; y++) 35 if(gcd(x, y) == 1) K++; 36 } 37 K = (K+1)*4; 38 printf("%.7f\n", (double)K / N); 39 } 40 41 return 0; 42 }
代码君二
时间: 2024-08-08 01:16:28