去年省选不会做,过了一年又要省选了还是不会做,滚粗既视感。
一维的时候我们知道,
∑ni=1d(i)=∑ni=1⌊ni⌋
,然后这个可以拓展到二维,然而我并不会证明
ans=∑i=1n∑j=1m⌊ni⌋⌊mj⌋,gcd(i,j)=1
然后莫比乌斯反演。
ans=∑i=1n∑j=1m⌊ni⌋⌊mj⌋∑d∣gcd(i,j)μ(d)
=∑d=1nμ(d)∑i=1⌊nd⌋∑j=1⌊md⌋⌊nid⌋⌊mjd⌋
=∑d=1nμ(d)∑i=1⌊nd⌋⌊nid⌋∑i=1⌊md⌋⌊njd⌋
=∑d=1nμ(d)f(⌊nd⌋)f(⌊md⌋)
这里的 f(x)=∑xi=1⌊xi⌋ ,可以发现就是 d(x) 的前缀和。
我们可以线性筛除 d(x) 然后就有了 f(x) 然后就分块就好辣。
#include<iostream>
#include<cstdio>
#define MAXN 50005
#define ll long long
using namespace std;
int mobius[MAXN],prime[MAXN],d[MAXN],f[MAXN];
bool flag[MAXN];
inline ll read()
{
ll a=0,f=1; char c=getchar();
while (c<'0'||c>'9') {if (c=='-') f=-1; c=getchar();}
while (c>='0'&&c<='9') {a=a*10+c-'0'; c=getchar();}
return a*f;
}
inline void prepare()
{
mobius[1]=d[1]=1;
for (int i=2;i<MAXN;i++)
{
if (!flag[i]) prime[++prime[0]]=i,mobius[i]=-1,d[i]=2,f[i]=1;
for (int j=1;j<=prime[0]&&i*prime[j]<MAXN;j++)
{
flag[i*prime[j]]=1;
if (i%prime[j]==0)
{
mobius[i*prime[j]]=0;
f[i*prime[j]]=f[i]+1;
d[i*prime[j]]=d[i]/(f[i]+1)*(f[i]+2);
break;
}
mobius[i*prime[j]]=-mobius[i];
f[i*prime[j]]=1;
d[i*prime[j]]=d[i]*2;
}
}
for (int i=1;i<MAXN;i++) d[i]+=d[i-1],mobius[i]+=mobius[i-1];
}
int main()
{
prepare();
int testcase=read();
while (testcase--)
{
ll n=read(),m=read();
if (n>m) swap(n,m);
ll ans=0;
for (int i=1,pos;i<=n;i=pos+1)
{
pos=min(n/(n/i),m/(m/i));
ans+=(ll)(mobius[pos]-mobius[i-1])*d[n/i]*d[m/i];
}
printf("%lld\n",ans);
}
return 0;
}