题意: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;
}