题意:
这道题和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 }