hdoj 4746 莫比乌斯反演 + 优化

hdoj 4746

题意：给三个整数n、m、p，问有多少数对满足gcd（n，m）的素因子（可重复）的个数小于等于p。

思路：

设函数A(x)为在n和m的所有数对中，gcd = x的数对个数；函数B(x)为在n和m的所有数对中，公约数为x的数对个数，显然B(x) = (n / x) * (m / x)。

显然A的数对包含于B的数对中。

根据容斥原理得A(x) = u(1) * B(x * 1) + u(2) * B(x * 2) + u(3) * B(x * 3)...

ans = sigma(A(i))(1<= i < min(n, m), i的素因子个数小于等于p)

根据规律可以发现u(1)、u(2)、u(3)...满足莫比乌斯函数。

到这里这个题目的理论部分基本就完成了，但是复杂度很高，而且这个题目时间卡的比较紧，重点在优化方面。

看这样一组样例：

A(1) = u(1) * B(1) + u(2) * B(2) + u(3) * B(3) + u(4) * B(4) + u(5) * B(5)...

A(2) = u(1) * B(2) + u(2) * B(4) + u(3) * B(6) + u(4) * B(8) + u(5) * B(10)...

A(3) = u(1) * B(3) + u(2) * B(6) + u(3) * B(9) + u(4) * B(12) + u(5) * B(15)...

我们发现如果计算m=n=3，p=1，那么结果就是

A(1) + A(2) + A(3)= u(1) * B(1) + (u(1) + u (2)) * B(2) + ...

其中A(1)是素因子个数为0的，A(2)、A(3)是个数为1的。

这样我们就可以把u(i)以前缀和的形式整理出来，并以素因子个数分组，即sum[i][j]是数0到i，素因子个数小于等于j的函数U(x)的和，然后遍历1到min(n, m)就可以求出答案。

这样我们的时间复杂度就降到了O(n)。

然而这样依然过不了这题！！

虽然预处理时间稍微有点长，但是这道题的test数量有点多，时间真的非常紧。

然后就需要进行分块。

我们可以发现，枚举的数可以分成若干段，这一段中的n/i和m/i都是不变的。

分段规则是[i, min(n / (n / i), m / (m / i))]。

通过这样的加速运算。。我的代码才压线过。。我真是太水了。。

#include <cstring>
#include <cstdio>
#include <algorithm>
using namespace std;
const int M = 500005;
int mu[M], prmdiv[M], presum[M][20];
bool isPrime[M];
int cal(int a, int b){
int ret = 0;
while(b % a == 0) {
b /= a;
ret++;
}
return ret;
}
void init(){
memset(mu, -1, sizeof mu);
memset(isPrime, 0, sizeof isPrime);
memset(prmdiv, 0, sizeof prmdiv);
memset(presum, 0, sizeof presum);
for(int i = 2; i < M; i++) {//筛素数 求素因子个数 并将有重复素因子个数的数字的莫比乌斯函数设为0
if(isPrime[i]) continue;
isPrime[i] = false, prmdiv[i] = 1;
for(long long j = i + i; j < M; j += i){
isPrime[j] = true;
int t = cal(i, j);
prmdiv[j] += t;
if(t > 1) mu[j] = 0;
}
}
for(int i = 1; i < M; i++)
if(mu[i] != 0 && prmdiv[i] % 2 == 0) mu[i] = 1;
for(int i = 1; i < M; i++){//莫比乌斯函数前缀和
if(mu[i] == 0) continue;
for(long long j = i; j < M; j += i) {
presum[j][prmdiv[j / i]] += mu[i];
}
}
for(int i = 1; i < M; i++)
for(int j = 1; j <= 18; j++){
presum[i][j] += presum[i][j - 1];
}
for(int j = 0; j <= 18; j++)
for(int i = 2; i <= M; i++)
presum[i][j] += presum[i - 1][j];
}
main() {
init();
long long n, m, t, p;
scanf("%I64d", &t);
while(t--) {
scanf("%I64d %I64d %I64d", &n, &m, &p);
if(p >= 18) {
printf("%I64d\n", n * m);
continue;
}
if(n > m) swap(n, m);
long long ans = 0;
for(int i = 1, j; i <= n; i = j + 1){
j = min(n / (n / i), m / (m / i));//分块加速
ans += (long long)(presum[j][p] - presum[i - 1][p])* (n / i) * (m / i);
}
printf("%I64d\n", ans);
}
}