spoj—7001Visible Lattice Points（莫比乌斯）-优快云博客

本文详细介绍了如何使用莫比乌斯反演公式解决特定的组合计数问题，即求解满足最大公约数为1的三元组(a,b,c)的数量，其中a、b、c均不超过给定的整数n。文章提供了两种求解莫比乌斯函数的方法，并附带了完整的代码实现。

题意：给定n(1<= n <= 1000000)，求满足gcd（a, b, c） = 1 (0 <= a, b, c <= n) 的(a, b, c)的个数。

分析：先说一下莫比乌斯反演：

莫比乌斯反演是用来简化运算的：如果F(n) = sigma(f(d)), d | n,F(n)已知，如何求f(n)?

莫比乌斯反演公式：f(k) = sigma(mu(d) * F(k / d)), d | k。

更常用的形式(在一定范围内) ：f(k) = sigma(mu(d | k) * F(d)), d | k;

其中mu(k) = 1 (k = 1)

(-1)^t (k = p1 *p2 * …… *pt, pi互不相同)

0 （k = 其他）

记g(k) 为gcd(a, b, c) = k 的（a, b, c）的个数，f(k) 为 k | gcd(a, b, c)的（a, b, c）的个数。

则f(k) = sigma(g(d)) ,d | k,

由莫比乌斯反演可知g(k) = sigma(mu[d / k] * f(d)), k | d。

f(k)是容易算出来的：f(k) = (n / k) * (n / k) * (n / k)

ans = g(1) = sigma(mu[d / 1] * f(d)), 1 | d, 即d = 1, 2 ……n

 1 #include <stdio.h>
 2 #include <string.h>
 3 #define maxn 1000010
 4 int mu[maxn], prim[maxn];
 5 
 6 void moblus()
 7 {
 8      memset(prim, 0, sizeof(prim));
 9      for(int i = 1; i < maxn; i++){
10           mu[i] = 1;
11      }
12      for(int i = 2; i < maxn; i++)
13      {
14           if(prim[i] == 0){
15                for(int j = i; j < maxn; j += i)
16                {
17                     prim[j] = 1;
18                     if(j % (i * i) == 0){
19                          mu[j] = 0;
20                     }else{
21                          mu[j] *= -1;
22                     }
23                }
24           }
25      }
26      return ;
27 }
28 
29 int main()
30 {
31      int n, T;
32      moblus();
33      scanf("%d", &T);
34      while(T--){
35           scanf("%d", &n);
36           long long ans = 3;
37           for(int i = 1; i <= n; i++){
38                ans += (long long)mu[i] * (n / i) * (n / i) * ((n / i) + 3);
39           }
40           printf("%lld\n", ans);
41      }
42      return 0;
43 }
44

View Code

上面代码求mu数组时间复杂度为O(n * log(n))，当然我们也可以线性求得：

 1 void moblus()
 2 {
 3      memset(check, 0, sizeof(check));
 4      int cnt = 0;
 5      mu[1] = 1;
 6      for(int i = 2; i < maxn; i++)
 7      {
 8           if(!check[i]){
 9                prim[++cnt] = i;
10                mu[i] = -1;
11           }
12           for(int j = 1; j <= cnt; j++)
13           {
14                if(i * prim[j] >= maxn){
15                     break ;
16                }
17                check[i * prim[j]] = 1;
18                if(i % prim[j] == 0){
19                     mu[i * prim[j]] = 0;
20                     break ;
21                }else{
22                     mu[i * prim[j]] = -mu[i];
23                }
24           }
25      }
26      return ;
27 }