RT，WA5分，人麻了

#include <bits/stdc++.h>
#define int long long
#define debug puts("IAKIOI")
#define endl '\n'
#define inf 1e18
using namespace std;
const int mod=20101009,inv2=10050505,inv3=13400673;
int pm[1000010],tot,mu[10000010],s[10000010];
bool vis[10000010];
void init()
{
    mu[1]=1;
    for ( int i = 2 ; i <= 10000000 ; i++ )
    {
        if(!vis[i])
        {
            mu[i]=-1;
            pm[++tot]=i;
        }
        for ( int j = 1 ; j <= tot ; j++ )
        {
            if(pm[j]*i>10000000)
            {
                break;
            }
            vis[pm[j]*i]=1;
            if(i%pm[j]==0)
            {
                break;
            }
            mu[i*pm[j]]=-mu[i];
        }
    }
    for ( int i = 1 ; i <= 10000000 ; i++ )
    {
        mu[i]+=mod;
        mu[i]%=mod;
    }
    for ( int i = 1 ; i <= 10000000 ; i++ )
    {
        s[i]=i*i%mod*mu[i];
        s[i]+=mod;
        s[i]%=mod;
    }
    for ( int i = 1 ; i <= 10000000 ; i++ )
    {
        s[i]+=s[i-1]+mod;
        // s[i]+=mod;
        s[i]%=mod;
    }
}
// int pfh(int n)
// {
//     return n*(n+1)/2*(2*n+1)/3; 
// }
int calc(int n,int m)
{
    return n*(n+1)%mod*inv2%mod*m%mod*(m+1)%mod*inv2%mod;
}
int sum(int n,int m)
{
    int l=1,ans=0;
    while(l<=n)
    {
        int r=n/(n/l);
        int x=s[r]-s[l-1]+mod;
        x%=mod;
        // x=x*(pfh(r)-pfh(l-1))%mod;
        ans=ans+x*calc(n/l,m/l)%mod;
        ans%=mod;
        l=r+1;
    }
    return ans;
}
signed main()
{
	init();
    int n,m;
    cin >> n >> m;
    if(n>m)
    {
        swap(n,m);
    }
    int l=1,ans=0;
    while(l<=n)
    {
        int r=n/(n/l);
        int x=(l+r)%mod*(r-l+1)%mod*inv2%mod;
        ans=ans+x*sum(n/l,m/l)%mod;
        ans%=mod;
        l=r+1;
    }
    cout << ans;
	return 0;
}