标签：
这道题看巨巨的题解看了好久，好久。。
本文转自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都是一样的。
另外我想再补充一下：
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 }
标签：
原文地址：http://www.cnblogs.com/AOQNRMGYXLMV/p/4175284.html
