思路：

莫比乌斯反演+整除分块

代码：

#pragma GCC optimize(2)#pragma GCC optimize(3)#pragma GCC optimize(4)#include
    
     using namespace std;#define y1 y11#define fi first#define se second#define pi acos(-1.0)#define LL long long//#define mp make_pair#define pb push_back#define ls rt<<1, l, m#define rs rt<<1|1, m+1, r#define ULL unsigned LL#define pll pair
     
      #define pli pair
      
       #define pii pair
       
        #define piii pair
        
         #define pdd pair
         
          #define mem(a, b) memset(a, b, sizeof(a))#define debug(x) cerr << #x << " = " << x << "\n";#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);//head const int MOD = 20101009;const int N = 1e7 + 10;int prime[N/10], mu[N], cnt;LL c[N], sum[N];bool not_p[N];int n, m;inline void seive(int n) { mu[1] = 1; for (int i = 2; i <= n; ++i) { if(!not_p[i]) prime[++cnt] = i, mu[i] = -1; for (int j = 1; j <= cnt && i*prime[j] <= n; ++j) { not_p[i*prime[j]] = true; if(i%prime[j] == 0) { mu[i*prime[j]] = 0; break; } mu[i*prime[j]] = -mu[i]; } } for (int i = 1; i <= n; ++i) sum[i] = (sum[i-1]+mu[i]*i*1LL*i)%MOD, c[i] = (c[i-1]+i)%MOD;}inline LL C(int x) { return (x*1LL*(x+1)/2) % MOD;}inline LL solve(int n, int m) { int up = min(n, m); LL ans = 0; for (int l = 1, r; l <= up; l = r+1) { r = min(n/(n/l), m/(m/l)); ans += (sum[r]-sum[l-1])*C(n/l)%MOD*C(m/l)%MOD; ans %= MOD; } return ans;}int main() { scanf("%d %d", &n, &m); int up = min(n, m); seive(up); LL ans = 0; for (int l = 1, r; l <= up; l = r+1) { r = min(n/(n/l), m/(m/l)); ans += (c[r]-c[l-1])*solve(n/l, m/l); ans %= MOD; } ans = (ans + MOD)%MOD; printf("%lld\n", ans); return 0;}