hdu4746 Mophues

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4746

题意：给出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)的贡献系数

代码：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>

using namespace std;
const int maxn=500010;

int n,m,p;
int vis[maxn];
int prime[maxn];
int cnt;
int num[maxn];
int mu[maxn];
int f[maxn][20];

void init()
{
    memset(vis,0,sizeof(vis));
    cnt=0;
    mu[1]=1;
    for(int i=2;i<maxn;i++)
    {
        if(!vis[i])
        {
            prime[cnt++]=i;
            mu[i]=-1;
            num[i]=1;
        }
        for(int j=0;j<cnt&&i*prime[j]<maxn;j++)
        {
            vis[i*prime[j]]=1;
            num[i*prime[j]]=num[i]+1;
            if(i%prime[j])
                mu[i*prime[j]]=-mu[i];
            else
            {
                mu[i*prime[j]]=0;
                break;
            }
        }
    }
    memset(f,0,sizeof(f));
    for(int i=1;i<maxn;i++)
        for(int j=i;j<maxn;j+=i)
            f[j][num[i]]+=mu[j/i];

    for(int i=0;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++)
            f[i][j]+=f[i-1][j];
}

int main()
{
    init();
    int T;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d",&n,&m,&p);
        if(p>=20)
        {
            printf("%lld
",(long long)n*m);
            continue;
        }
        if(n>m)
            swap(n,m);
        int last;
        long long ans=0;
        for(int i=1;i<=n;i=last+1)
        {
            last=min(n/(n/i),m/(m/i));
            ans+=(long long)(f[last][p]-f[i-1][p])*(n/i)*(m/i);
        }
        printf("%lld
",ans);
    }
    return 0;
}