给(n,m),让你求
[sumlimits_{i=1}^n sumlimits_{j=1}^m [gcd(i,j) in prime]
]
有(T)组询问((T le 10^4,n,mle 10^7))。
枚举质数(p),然后柿子变成
[sumlimits_{p} sumlimits_{i=1}^nsumlimits_{j=1}^m [gcd(i,j)=p]
]
等价于
[sumlimits_{p} sumlimits_{i=1}^{leftlfloor n/p
ight
floor}sumlimits_{j=1}^{leftlfloor m/p
ight
floor} [gcd(i,j)=1]
]
因为把(leftlfloor n/p ight floor,leftlfloor m/p ight floor)以内互质的一对数乘上(p)就是(gcd=p)的数了。
令(S(n,m)=sumlimits_{i=1}^nsumlimits_{j=1}^m [gcd(i,j)=1]),由于(mu)的性质([n=1]=sumlimits_{dmid n} mu(d)),所以
[egin{aligned}S(n,m)&=sumlimits_{i=1}^n sumlimits_{j=1}^m sumlimits_{dmid gcd} mu(d) \&= sumlimits_{d}mu(d)sumlimits_{i=1}^n [dmid i]sumlimits_{j=1}^m [dmid j] \&= sumlimits_{d} mu(d)leftlfloor frac{n}{d}
ight
floorleftlfloor frac{m}{d}
ight
floorend{aligned}
]
又因为(leftlfloorfrac{leftlfloorfrac{n}{p} ight floor}{d} ight floor=leftlfloorfrac{n}{dp} ight floor),所以柿子是
[sumlimits_{p} sumlimits_{d} mu(d) leftlfloorfrac{n}{dp}
ight
floor leftlfloorfrac{m}{dp}
ight
floor
]
可以枚举(k=dp),有
[sumlimits_{k} leftlfloorfrac{n}{k}
ight
floor leftlfloorfrac{m}{k}
ight
floorsumlimits_{pmid k} muleft(frac{k}{p}
ight)
]
令(f(n)=sumlimits_{pmid n} mu(frac{n}{p})),预处理出(f)的前缀和,数论分块就好了。
并不会算预处理的复杂度qwq...
#include <bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for (int i=(a);i<(b);++i)
#define per(i,a,b) for (int i=(a)-1;i>=(b);--i)
#define mp make_pair
#define pb push_back
#define fi first
#define se second
#define all(x) (x).begin(),(x).end()
#define SZ(x) ((int)(x).size())
typedef double db;
typedef long long ll;
typedef pair<int,int> PII;
typedef vector<int> VI;
const int maxn=1e7,N=maxn+10;
int vis[N],p[N],pn,mu[N],sum[N];
#define ss(l,r) (sum[r]-sum[l-1])
void init(int n) {
mu[1]=1;
rep(i,2,n+1) {
if(!vis[i]) {p[pn++]=i;mu[i]=-1;}
for(int j=0;j<pn&&i*p[j]<=n;j++) {
vis[i*p[j]]=1;
if(i%p[j]==0) {mu[i*p[j]]=0;break;}
else mu[i*p[j]]=-mu[i];
}
}
rep(i,0,pn) for(int j=p[i];j<=n;j+=p[i])
sum[j]+=mu[j/p[i]];
rep(i,1,n+1) sum[i]+=sum[i-1];
}
ll solve(int n,int m) {
int tn=min(n,m); ll ans=0;
for(int l=1,r=0;l<=tn;l=r+1) {
r=min(n/(n/l),m/(m/l));
ans+=1ll*(n/l)*(m/l)*ss(l,r);
}
return ans;
}
int main() {
#ifdef LOCAL
freopen("a.in","r",stdin);
#endif
init(maxn);
int _,n,m;for(scanf("%d",&_);_;_--) {
scanf("%d%d",&n,&m);
printf("%lld
",solve(n,m));
}
return 0;
}