【莫比乌斯】 HDU 4746 Mophues

通道

题意：给出n, m, p，求有多少对a, b满足gcd(a, b)的素因子个数<=p(1<=a<=n, 1<=b<=m)

思路：转自http://blog.csdn.net/acdreamers/article/details/12871643

设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 <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

typedef long long ll;

template <class T>
inline bool rd(T &ret) {
    char c; int sgn;
    if(c = getchar() , c == EOF) return false;
    while(c != '-' && (c < '0' || c > '9')) c = getchar();
    sgn = (c == '-') ? -1 : 1;
    ret = (c == '-') ? 0 : (c - '0');
    while(c = getchar(), c >= '0' && c <= '9') ret = ret * 10 + (c - '0');
    ret *= sgn;
    return true;
}

const int MAX_N = 500007;
const int MAX_M = 21;

bool vis[MAX_N];
int mu[MAX_N], cnt[MAX_N];
int dp[MAX_N][MAX_M];

void Mobius(int n) {
    vis[1] = mu[1] = 1;
    for(int i = 1;i <= n; i++)  mu[i] = 1;
    for(int i = 2;i <= n; i++) {
        if(!vis[i]) {
            for(int j = i;j <= n;j += i) {
                vis[j] = 1;
                mu[j] *= -1;
                if((j/i) % i == 0)  mu[j] = 0;
                int t = j;
                while(t % i == 0) t /= i, cnt[j]++;
            }
        }
    }
}

int n, m, p;

int main() {
    Mobius(MAX_N - 2);
    for (int i = 1; i < MAX_N; ++i) 
        for (int j = i; j < MAX_N; j += i) 
            dp[j][cnt[i]] += mu[j / i];
    for (int i = 1; i < MAX_N; ++i) 
        for (int j = 1; j < MAX_M; ++j)
            dp[i][j] += dp[i][j - 1];
    for (int i = 1; i < MAX_N; ++i) 
        for (int j = 0; j < MAX_M; ++j)
            dp[i][j] += dp[i - 1][j];
    int T; rd(T);
    while (T-- > 0) {
        rd(n), rd(m), rd(p);
        ll ans = 0;
        if (p >= MAX_M) ans = (ll) n * m;
        else {
            if (n > m) swap(n, m);
            for (int i = 1, j = 0; i <= n; i = j + 1) {
                j = min(n / (n / i), m / (m / i));
                ans += ((ll)dp[j][p] - dp[i - 1][p]) * (n / i) * (m / i);
            }
        }
        printf("%I64d\n", ans);
    }
    return 0;
}