#include <bits/stdc++.h>
using namespace std;
int n, m;
int xmod[] = {0, 0, 1, 1};
int ymod[] = {0, 1, 0, 1};
long long cnt[] = {0, 0, 0, 0};
long long mod = 1000000007LL;
long long a[5001], s[5001];
long long A3(long long x) { return (((x * (x - 1)) % mod) * (x - 2)) % mod; }
int gcd(int x, int y) {
  if (x == 0 || y == 0) return x + y;
  return gcd(y, x % y);
}
long long grid(long long x, long long y) {
  if (x <= 0 || y <= 0) return 0;
  if (x > y) swap(x, y);
  return (s[x] + s[x - 1] + a[x] * (y - x)) % mod;
}
long long calc(long long x, long long y) {
  long long small_x = (n + 1) / x;
  long long big_x = small_x + 1;
  long long cnt_big_x = (n + 1) % x;
  long long cnt_small_x = x - cnt_big_x;
  long long small_y = (m + 1) / y;
  long long big_y = small_y + 1;
  long long cnt_big_y = (m + 1) % y;
  long long cnt_small_y = y - cnt_big_y;
  long long ans = 0;
  ans = (ans + ((cnt_small_x * cnt_small_y) % mod) * grid(small_x, small_y)) %
        mod;
  ans = (ans + ((cnt_big_x * cnt_small_y) % mod) * grid(big_x, small_y)) % mod;
  ans = (ans + ((cnt_small_x * cnt_big_y) % mod) * grid(small_x, big_y)) % mod;
  ans = (ans + ((cnt_big_x * cnt_big_y) % mod) * grid(big_x, big_y)) % mod;
  return ans;
}
int MAIN() {
  memset(a, 0, sizeof(a));
  memset(s, 0, sizeof(s));
  for (int i = 1; i <= 5000; i++) {
    a[i] = A3(i);
    s[i] = (s[i - 1] + a[i]) % mod;
  }
  cin >> n >> m;
  for (int i = 0; i <= n; i++)
    for (int j = 0; j <= m; j++)
      for (int k = 0; k < 4; k++)
        if (xmod[k] == i % 2 && ymod[k] == j % 2) cnt[k]++;
  long long ans = 0;
  for (int i = 0; i <= 1; i++)
    for (int j = 0; j <= 1; j++) {
      long long f = 0;
      for (int k = 0; k < 4; k++)
        if (xmod[k] == i % 2 && ymod[k] == j % 2) f = cnt[k], cnt[k]--;
      for (int k1 = 0; k1 < 4; k1++)
        for (int k2 = 0; k2 < 4; k2++) {
          int x1 = (i - xmod[k1] + 2) % 2;
          int y1 = (j - ymod[k1] + 2) % 2;
          int x2 = (i - xmod[k2] + 2) % 2;
          int y2 = (j - ymod[k2] + 2) % 2;
          if ((x1 * y2) % 2 == (x2 * y1) % 2) {
            ans = (ans + ((f * cnt[k1]) % mod) * cnt[k2]) % mod;
            if (k1 == k2) ans = ((ans - f * cnt[k1]) % mod + mod) % mod;
          }
        }
      for (int k = 0; k < 4; k++)
        if (xmod[k] == i % 2 && ymod[k] == j % 2) cnt[k]++;
    }
  ans -= (n + 1) * A3(m + 1);
  ans %= mod, ans += mod, ans %= mod;
  ans -= (m + 1) * A3(n + 1);
  ans %= mod, ans += mod, ans %= mod;
  for (int i = 1; i + i <= n; i++)
    for (int j = 1; j + j <= m; j++)
      if (gcd(i, j) == 1) ans = ((ans - calc(i, j) * 2) % mod + mod) % mod;
  cout << ans << endl;
  return 0;
}
int main() {
  ios ::sync_with_stdio(false);
  cout << fixed << setprecision(16);
  return MAIN();
}

### Prompt

Please create a solution in CPP to the following problem:

The Little Elephant is playing with the Cartesian coordinates' system. Most of all he likes playing with integer points. The Little Elephant defines an integer point as a pair of integers (x; y), such that 0 ≤ x ≤ w and 0 ≤ y ≤ h. Thus, the Little Elephant knows only (w + 1)·(h + 1) distinct integer points.

The Little Elephant wants to paint a triangle with vertexes at integer points, the triangle's area must be a positive integer. For that, he needs to find the number of groups of three points that form such triangle. At that, the order of points in a group matters, that is, the group of three points (0;0), (0;2), (2;2) isn't equal to the group (0;2), (0;0), (2;2).

Help the Little Elephant to find the number of groups of three integer points that form a nondegenerate triangle with integer area.

Input

A single line contains two integers w and h (1 ≤ w, h ≤ 4000).

Output

In a single output line print an integer — the remainder of dividing the answer to the problem by 1000000007 (109 + 7).

Examples

Input

2 1


Output

36


Input

2 2


Output

240

### Response

```cpp
#include <bits/stdc++.h>
using namespace std;
int n, m;
int xmod[] = {0, 0, 1, 1};
int ymod[] = {0, 1, 0, 1};
long long cnt[] = {0, 0, 0, 0};
long long mod = 1000000007LL;
long long a[5001], s[5001];
long long A3(long long x) { return (((x * (x - 1)) % mod) * (x - 2)) % mod; }
int gcd(int x, int y) {
  if (x == 0 || y == 0) return x + y;
  return gcd(y, x % y);
}
long long grid(long long x, long long y) {
  if (x <= 0 || y <= 0) return 0;
  if (x > y) swap(x, y);
  return (s[x] + s[x - 1] + a[x] * (y - x)) % mod;
}
long long calc(long long x, long long y) {
  long long small_x = (n + 1) / x;
  long long big_x = small_x + 1;
  long long cnt_big_x = (n + 1) % x;
  long long cnt_small_x = x - cnt_big_x;
  long long small_y = (m + 1) / y;
  long long big_y = small_y + 1;
  long long cnt_big_y = (m + 1) % y;
  long long cnt_small_y = y - cnt_big_y;
  long long ans = 0;
  ans = (ans + ((cnt_small_x * cnt_small_y) % mod) * grid(small_x, small_y)) %
        mod;
  ans = (ans + ((cnt_big_x * cnt_small_y) % mod) * grid(big_x, small_y)) % mod;
  ans = (ans + ((cnt_small_x * cnt_big_y) % mod) * grid(small_x, big_y)) % mod;
  ans = (ans + ((cnt_big_x * cnt_big_y) % mod) * grid(big_x, big_y)) % mod;
  return ans;
}
int MAIN() {
  memset(a, 0, sizeof(a));
  memset(s, 0, sizeof(s));
  for (int i = 1; i <= 5000; i++) {
    a[i] = A3(i);
    s[i] = (s[i - 1] + a[i]) % mod;
  }
  cin >> n >> m;
  for (int i = 0; i <= n; i++)
    for (int j = 0; j <= m; j++)
      for (int k = 0; k < 4; k++)
        if (xmod[k] == i % 2 && ymod[k] == j % 2) cnt[k]++;
  long long ans = 0;
  for (int i = 0; i <= 1; i++)
    for (int j = 0; j <= 1; j++) {
      long long f = 0;
      for (int k = 0; k < 4; k++)
        if (xmod[k] == i % 2 && ymod[k] == j % 2) f = cnt[k], cnt[k]--;
      for (int k1 = 0; k1 < 4; k1++)
        for (int k2 = 0; k2 < 4; k2++) {
          int x1 = (i - xmod[k1] + 2) % 2;
          int y1 = (j - ymod[k1] + 2) % 2;
          int x2 = (i - xmod[k2] + 2) % 2;
          int y2 = (j - ymod[k2] + 2) % 2;
          if ((x1 * y2) % 2 == (x2 * y1) % 2) {
            ans = (ans + ((f * cnt[k1]) % mod) * cnt[k2]) % mod;
            if (k1 == k2) ans = ((ans - f * cnt[k1]) % mod + mod) % mod;
          }
        }
      for (int k = 0; k < 4; k++)
        if (xmod[k] == i % 2 && ymod[k] == j % 2) cnt[k]++;
    }
  ans -= (n + 1) * A3(m + 1);
  ans %= mod, ans += mod, ans %= mod;
  ans -= (m + 1) * A3(n + 1);
  ans %= mod, ans += mod, ans %= mod;
  for (int i = 1; i + i <= n; i++)
    for (int j = 1; j + j <= m; j++)
      if (gcd(i, j) == 1) ans = ((ans - calc(i, j) * 2) % mod + mod) % mod;
  cout << ans << endl;
  return 0;
}
int main() {
  ios ::sync_with_stdio(false);
  cout << fixed << setprecision(16);
  return MAIN();
}

```