( ext{Solution:})
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int MAXN=1e7+10;
const int mod=20101009;
bitset<MAXN+10>vis;
int p[MAXN+10],mu[MAXN+10],cnt,n,m,inv2,inv6;
void predo(int g){
mu[1]=1;
for(int i=2;i<=g;++i){
if(!vis[i])mu[i]=-1,p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=g;++j){
vis[i*p[j]]=1;
if(i%p[j]==0)break;
mu[i*p[j]]=-mu[i];
}
}
for(int i=1;i<=g;++i)mu[i]=mu[i-1]+mod+(i*i*mu[i]),mu[i]=(mu[i]+mod)%mod;
}
inline int mul(int a,int b){return 1ll*a*b%mod;}
inline int add(int a,int b){return (a+b+mod)%mod;}
inline int qpow(int a,int b){
int res=1;
while(b){
if(b&1)res=mul(res,a);
a=mul(a,a);b>>=1;
}
return res;
}
inline int G(int a,int b){return mul(mul(mul(a,a+1),inv2),mul(mul(b,b+1),inv2));}
inline int calc(int N,int M){
int res=0;
for(int l=1,r;l<=min(N,M);l=r+1){
r=min(N/(N/l),M/(M/l));
res=add(res,mul(mu[r]-mu[l-1]+mod,G(N/l,M/l)));
}
return res;
}
inline int solve(int N,int M){
int res=0;
for(int l=1,r;l<=min(N,M);l=r+1){
r=min(N/(N/l),M/(M/l));
res=add(res,mul(calc(N/l,M/l),mul(r-l+1,mul(l+r,inv2))));
}
return res;
}
signed main(){
scanf("%lld%lld",&n,&m);
predo(min(n,m));inv2=qpow(2ll,mod-2);
printf("%lld
",solve(n,m));
return 0;
}
注意(k^2*mu(i))是要一起处理的,不是两个前缀相乘。以及,记得处理(2)的逆元。模数是质数。