HDU4746(莫比乌斯反演，分块加速）

题意：给出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<algorithm>
#include<cstdio>
#include<cstring>
#include<iostream>
#define N 500009
#define M 19
#define ll long long
using namespace std;
//F[i][j]: 表示素因子个数<=j条件下的莫比乌斯前缀和：μ(1)+μ(2)+...+μ(i)  
int h[N],num[N],F[N][M];//h用于处理mobpus函数，num用于统计质因数
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 cal(int n, int x) 
{  
    int res = 0;  
    do {  
        ++res;  
        n /= x;  
    } while (n % x == 0);  
    return res;  
}  
void init()
{
	for(int i=2;i<N;i++)
	{
		if(num[i]) continue;//说明不是质数 
		int cnt;
		for(int j=i;j<N;j+=i)
		{
			cnt=cal(j, i);
			num[j]+=cnt;
			if(cnt>1)
			{
				h[j]=-1;
			}else 
			if(h[j]>=0)
			{
				h[j]++; 
			}
		}
	}
	for(int i=1;i<N;i++)
	{
		for(int j=i;j<N;j+=i)
		{
			F[j][num[i]]+=mob(j/i);
		}
	}
    for (int i = 1; i < N; i++) 
	{  
        for (int j = 1; j < M; j++) 
		{  
            F[i][j] += F[i][j-1];  
        }  
    }
    //为了分组加速求解，求i的前缀和  
    for (int i = 1; i < N; i++) 
	{  
        for (int j = 0; j < M; j++) 
		{  
            F[i][j] += F[i-1][j];  
        }  
    }    
}
int main()
{
	//freopen("t.txt","r",stdin);
	int T,n,m,p;
	init();
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d%d%d",&n,&m,&p);
		ll ans=0;
		if(p>=M)
		{
			ans =(ll)n*m;  
		}else
		{
			if(n>m) 
			{
				swap(n,m);
			}
			int j;
			for(int 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;
 }