#include <bits/stdc++.h>
using namespace std;
const int MOD = 1000000007;
int totalCases, testNum;
int W, H;
int DP[2][2];
int gcd(int a, int b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}
int coll(int x1, int x2) {
  int g = gcd(x1, x2);
  x2 /= g;
  x1 /= g;
  if (x1 == 0) x2 = 1;
  if (x2 < x1) swap(x1, x2);
  return 2 * (H / x2) * (H + 1) - x2 * (1 + H / x2) * (H / x2) + H + 1;
}
void preprocess() {}
bool input() {
  scanf("%d", &W);
  scanf("%d", &H);
  return true;
}
void solve() {
  memset(DP, 0, sizeof(DP));
  for (int y1 = -H; y1 <= H; y1++)
    for (int y2 = -H; y2 <= H; y2++) {
      int yext = max(abs(y1), max(abs(y2), abs(y1 - y2)));
      if (yext > H) continue;
      if (abs(y1) & 1) {
        if (abs(y2) & 1) {
          DP[0][0] = (DP[0][0] + (H + 1) - yext) % MOD;
          if (DP[0][0] < 0) DP[0][0] += MOD;
          DP[1][1] = (DP[1][1] + (H + 1) - yext) % MOD;
          if (DP[1][1] < 0) DP[1][1] += MOD;
        } else {
          for (int i = 0; i < 2; i++) {
            DP[i][0] += H + 1 - yext;
            if (DP[i][0] >= MOD) DP[i][0] -= MOD;
          }
        }
      } else {
        if (abs(y2) & 1) {
          for (int i = 0; i < 2; i++) {
            DP[0][i] += H + 1 - yext;
            if (DP[0][i] >= MOD) DP[0][i] -= MOD;
          }
        } else {
          for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++) {
              DP[i][j] += H + 1 - yext;
              if (DP[i][j] >= MOD) DP[i][j] -= MOD;
            }
        }
      }
    }
  int ans = 0;
  for (int x1 = 1; x1 <= W; x1++)
    for (int x2 = 1; x2 <= W; x2++) {
      int term = DP[x1 & 1][x2 & 1];
      term -= coll(x1, x2);
      if (term < 0) term += MOD;
      ans += (term * 1ll * (W + 1 - max(x1, x2))) % MOD;
      if (ans >= MOD) ans -= MOD;
    }
  for (int x2 = 1; x2 <= W; x2++) {
    int term = DP[0][x2 & 1];
    term -= coll(0, x2);
    ans += (term * 1ll * (W + 1 - x2)) % MOD;
    if (ans >= MOD) ans -= MOD;
  }
  ans = (3ll * ans) % MOD;
  printf("%d\n", ans);
}
int main() {
  preprocess();
  totalCases = 1;
  for (testNum = 1; testNum <= totalCases; testNum++) {
    if (!input()) break;
    solve();
  }
}

### Prompt

Construct a Cpp code solution to the problem outlined:

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;
const int MOD = 1000000007;
int totalCases, testNum;
int W, H;
int DP[2][2];
int gcd(int a, int b) {
  if (b == 0) return a;
  return gcd(b, a % b);
}
int coll(int x1, int x2) {
  int g = gcd(x1, x2);
  x2 /= g;
  x1 /= g;
  if (x1 == 0) x2 = 1;
  if (x2 < x1) swap(x1, x2);
  return 2 * (H / x2) * (H + 1) - x2 * (1 + H / x2) * (H / x2) + H + 1;
}
void preprocess() {}
bool input() {
  scanf("%d", &W);
  scanf("%d", &H);
  return true;
}
void solve() {
  memset(DP, 0, sizeof(DP));
  for (int y1 = -H; y1 <= H; y1++)
    for (int y2 = -H; y2 <= H; y2++) {
      int yext = max(abs(y1), max(abs(y2), abs(y1 - y2)));
      if (yext > H) continue;
      if (abs(y1) & 1) {
        if (abs(y2) & 1) {
          DP[0][0] = (DP[0][0] + (H + 1) - yext) % MOD;
          if (DP[0][0] < 0) DP[0][0] += MOD;
          DP[1][1] = (DP[1][1] + (H + 1) - yext) % MOD;
          if (DP[1][1] < 0) DP[1][1] += MOD;
        } else {
          for (int i = 0; i < 2; i++) {
            DP[i][0] += H + 1 - yext;
            if (DP[i][0] >= MOD) DP[i][0] -= MOD;
          }
        }
      } else {
        if (abs(y2) & 1) {
          for (int i = 0; i < 2; i++) {
            DP[0][i] += H + 1 - yext;
            if (DP[0][i] >= MOD) DP[0][i] -= MOD;
          }
        } else {
          for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++) {
              DP[i][j] += H + 1 - yext;
              if (DP[i][j] >= MOD) DP[i][j] -= MOD;
            }
        }
      }
    }
  int ans = 0;
  for (int x1 = 1; x1 <= W; x1++)
    for (int x2 = 1; x2 <= W; x2++) {
      int term = DP[x1 & 1][x2 & 1];
      term -= coll(x1, x2);
      if (term < 0) term += MOD;
      ans += (term * 1ll * (W + 1 - max(x1, x2))) % MOD;
      if (ans >= MOD) ans -= MOD;
    }
  for (int x2 = 1; x2 <= W; x2++) {
    int term = DP[0][x2 & 1];
    term -= coll(0, x2);
    ans += (term * 1ll * (W + 1 - x2)) % MOD;
    if (ans >= MOD) ans -= MOD;
  }
  ans = (3ll * ans) % MOD;
  printf("%d\n", ans);
}
int main() {
  preprocess();
  totalCases = 1;
  for (testNum = 1; testNum <= totalCases; testNum++) {
    if (!input()) break;
    solve();
  }
}

```