码迷,mamicode.com
首页 > 其他好文 > 详细

HDU 4746 (莫比乌斯反演) Mophues

时间:2014-12-20 14:16:33      阅读:199      评论:0      收藏:0      [点我收藏+]

标签:

这道题看巨巨的题解看了好久,好久。。

本文转自hdu4746(莫比乌斯反演)

题意:给出n, m, p,求有多少对a, b满足gcd(a, b)的素因子个数<=p,(其中1<=a<=n, 1<=b<=m)

分析:设A(d):gcd(a, b)=d的有多少种

     设B(j): gcd(a, b)是j的倍数的有多少种,易知B(j) = (n/j)*(m/j)

     则由容斥原理得:(注:不同行的μ是不相同的,μ为莫比乌斯函数)

     A(1) = μ(1)*B(1) + μ(2)*B(2) + μ(3)*B(3) + ... + μ(p1*p2...)*B(p1*p2...)

     A(2) = μ(1)*B(1*2) + μ(2)*B(2*2) + μ(3)*B(3*2) + ... + μ(p1*p2..)*B(p1*p2..*2)

     ...

     A(d) = μ(1)*B(1*d) + μ(2)*B(2*d) + μ(3)*B(3*d) + ... + μ(p1*p2..)*B(p1*p2..*d)

 

     ans = A(1)+A(2)+...+A(d) = F(1)*B(1) + F(2)*B(2) + ... + F(p1*p2..)*B(p1*p2..)

     于是可以枚举公约数i{表示A(i)},利用筛法找出i的倍数j,i对B(j)的贡献系数为:F(j)+=μ(j/i)

     总之,求出B(j)的总贡献系数F(j)即可得答案:F(1)*B(1)+F(2)*B(2)+...+F(n)*B(n)

     上面没有限制gcd的素因子个数,要限制其实不难,给系数加多一维即可:

     F(d)(p)表示:素因子个数<=p时,对B(d)的贡献系数

   

     分块加速思想

     你可以再纸上模拟一下:设d在[i, n/(n/i)]的区间上,则该区间内所有的n/d都是一样的。

另外我想再补充一下:

  • 根据上面A、B的含义,有技术分享,然后根据莫比乌斯反演公式,反解出A(n),得技术分享
  • 代码中const int N = 19;以及特判if(p >= N)。为什么要定义N为19呢,因为如果一个正整数的素因子的个数大于等于19的话,那么这个数一定要比5×105要大,因为素因子个数为19的最小整数为219>5×105
  • 分块加速还想再啰嗦两句,因为技术分享,在计算B(i)时可以不用枚举每个i计算B(i)。举个栗子,技术分享。正如上面所说,d在区间[i, n/(n/i)]中,所有n/d的值都是一样的。这样就避免了重复计算B(i),在计算答案的时候预处理F(i, p)的前缀和即可。
技术分享
 1 #include <cstdio>
 2 #include <algorithm>
 3 typedef long long LL;
 4 
 5 const int M = 500000 + 10;
 6 const int N = 19;
 7 int F[M][N], num[M], h[M];//num记录素因子的个数,h如果含平凡因子则为-1,否则记录素因子的种类
 8 
 9 int Mob(int n)
10 {
11     if(h[n] == -1) return 0;
12     if(h[n] & 1) return -1;
13     return 1;
14 }
15 
16 void Init()
17 {
18     for(int i = 2; i < M; ++i)
19     {
20         if(num[i]) continue;
21         for(int j = i; j < M; j += i)
22         {
23             int cnt = 0, temp = j;
24             while(temp % i == 0)
25             {
26                 cnt++;
27                 temp /= i;
28             }
29             num[j] += cnt;
30             if(cnt > 1) h[j] = -1;
31             else if(h[j] >= 0) ++h[j];
32         }
33     }
34 
35     for(int i = 1; i < M; ++i)
36         for(int j = i; j < M; j += i)
37             F[j][num[i]] += Mob(j / i);
38     //求j的前缀和,使F表示素因子个数<=j的含义
39     for(int i = 1; i < M; ++i)
40         for(int j = 1; j < N; ++j)
41             F[i][j] += F[i][j-1];
42     //求i的前缀和,用于分块加速
43     for(int i = 1; i < M; ++i)
44         for(int j = 0; j < N; ++j)
45             F[i][j] += F[i-1][j];
46 }
47 
48 int main()
49 {
50     Init();
51 
52     int T;
53     scanf("%d", &T);
54     while(T--)
55     {
56         int n, m, p;
57         scanf("%d%d%d", &n, &m, &p);
58         LL ans = 0;
59         if(p >= N)
60         {
61             ans = (LL)n * m;
62         }
63         else
64         {
65             if(n > m) std::swap(n, m);
66             for(int i = 1, j; i <= n; i = j + 1)
67             {
68                 j = std::min(n/(n/i), m/(m/i));
69                 ans += ((LL)F[j][p] - F[i-1][p]) * (n/i) * (m/i);
70             }
71         }
72 
73         printf("%I64d\n", ans);
74     }
75 
76     return 0;
77 }
代码君

 

HDU 4746 (莫比乌斯反演) Mophues

标签:

原文地址:http://www.cnblogs.com/AOQNRMGYXLMV/p/4175284.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!