解题思路
[ans=sumlimits_{i=1}^nsumlimits_{j=1}^mlcm(i,j)
]
[ans=sumlimits_{i=1}^nsumlimits_{j=1}^mfrac{i*j}{gcd(i,j)}
]
[ans=sumlimits_{d=1}^ndsumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}}i*j*[gcd(i,j)=1]
]
[ans=sumlimits_{d=1}^ndsumlimits_{d'=1}^{frac{n}{d}}mu(d')*d'^2sumlimits_{i=1}^{frac{n}{dd'}}sumlimits_{j=1}^{frac{n}{dd'}}i*j
]
发现后面是一个等差数列求和，记为(sum)
[ans=sumlimits_{T=1}^nsum(n/T,m/T)sumlimits_{d=1}^nmu(d)*d^2*frac{T}{d}
]
前面的可以数论分块，主要是求后面的。注意(d)和(frac{T}{d})不能直接约分，因为后者为向下取整。发现后面的东西其实是一个积性函数，可以晒出来的。设后面的为(f)，首先(f(p)=p-p*p)，然后考虑(f(i*p))，如果(i)中有(p)这个因子，那么(p)对(i)没有新的约数的贡献，因为新的约数中一定含有(p^2)这个因子，而对应的(mu)为(0)，(p)对(i)的贡献为(i)中(frac{i}{d})的(i)扩大，所以这种情况下(f(i*p)=f(i)*p)。否则(f(i*p)=f(i)*f(p))。
代码
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
using namespace std;
typedef long long LL;
const int N=1e7+5;
const int MOD=1e8+9;
inline int rd(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)) f=ch=='-'?0:1,ch=getchar();
while(isdigit(ch)) x=(x<<1)+(x<<3)+ch-'0',ch=getchar();
return f?x:-x;
}
inline int min(int x,int y){return x<y?x:y;}
int T,n,m,prime[N/10],cnt,f[N];
bool vis[N];
inline void prework(){
f[1]=1;int lim=1e7;
for(int i=2;i<=lim;i++){
if(!vis[i]) {prime[++cnt]=i; f[i]=(i-(LL)i*i%MOD+MOD)%MOD;}
for(int j=1;j<=cnt && (LL)prime[j]*i<=lim;j++){
vis[i*prime[j]]=1;
if(!(i%prime[j])) {
f[i*prime[j]]=(LL)f[i]*prime[j]%MOD;
break;
}
f[i*prime[j]]=(LL)f[i]*f[prime[j]]%MOD;
}
}
for(int i=1;i<=lim;i++) f[i]=(f[i-1]+f[i])%MOD;
}
inline int calc(int x,int y){
return (LL)((LL)(x+1)*x/2%MOD)*((LL)(y+1)*y/2%MOD)%MOD;
}
inline void work(){
int ans=0;
n=rd(),m=rd();if(n>m) swap(n,m);
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans=(ans+(LL)calc(n/l,m/l)*(f[r]-f[l-1]+MOD)%MOD)%MOD;
}
printf("%d
",ans);
}
int main(){
prework();T=rd();
while(T--) work();
return 0;
}