莫比乌斯反演推柿子，数论分块降复杂度，最后时间复杂度为O(n)。
1 #include<bits/stdc++.h> 2 #define ll long long 3 using namespace std; 4 const int N=1e7+10,mod=20101009; 5 int n,m,mi; 6 ll p[N],tot; 7 ll u[N]; 8 bool isp[N]; 9 void getu(int lim){ 10 u[1]=1; 11 for(int i=2;i<=lim;i++){ 12 if(!isp[i]) u[i]=-1,p[++tot]=i; 13 for(int j=1;j<=tot&&1LL*i*p[j]<=lim;j++){ 14 isp[p[j]*i]=1,u[i*p[j]]=-u[i]; 15 if(!(i%p[j])){u[i*p[j]]=0;break;} 16 } 17 } 18 for(int i=2;i<=lim;i++) u[i]=((u[i]*i*i+u[i-1])%mod+mod)%mod; 19 return; 20 } 21 ll solve(int d){ 22 ll x=n/d,y=m/d,l=1,r,lim=mi/d,res=0; 23 while(l<=lim){ 24 r=min(x/(x/l),y/(y/l)); 25 ll temp=((x/l+1)*(x/l)/2%mod)*((y/l+1)*(y/l)/2%mod)%mod; 26 res=((res+(u[r]-u[l-1])*temp)%mod+mod)%mod; 27 l=r+1; 28 } 29 return res; 30 } 31 int main(){ 32 scanf("%d%d",&n,&m); 33 mi=min(n,m); 34 getu(mi); 35 ll ans=0; 36 int l=1,r; 37 while(l<=mi){ 38 r=min(n/(n/l),m/(m/l)); 39 ans=(ans+solve(l)*(l+r)*(r-l+1)/2+mod)%mod; 40 l=r+1; 41 } 42 printf("%lld",ans); 43 return 0; 44 }