Description
求(sumlimits_{i=1}^nsumlimits_{j=1}^m lcm(i,j))，答案模1e9+9输出，多组询问
Input
一个正整数T表示数据组数
接下来T行 每行两个正整数 表示N、M
Output
T行 每行一个整数 表示第i组数据的结果
Sample Input
1
4 5
Sample Output
122
HINT
T <= 10000
N, M<=10000000
我们令n<m，然后将柿子化简
[sumlimits_{i=1}^nsumlimits_{j=1}^m lcm(i,j)
]
[sumlimits_{i=1}^nsumlimits_{j=1}^m dfrac{i imes j}{gcd(i,j)}
]
[sumlimits_{d=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^mdfrac{i imes j}{d}[gcd(i,j)=d]
]
[sumlimits_{d=1}^nsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor}dfrac{d^2 imes i imes j}{d}[gcd(i,j)=1]
]
[sumlimits_{d=1}^n dsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor}i imes jsumlimits_{x|i,x|j}mu(x)
]
[sumlimits_{d=1}^n dsumlimits_{x=1}^{lfloorfrac{n}{d}
floor}mu(x) imes x^2sumlimits_{i=1}^{lfloorfrac{n}{dx}
floor}sumlimits_{j=1}^{lfloorfrac{m}{dx}
floor}i imes j
]
我们发现最后那个是等差数列，继续化简
[sumlimits_{d=1}^n dsumlimits_{x=1}^{lfloorfrac{n}{d}
floor}mu(x) imes x^2dfrac{lfloorfrac{n}{dx}
floor(lfloorfrac{n}{dx}
floor+1)lfloorfrac{m}{dx}
floor(lfloorfrac{m}{dx}
floor+1)}{4}
]
然后我们令(T=dx)，那么得到
[sumlimits_{T=1}^ndfrac{lfloorfrac{n}{T}
floor(lfloorfrac{n}{T}
floor+1)lfloorfrac{m}{T}
floor(lfloorfrac{m}{T}
floor+1)}{4}Tsumlimits_{x|T}mu(x)x
]
我们设(f(T)=sumlimits_{x|T}mu(x)x)，预处理出f，就可以分块了
设(g(x)=mu(x)x)，当a,b互质，(g(a) imes g(b)=abmu(a)mu(b)=abmu(ab)=g(ab))，所以g是积性函数，根据莫比乌斯反演的性质，f也是积性函数
令(T=prodlimits_{i=1}^k P_i^{x_i})，(f(P_i^{x_i})=(1-P_i))，那么
[f(T)=prodlimits_{i=1}^k(1-P_i)
]
这样子我们可以在(O(n))时间内线筛出来，然后维护一下(T imes f(T))的前缀和，然后就可以在(O(sqrt N))的时间内完成每次询问
/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline int read(){
int x=0,f=1;char ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';
return x*f;
}
inline void print(int x){
if (x>=10) print(x/10);
putchar(x%10+'0');
}
const int N=1e7,p=1e8+9;
int prime[N+10],f[N+10],tot;
bool inprime[N+10];
void prepare(){
f[1]=1;
for (int i=2;i<=N;i++){
if (!inprime[i]) prime[++tot]=i,f[i]=1-i+p;
for (int j=1;j<=tot&&i*prime[j]<=N;j++){
inprime[i*prime[j]]=1;
if (i%prime[j]==0){
f[i*prime[j]]=f[i];
break;
}
f[i*prime[j]]=1ll*f[i]*f[prime[j]]%p;
}
}
for (int i=1;i<=N;i++) f[i]=(f[i-1]+1ll*i*f[i]%p)%p;
}
int get(int x){return (1ll*x*(x+1)>>1)%p;}
int main(){
prepare();
for (int Data=read();Data;Data--){
int n=read(),m=read(),pos,Ans=0;
if (n>m) swap(n,m);
for (int T=1;T<=n;T=pos+1){
pos=min(n/(n/T),m/(m/T));
Ans=(Ans+1ll*get(n/T)*get(m/T)%p*(f[pos]-f[T-1]+p)%p)%p;
}
printf("%d
",Ans);
}
return 0;
}