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);    }}