题解都在论文里了

#include
     
      using namespace std;#define maxn 10000005#define ll long long #define mod 20101009bool vis[maxn];int sum[maxn],prime[maxn],mm,mu[maxn];void primes(){    mu[1]=1;    for(int i=2;i
      
       =maxn)break;            vis[i*prime[j]]=1;            if(i%prime[j]==0){                mu[i*prime[j]]=0;                break;            }            else mu[i*prime[j]]=-mu[i];        }    }     for(int i=1;i
       
        m)swap(n,m);    for(int l=1,r;l<=n;l=r+1){        r=min(n/(n/l),m/(m/l));        ll tmp=((sum[r]-sum[l-1])%mod+mod)%mod;         res=(res+tmp*Sum(n/l,m/l)%mod)%mod;    }    return res;}int main(){    primes();        scanf("%lld%lld",&n,&m);    if(n>m)swap(n,m);    ll ans=0;    for(int l=1,r;l<=n;l=r+1){        r=min(n/(n/l),m/(m/l));        ll tmp=(ll)(l+r)*(r-l+1)/2%mod;        ans=(ans+tmp*F(n/l,m/l)%mod)%mod;    }    cout<
        
         <