#include "bits/stdc++.h"
#include "atcoder/all"
using namespace std;
using namespace atcoder;
using mint = modint1000000007;
const int mod = 1000000007;
//using mint = modint998244353;
//const int mod = 998244353;
//const int INF = 1e9;
//const long long LINF = 1e18;
//const bool debug = false;
#define rep(i, n) for (int i = 0; i < (n); ++i)
#define rep2(i,l,r)for(int i=(l);i<(r);++i)
#define rrep(i, n) for (int i = (n-1); i >= 0; --i)
#define rrep2(i,l,r)for(int i=(r-1);i>=(l);--i)
#define all(x) (x).begin(),(x).end()
#define allR(x) (x).rbegin(),(x).rend()
#define endl "\n"
#define P pair<int,int>
template<typename A, typename B> inline bool chmax(A & a, const B & b) { if (a < b) { a = b; return true; } return false; }
template<typename A, typename B> inline bool chmin(A & a, const B & b) { if (a > b) { a = b; return true; } return false; }
// f(k) = sum_{GCD(i, j) = k} g(i)h(j)
// F(k) = sum_{k | i} f(i)
template<class T> struct FastGCDConvolution {
	int n;
	FastGCDConvolution(int n) : n(n) {}

	void zeta(vector<mint> &v) {
		int n = v.size();
		for (int k = 1; k < n; ++k) {
			for (int i = k * 2; i < n; i += k) {
				v[k] += v[i];
			}
		}
	}

	void mebius(vector<T> &v) {
		int n = v.size();
		for (int k = n - 1; k >= 1; --k) {
			for (int i = k * 2; i < n; i += k) {
				v[k] -= v[i];
			}
		}
	}

	vector<T> mult(vector<T> &a, vector<T> &b) {
		vector<T> c(n, 0);
		//vector<T> subc(n, 0);
		//rep2(i, 1, a.size())subc[i] += a[i] * b[0];
		//rep2(i, 1, b.size())subc[i] += a[0] * b[i];
		zeta(a); zeta(b);
		for (int i = 0; i < min(a.size(), b.size()); ++i) c[i] = a[i] * b[i];
		mebius(c);
		//rep(i, n)c[i] += subc[i];
		return c;
	}
};
long long gcd(long long a, long long b) {
	if (b == 0) return a;
	else return gcd(b, a%b);
}
int main() {
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n, m; cin >> n >> m;
	mint nm = ((mint)n) * m;
	mint ans = nm * (nm - 1)*(nm - 2) / 6;
	vector<mint>v0(n), v1(m);
	rep(i, n) v0[i] = n - i;
	rep(i, m) v1[i] = m - i;
	FastGCDConvolution<mint> fz(max(n, m));
	auto v = fz.mult(v0, v1);
	rep2(i, 1, max(n, m)) ans -= v[i] * (i - 1) * 2;
	ans -= ((mint)n)*m*(m - 1)*(m - 2) / 6;
	ans -= ((mint)m)*n*(n - 1)*(n - 2) / 6;
	cout << ans.val() << endl;
	return 0;
}