題目：http://acm.hdu.edu.cn/showproblem.php?pid=4746

?

題意：給出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都是一樣的。

?

#include <iostream> #include <cstdio> #include <cmath> #include <algorithm> using namespace std; #define LL long long #define M 500005 #define N 19//返回n中有多少個x因子 int cal(int n, int x) {int res = 0;do{++res;n /= x;}while (n % x == 0);return res; }//備注：分塊加速求解需要求前綴和 //F[i][j]: 表示素因子個數<=j條件下的莫比烏斯前綴和：μ(1)+μ(2)+...+μ(i) int F[M][N]; int num[M]; //num[i]: i中含有多少個素因子 int h[M]; //h[i]: -1表示存在平方因子，否則表示有多少種素因子//莫比烏斯函數的定義 int mob(int n) {if (h[n] == -1) return 0; //存在平方因子時，μ(n)=0if (h[n] & 1) return -1; //奇數個不同素數之積，μ(n)=-1return 1; //偶數個不同素數之積，μ(n)=1 }int main() {int t, n, m, p, i, j;//篩法算出num[]以及h[]for (i = 2; i < M; i++){if (num[i]) continue;for (j = i; j < M; j+=i){int tp = cal(j, i);num[j] += tp;if (tp > 1) //j中含有多個i，必然存在平方因子{h[j] = -1;}else if (h[j] >= 0){++h[j];}}}//枚舉i作為公因子，對B(j)的貢獻值為：mob(j/i)for (i = 1; i < M; i++){for (j = i; j < M; j+=i){F[j][num[i]] += mob(j/i);}}//為了表示素因子數<=j的意義，求j的前綴和for (i = 1; i < M; i++){for (j = 1; j < N; j++){F[i][j] += F[i][j-1];}}//為了分組加速求解，求i的前綴和for (i = 1; i < M; i++){for (j = 0; j < N; j++){F[i][j] += F[i-1][j];}}scanf("%d", &t);while (t--){scanf("%d%d%d", &n, &m, &p);LL ans = 0;if (p >= N){ans = (LL)n*m;}else{if (n > m){n ^= m;m ^= n;n ^= m;}for (i = 1; i <= n; i = j + 1){j = min(n/(n/i), m/(m/i));ans += ((LL)F[j][p]-F[i-1][p])*(n/i)*(m/i);}}printf("%I64d\n", ans);}return 0; }

?

總結

以上是生活随笔為你收集整理的hdu4746(莫比乌斯反演)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。