MOD = 10**9 + 7 inv2 = pow(2, MOD - 2, MOD) def compute_mobius(max_d): max_d += 1 is_prime = [True] * max_d mu = [1] * max_d min_prime = [0] * max_d for i in range(2, max_d): if is_prime[i]: min_prime[i] = i for j in range(i, max_d, i): is_prime[j] = False if min_prime[j] == 0: min_prime[j] = i for i in range(2, max_d): if i == 1: continue p = min_prime[i] if p == 0: continue cnt = 0 tmp = i while tmp % p == 0: tmp //= p cnt += 1 if cnt >= 2: mu[i] = 0 else: mu[i] = -mu[tmp] return mu n, m = map(int, input().split()) def comb(k): if k < 3: return 0 return k * (k - 1) * (k - 2) // 6 % MOD total = comb(n * m) horiz = (n * comb(m)) % MOD vert = (m * comb(n)) % MOD G = min(n - 1, m - 1) if (n >= 1 and m >= 1) else 0 max_mobius = max(n, m) if G == 0 else max((n - 1) // 1, (m - 1) // 1) mu = compute_mobius(max_mobius) c_diag = 0 for g in range(1, G + 1): d_max = min((m - 1) // g, (n - 1) // g) if d_max < 1: continue S_g = 0 for d in range(1, d_max + 1): mud = mu[d] if mud == 0: continue gd = g * d A = (m - 1) // gd B = (n - 1) // gd sum_a_part1 = (m % MOD) * (A % MOD) % MOD sum_a_part2 = (gd % MOD) * (A % MOD) % MOD sum_a_part2 = sum_a_part2 * ((A + 1) % MOD) % MOD sum_a_part2 = sum_a_part2 * inv2 % MOD sum_a = (sum_a_part1 - sum_a_part2) % MOD sum_b_part1 = (n % MOD) * (B % MOD) % MOD sum_b_part2 = (gd % MOD) * (B % MOD) % MOD sum_b_part2 = sum_b_part2 * ((B + 1) % MOD) % MOD sum_b_part2 = sum_b_part2 * inv2 % MOD sum_b = (sum_b_part1 - sum_b_part2) % MOD term = mud * sum_a % MOD term = term * sum_b % MOD S_g = (S_g + term) % MOD contribution = (g - 1) * S_g % MOD contribution = contribution * 2 % MOD c_diag = (c_diag + contribution) % MOD collinear = (horiz + vert + c_diag) % MOD ans = (total - collinear) % MOD print(ans if ans >= 0 else ans + MOD)