慢慢化柿子吧
要求的是这个
[sum_{i=1}^Nsum_{j=1}^Md(ij)
]
神奇的约数个数函数有一个这样的性质
[d(ij)=sum_{x|i}sum_{y|j}[(x,y)=1]
]
试着从唯一分解定理的角度去理解，将(i,j)分别分解质因数
显然(d(ij))应该等于每一个(p)在(i,j)中分解出来的指数加起来加1再相乘
所以分别枚举所有约数的话，保证这对约数互质就可以求出所有约数了
之后现在的答案变成了
[sum_{i=1}^Nsum_{j=1}^Msum_{x|i}sum_{y|j}[(x,y)=1]
]
一看就很烦人啊，把(x,y)放到前面来枚举
[sum_{x=1}^Nsum_{y=1}^M[(x,y)=1]sum_{x|i}sum_{y|j}1
]
显然后面那个东西就是
[left lfloor frac{N}{x}
ight
floor imes left lfloor frac{M}{y}
ight
floor
]
于是答案变成了
[sum_{x=1}^Nsum_{y=1}^M[(x,y)=1]left lfloor frac{N}{x}
ight
floor imes left lfloor frac{M}{y}
ight
floor
]
开始套路了
[f(n)=sum_{i=1}^Nsum_{j=1}^M[(i,j)=n]left lfloor frac{N}{i}
ight
floor imes left lfloor frac{M}{j}
ight
floor
]
[F(n)=sum_{i=1}^Nsum_{j=1}^M[n|(i,j)]left lfloor frac{N}{i}
ight
floor imes left lfloor frac{M}{j}
ight
floor
]
[=sum_{i=1}^{ left lfloor frac{N}{i}
ight
floor}sum_{j=1}^{ left lfloor frac{M}{j}
ight
floor} left lfloor frac{N}{in}
ight
floor left lfloor frac{M}{jn}
ight
floor
]
显然就有
[F(n)=sum_{n|d}f(d)
]
[f(n)=sum_{n|d}mu(frac{d}{n})F(d)
]
因为我们要求的只有(f(1))
所以
[Ans=sum_{d=1}^Nmu(d)sum_{i=1}^{left lfloor frac{N}{i}
ight
floor}sum_{j=1}^{left lfloor frac{M}{j}
ight
floor}left lfloor frac{N}{id}
ight
floor left lfloor frac{M}{jd}
ight
floor
]
我们发现如果(frac{N}{d})和(frac{M}{d})固定了，后面这个柿子就非常好求了
就是(frac{N}{d})和(frac{M}{d})的约数个数前缀和
于是我们可以线筛约数个数函数和(mu)之后就可以了
代码
#include<iostream>
#include<cstdio>
#include<cstring>
#define LL long long
#define min(a,b) ((a)<(b)?(a):(b))
#define max(a,b) ((a)>(b)?(a):(b))
#define re register
#define maxn 50005
inline int read()
{
char c=getchar();
int x=0;
while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();
return x;
}
int p[maxn>>1],f[maxn],pre[maxn],mu[maxn],d[maxn],t[maxn];
int T,n,m;
int main()
{
f[1]=1,mu[1]=1,d[1]=1;
for(re int i=2;i<=50000;i++)
{
if(!f[i]) p[++p[0]]=i,mu[i]=-1,d[i]=2,t[i]=1;
for(re int j=1;j<=p[0]&&p[j]*i<=50000;j++)
{
f[p[j]*i]=1;
if(i%p[j]==0)
{
d[p[j]*i]=d[i]/(t[i]+1)*(t[i]+2);
t[p[j]*i]=t[i]+1;
break;
}
d[p[j]*i]=d[p[j]]*d[i];
mu[p[j]*i]=-1*mu[i];
t[p[j]*i]=1;
}
}
for(re int i=1;i<=50000;i++) pre[i]=pre[i-1]+mu[i],d[i]+=d[i-1];
T=read();
while(T--)
{
n=read(),m=read();
if(n>m) std::swap(n,m);
LL ans=0;
for(re int l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans+=(LL)d[n/l]*d[m/l]*(pre[r]-pre[l-1]);
}
printf("%lld
",ans);
}
return 0;
}