这题美好体验就是卡常
Description
求
[sum^{n}_ {i=1} sum^m_{j=1} [gcd(i,j)in prime]
]
其中：(1leq n,m leq 10^7)，多组询问
Solution
对于这种与(gcd)相关的反演题，有一个好的套路
设(f(d)=[gcd(i,j)=d]),(F(n))为(gcd(i,j)=d)和(d)的倍数的个数，即：
[f(d)=[gcd(i,j)=d]
]
[F(n)=sum_{d|n} f(d)= lfloor frac{N}{n}
floor lfloor frac{M}{m}
floor
]
[f(n)=sum_{n|d} mu(lfloor frac{d}{n}
floor)F(d)
]
由着这个套路,我们开始化简这个式子
[sum_{p in prime} sum_{i=1}^n sum_{j=1}^m [gcd(i,j)==p]
]
将(f(p))带入：
[sum_{p in prime}f(p)
]
把(f(x))换成(F(x))
[sum_{p in prime}sum_{p|d} mu(lfloor frac{d}{p}
floor)F(d)
]
我们枚举(lfloor frac{d}{p} floor)
[sum_{p in prime} sum_{d=1}^{min(lfloor frac{n}{p}
floor,lfloor frac{m}{p}
floor} mu(d)F(dp)
]
再把(F(dp))换成最终式：
[sum_{p in prime} sum_{d=1}^{min(lfloor frac{n}{p}
floor,lfloor frac{m}{p}
floor)} mu(d)lfloor frac{N}{n}
floor lfloor frac{M}{m}
floor
]
令(T=dp)，则有：
[sum^{min(n,m)}_ {T=1} sum_{t|T,tin prime} mu(lfloor frac{T}{t}
floor)lfloor frac{n}{T}
floorlfloor frac{m}{T}
floor
]
[sum^{min(n,m)}_ {T=1} lfloor frac{n}{T}
floorlfloor frac{m}{T}
floor(sum_{t|T} mu(lfloor frac{T}{t}
floor))
]
推到这里，我们就都可以做了
(mu(space))可以线性筛，其他的可以整除分块
CODE
#include<bits/stdc++.h>
#define reg register
using namespace std;
namespace yspm{
inline int read()
{
int res=0,f=1; char k;
while(!isdigit(k=getchar())) if(k=='-') f=-1;
while(isdigit(k)) res=res*10+k-'0',k=getchar();
return res*f;
}
const int N=10000010;
bool vis[N]; int pri[N],mu[N],g[N],cnt;
#define ll long long
ll sum[N];
inline void prework()
{
mu[1]=1;
for(reg int i=2;i<N;++i)
{
if(!vis[i]){mu[i]=-1;pri[++cnt]=i;}
for(reg int j=1;j<=cnt&&pri[j]*i<N;++j)
{
vis[i*pri[j]]=1;
if(i%pri[j]==0) break;
else mu[pri[j]*i]-=mu[i];
}
}
for(reg int j=1;j<=cnt;++j)
{
for(reg int i=1;i*pri[j]<N;++i) g[i*pri[j]]+=mu[i];
}
for(reg int i=1;i<N;++i) sum[i]=sum[i-1]+(ll)g[i];
return ;
}
inline void work()
{
int n=read(),m=read(); if(n>m) swap(n,m);
ll ans=0;
for(reg int l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
}printf("%lld
",ans);
return ;
}
signed main()
{
prework(); int T=read(); while(T--) work();
return 0;
}
}
signed main(){return yspm::main();}