HDU 4746 Mophues（莫比乌斯反演）

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

思路：不难看出可以看出p最多不超过19.

首先考虑这个问题的简化版，求gcd(a,b)小于等于d的数对的数量。

对于这个问题，由莫比乌斯反演我们可以求出gcd为i的数对数量然后从1到j累加起来，但是d太大的话时间复杂度无法承受。

但是我们注意到，对于每个gcd = i，设A(d)：gcd(a, b)=d的有多少种，设B(j): gcd(a, b)是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..)。

而每项的系数和要求的d无关，我们可以用类似筛法nlogn的时间预处理出来，那么这个问题就可以解决了。

考虑本题，只是多了一个p这个参数，那么我们只需要给系数再加一维，表示素因子个数小于等于p时的系数值，这个可以通过一个累加求出，这样我们就预处理出了所要用到的系数，但是这样还是会t,在计算答案时还要加上一个分段优化，即把(n/i)*(m/i)值相同的合并起来并对系数预处理出一个前缀和，也就是将[i, n/(n/i), m/(m/i)]合并起来。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#include<vector>
#include<map>
#include<queue>
#include<stack>
#include<string>
#include<map>
#include<set>
#include<ctime>
#define eps 1e-6
#define LL long long
#define pii pair<int, int>
//#pragma comment(linker, "/STACK:1024000000,1024000000")
using namespace std;
const int MAXN = 500000+50;
bool check[MAXN+10];
int prime[MAXN+10];
int mu[MAXN+10], pri_fac[MAXN+10], F[MAXN][20];
int Sum[MAXN][30];
void Moblus()
{
    mu[1] = 1;
    int tot = 0;
    for(int i = 2; i <= MAXN; i++)
    {
        if( !check[i] )
        {
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < tot; j++)
        {
            if(i * prime[j] > MAXN) break;
            check[i * prime[j]] = true;
            if( i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else
            {
                mu[i * prime[j]] = -mu[i];
            }
        }
    }
    for(int i = 0; i < tot; i++) {
		for(int j = prime[i]; j < MAXN; j += prime[i]) {
			int t = j;
			while(t && t%prime[i]==0) pri_fac[j]++, t /= prime[i];
		}
    }
} 
void init() {
	for(int i = 1; i < MAXN; i++) {
		for(int j = i; j < MAXN; j += i) {
			F[j][pri_fac[i]] += mu[j/i]; 
		} 
	}
	for(int i = 1; i < MAXN; i++) {
		for(int j = 1; j < 20; j++) F[i][j] += F[i][j-1];
	}
	for(int i = 1; i < MAXN; i++) {
		for(int j = 0; j < 20; j++) {
			Sum[i][j] = Sum[i-1][j] + F[i][j]; 
		}
	}
}
int main() {
    //freopen("input.txt", "r", stdin);
	Moblus(); 
	init();
	int q, n, m, p; 
	cin >> q;
	while(q--) {
		scanf("%d%d%d", &n, &m, &p);
		LL ans = 0;
		int t = min(n, m), last;
		p = min(19, p);
		for(int i = 1; i <= t; i = last+1) {
			last = min(n/(n/i), m/(m/i));
			ans += (LL)(Sum[last][p]-Sum[i-1][p])*(n/i)*(m/i);
		}
		cout << ans << endl;	
	}
    return 0;
}