莫比乌斯函数完整定义的通俗表达:
1)莫比乌斯函数μ(n)的定义域是N
2)μ(1)=1
3)当n存在平方因子时,μ(n)=0
4)当n是素数或奇数个不同素数之积时,μ(n)=-1
5)当n是偶数个不同素数之积时,μ(n)=1
/* * [题意] * 给出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)=0 if (h[n] & 1) return -1; //奇数个不同素数之积,μ(n)=-1 return 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; }
相关推荐
HDU的1250,主要是利用高精度加法,但是代码有点繁琐,效率不是很高
HDU1059的代码
杭电ACMhdu1163
HDU的一题........HDU DP动态规
hdu1001解题报告
hdu 1574 passed sorce
hdu acm 教案 搜索入门 hdu acm 教案 搜索入门
hdu2101AC代码
ACM HDU题目分类,我自己总结的大概只有十来个吧
搜索 dfs 解题代码 hdu1241
hdu acm 教案 动态规划(1) hdu acm 教案 动态规划(1)
hdu 5007 Post Robot 字符串枚举。 暴力一下就可以了。
hdu-acm源代码(上百题)hdu-acm源代码、hdu-acm源代码hdu-acm源代码
HDU最全ac代码
hdu动态规划算法集锦
hdu 1166线段树代码
hdu 1005.比较简单的一道题,有兴趣的可以看看。
hdu题目分类
HDU图论题目分类
ACM HDU 2000->2099 解题报告 ACM HDU 2000->2099 解题报告 ACM HDU 2000->2099 解题报告