题意:
这道题和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可以用莫比乌斯反演来做,我还写了普通的和分块加速的两份代码,交上去发现运行时间相差并不是太多。
#include <cstdio>
#include <algorithm>
typedef long long LL; const int maxn = ;
int mu[maxn+], vis[maxn+], prime[maxn]; void Mobius()
{
mu[] = ;
int cnt = ;
for(int i = ; i <= maxn; ++i)
{
if(!vis[i])
{
mu[i] = -;
prime[cnt++] = i;
}
for(int j = ; j < cnt && (LL)i*prime[j] <= maxn; ++j)
{
vis[i*prime[j]] = ;
if(i % prime[j] != ) mu[i*prime[j]] = -mu[i];
else
{
mu[i*prime[j]] = ;
break;
}
}
}
//计算前缀和,用于分块加速
for(int i = ; i <= ; ++i) mu[i] += mu[i-]; } int main()
{
Mobius();
int a, b;
while(scanf("%d%d", &a, &b) == )
{
if(a == && b == ) break;
LL K = , N = (LL)(a*+) * (b*+) - ;
if(a > b) std::swap(a, b);
for(int i = , j; i <= a; i = j + )
{
j = std::min(a/(a/i), b/(b/i));
K += (LL)(mu[j] - mu[i-]) * (a/i) * (b/i);
}
//for(int i = 1; i <= a; ++i) K += (LL) mu[i] * (a/i) * (b/i);
K = (K+)*; printf("%.7f\n", (double) K / N);
} return ;
}
代码君
也可以按照紫书上的思路求欧拉函数,因为a的范围比较小,所以可以逐列统计。不过这个方法要比上面的莫比乌斯反演慢得多,试验了下,不管是打表还是单独计算时间都差不多
#include <cstdio>
#include <cmath>
typedef long long LL; int phi(int n)
{
int m = sqrt(n + 0.5);
int ans = n;
for(int i = ; i <= m; ++i) if(n % i == )
{
ans = ans / i * (i-);
while(n % i == ) n /= i;
}
if(n > ) ans = ans / n * (n-);
return ans;
} int gcd(int a, int b)
{
return b == ? a : gcd(b, a%b);
} int main()
{
int a, b;
while(scanf("%d%d", &a, &b) == && a)
{
LL N = (LL)(a*+) * (b*+) - ;
LL K = ;
for(int x = ; x <= a; ++x)
{
int k = b / x;
K += phi(x) * k;
for(int y = k*x+; y <= b; y++)
if(gcd(x, y) == ) K++;
}
K = (K+)*;
printf("%.7f\n", (double)K / N);
} return ;
}
代码君二