HDU 4746 Mophues 莫比乌斯第三弹

题意：1<=x<=n,1<=y<=m,使得gcd(x,y)=k，k的素因数个数小于等于p

例：24=2*2*2*3，k=4

解：设f[n]为gcd(a,b)=n的对数

      F[d]为d|gcd(a,b)的对数

    f[n]=sigema(mu[i],F[i*n])：

    f[1]=mu[1]*F[1]+mu[2]*F[1*2]+...+mu[n]*F[1*n]

    f[2]=mu[2]*F[2]+mu[2]*F[2*2]+...+mu[n]*F[2*n]

    ......

    sum=f[1]+f[2]+...+f[n]=G[1]*F[1]+G[2]*F[2]+...+G[n]*F[n]

    枚举每一个i，则i的倍数j为G[j]提供了mu[j/i]的贡献，即G[j]+=mu[j/i]

    因为所求为k的素因数个数，所以将G[j]开成二维数组G[j][p]表示j对素因数个数为p的贡献

    需要使用分块加速的方法，否则还是要爆炸

#include <stdio.h>
#include <string.h>
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)<(b)?(b):(a))
#define ll __int64
const int maxn=500005;
int num[maxn];
int prime[maxn];
int mu[maxn];
int factor[maxn];
int mbs[maxn][20];
void mobius()
{
    memset(num,0,sizeof(num));
    int all=0;
    mu[1]=1;
    factor[1]=0;
    for(int i=2;i<maxn;i++)
    {
        if(!num[i])
        {
            prime[all++]=i;
            mu[i]=-1;
            factor[i]=1;                 //记录素因数个数
        }
        for(int j=0;j<all&&i*prime[j]<maxn;j++)
        {
            num[i*prime[j]]=1;
            factor[i*prime[j]]=factor[i]+1;
            if(i%prime[j])
            {
                mu[i*prime[j]]=-mu[i];
            }
            else
            {
                mu[i*prime[j]]=0;
                break;
            }
        }
    }
    return ;
}
void inti()
{
    memset(mbs,0,sizeof(mbs));
    for(int i=1;i<maxn;i++)                  
        for(int j=i;j<maxn;j+=i)
            mbs[j][factor[i]]+=mu[j/i];        //每个j在factor[i]个素因数中的贡献

/*下面是为了分块加速求和*/
    for(int i=1;i<maxn;i++)
        for(int j=0;j<19;j++)
            mbs[i][j]+=mbs[i-1][j];

    for(int i=0;i<maxn;i++)
        for(int j=1;j<19;j++)
            mbs[i][j]+=mbs[i][j-1];
    return ;
}
int main()
{
    int t,n,m,p;
    mobius();
    inti();
    ll sum;
    while(scanf("%d",&t)!=-1)
    {
        while(t--)
        {
            scanf("%d%d%d",&n,&m,&p);
            if(p>=19)                    //因为2^19>500000，所以超过了19就是全体均满足要求
            {
                printf("%I64d\n",(ll)n*m);
                continue;
            }
            if(n>m)
            {
                int te=n;
                n=m;
                m=te;
            }
            sum=0;
            for(int i=1,last;i<n;i=last+1)
            {
                last=MIN(n/(n/i),m/(m/i));           

//分块加速，因为[n/i][m/i]在i递加过程中具有重复部分，跳掉这些i可是简化计算，是复杂度降低至sqrt(n)
                sum+=((ll)(n/i)*(m/i)*(mbs[last][p]-mbs[i-1][p]));
            }
            printf("%I64d\n",sum);
        }
    }
    return 0;
}