UVa 10214 (莫比乌斯反演 or 欧拉函数) Trees in a Wood.

 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 }

 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 }