#include <bits/stdc++.h>
#pragma GCC diagnostic ignored "-Wunused-result"
using namespace std;
template <class T>
void _R(T &x) {
  cin >> x;
}
void _R(int &x) { scanf("%d", &x); }
void _R(int64_t &x) { scanf("%" PRId64, &x); }
void _R(double &x) { scanf("%lf", &x); }
void _R(char &x) { scanf(" %c", &x); }
void _R(char *x) { scanf("%s", x); }
void R() {}
template <class T, class... U>
void R(T &head, U &...tail) {
  _R(head);
  R(tail...);
}
template <class T>
void _W(const T &x) {
  cout << x;
}
void _W(const int &x) { printf("%d", x); }
void _W(const int64_t &x) { printf("%" PRId64, x); }
void _W(const double &x) { printf("%.16f\n", x); }
void _W(const char &x) { putchar(x); }
void _W(const char *x) { printf("%s", x); }
template <class T>
void _W(const vector<T> &x) {
  for (auto i = x.begin(); i != x.end(); _W(*i++))
    if (i != x.cbegin()) putchar(' ');
}
void W() {}
template <class T, class... U>
void W(const T &head, const U &...tail) {
  _W(head);
  putchar(sizeof...(tail) ? ' ' : '\n');
  W(tail...);
}
template <class T>
inline bool chkmax(T &a, const T &b) {
  return b > a ? a = b, true : false;
}
template <class T>
inline bool chkmin(T &a, const T &b) {
  return b < a ? a = b, true : false;
}
template <class T, class F = less<T>>
void sort_uniq(vector<T> &v, F f = F()) {
  sort(begin(v), end(v), f);
  v.resize(unique(begin(v), end(v)) - begin(v));
}
template <class T>
using MaxHeap = priority_queue<T>;
template <class T>
using MinHeap = priority_queue<T, vector<T>, greater<T>>;
template <int MOD>
struct ModInt {
  int x;
  ModInt() : x(0) {}
  ModInt(int _x) {
    x = _x % MOD;
    if (x < 0) x += MOD;
  }
  ModInt(long long _x) {
    x = _x % MOD;
    if (x < 0) x += MOD;
  }
  ModInt operator-() const { return {x == 0 ? 0 : MOD - x}; }
  ModInt &operator+=(ModInt rhs) {
    x += rhs.x;
    if (x >= MOD) x -= MOD;
    return *this;
  }
  ModInt &operator-=(ModInt rhs) {
    x -= rhs.x;
    if (x < 0) x += MOD;
    return *this;
  }
  ModInt &operator*=(ModInt rhs) {
    x = (long long)x * rhs.x % MOD;
    return *this;
  }
  ModInt &operator/=(ModInt rhs) { return *this *= rhs.inv(); }
  ModInt operator+(ModInt rhs) const { return ModInt(*this) += rhs; }
  ModInt operator-(ModInt rhs) const { return ModInt(*this) -= rhs; }
  ModInt operator*(ModInt rhs) const { return ModInt(*this) *= rhs; }
  ModInt operator/(ModInt rhs) const { return ModInt(*this) /= rhs; }
  ModInt inv() const {
    int a = x, b = MOD, u = 1, v = 0;
    while (b != 0) {
      int t = a / b;
      a -= t * b;
      u -= t * v;
      swap(a, b);
      swap(u, v);
    }
    return u;
  }
  ModInt pow(long long e) {
    ModInt r = 1, p = *this;
    while (e) {
      if (e & 1) r *= p;
      p *= p;
      e >>= 1;
    }
    return r;
  }
  bool operator==(ModInt rhs) { return x == rhs.x; }
  bool operator!=(ModInt rhs) { return x != rhs.x; }
  bool operator<(ModInt rhs) { return x < rhs.x; }
  bool operator<=(ModInt rhs) { return x <= rhs.x; }
  bool operator>(ModInt rhs) { return x > rhs.x; }
  bool operator>=(ModInt rhs) { return x >= rhs.x; }
  friend ostream &operator<<(ostream &os, ModInt i) { return os << i.x; }
};
const int MOD = 1e9 + 7;
using mint = ModInt<MOD>;
const int N = 310;
mint inv[N], fac[N], ifac[N], C[N][N];
void predo() {
  inv[1] = 1;
  for (int i = (2); i <= int(N - 1); i++)
    inv[i] = inv[MOD % i] * (MOD - MOD / i);
  fac[0] = 1;
  for (int i = (1); i <= int(N - 1); i++) fac[i] = fac[i - 1] * i;
  ifac[0] = 1;
  for (int i = (1); i <= int(N - 1); i++) ifac[i] = ifac[i - 1] * inv[i];
  C[0][0] = 1;
  for (int i = (1); i <= int(N - 1); i++) {
    C[i][0] = 1;
    for (int j = (1); j <= int(i); j++) C[i][j] = C[i - 1][j] + C[i - 1][j - 1];
  }
}
int n, a[N], b[N], c[N];
int f(int x) {
  static map<vector<int>, int> mp;
  vector<int> v;
  for (int i = 2; i * i <= x; i++) {
    int cnt = 0;
    while (x % i == 0) {
      x /= i;
      cnt ^= 1;
    }
    if (cnt) v.push_back(i);
  }
  if (x > 1) v.push_back(x);
  if (mp.count(v)) return mp[v];
  auto id = ((int)(mp).size()) + 1;
  return mp[v] = id;
}
mint _dp[2][N];
auto dp = _dp[0], dq = _dp[1];
int main() {
  predo();
  R(n);
  for (int i = 0; i < int(n); i++) R(a[i]);
  for (int i = 0; i < int(n); i++) b[i] = f(a[i]);
  sort(b, b + n);
  int m = 0;
  for (int i = 0, j; i < n; i = j) {
    for (j = i; j < n && b[i] == b[j]; j++)
      ;
    c[m++] = j - i;
  }
  sort(c, c + m);
  dp[0] = 1;
  int s = 1;
  for (int i = 0; i < int(m); i++) {
    memset(dq, 0, N * sizeof(mint));
    for (int j = (0); j <= int(s); j++) {
      mint me = dp[j];
      if (me.x == 0) continue;
      int uk = min(c[i], s);
      for (int k = (1); k <= int(uk); k++) {
        int lx = k - j;
        int ux = min(k, s - j);
        for (int x = (lx); x <= int(ux); x++) {
          int y = k - x;
          dq[j - y + (c[i] - k)] +=
              C[s - j][x] * C[j][y] * C[c[i] - 1][k - 1] * me;
        }
      }
    }
    s += c[i];
    swap(dp, dq);
  }
  mint ans = dp[0];
  for (int i = 0; i < int(m); i++) ans *= fac[c[i]];
  W(ans);
  return 0;
}

### Prompt

Generate a Cpp solution to the following problem:

A year ago on the bench in public park Leha found an array of n numbers. Leha believes that permutation p is right if for all 1 ≤ i < n condition, that api·api + 1 is not perfect square, holds. Leha wants to find number of right permutations modulo 109 + 7.

Input

First line of input data contains single integer n (1 ≤ n ≤ 300) — length of the array.

Next line contains n integers a1, a2, ... , an (1 ≤ ai ≤ 109) — found array.

Output

Output single integer — number of right permutations modulo 109 + 7.

Examples

Input

3
1 2 4


Output

2


Input

7
5 2 4 2 4 1 1


Output

144

Note

For first example:

[1, 2, 4] — right permutation, because 2 and 8 are not perfect squares.

[1, 4, 2] — wrong permutation, because 4 is square of 2.

[2, 1, 4] — wrong permutation, because 4 is square of 2.

[2, 4, 1] — wrong permutation, because 4 is square of 2.

[4, 1, 2] — wrong permutation, because 4 is square of 2.

[4, 2, 1] — right permutation, because 8 and 2 are not perfect squares.

### Response

```cpp
#include <bits/stdc++.h>
#pragma GCC diagnostic ignored "-Wunused-result"
using namespace std;
template <class T>
void _R(T &x) {
  cin >> x;
}
void _R(int &x) { scanf("%d", &x); }
void _R(int64_t &x) { scanf("%" PRId64, &x); }
void _R(double &x) { scanf("%lf", &x); }
void _R(char &x) { scanf(" %c", &x); }
void _R(char *x) { scanf("%s", x); }
void R() {}
template <class T, class... U>
void R(T &head, U &...tail) {
  _R(head);
  R(tail...);
}
template <class T>
void _W(const T &x) {
  cout << x;
}
void _W(const int &x) { printf("%d", x); }
void _W(const int64_t &x) { printf("%" PRId64, x); }
void _W(const double &x) { printf("%.16f\n", x); }
void _W(const char &x) { putchar(x); }
void _W(const char *x) { printf("%s", x); }
template <class T>
void _W(const vector<T> &x) {
  for (auto i = x.begin(); i != x.end(); _W(*i++))
    if (i != x.cbegin()) putchar(' ');
}
void W() {}
template <class T, class... U>
void W(const T &head, const U &...tail) {
  _W(head);
  putchar(sizeof...(tail) ? ' ' : '\n');
  W(tail...);
}
template <class T>
inline bool chkmax(T &a, const T &b) {
  return b > a ? a = b, true : false;
}
template <class T>
inline bool chkmin(T &a, const T &b) {
  return b < a ? a = b, true : false;
}
template <class T, class F = less<T>>
void sort_uniq(vector<T> &v, F f = F()) {
  sort(begin(v), end(v), f);
  v.resize(unique(begin(v), end(v)) - begin(v));
}
template <class T>
using MaxHeap = priority_queue<T>;
template <class T>
using MinHeap = priority_queue<T, vector<T>, greater<T>>;
template <int MOD>
struct ModInt {
  int x;
  ModInt() : x(0) {}
  ModInt(int _x) {
    x = _x % MOD;
    if (x < 0) x += MOD;
  }
  ModInt(long long _x) {
    x = _x % MOD;
    if (x < 0) x += MOD;
  }
  ModInt operator-() const { return {x == 0 ? 0 : MOD - x}; }
  ModInt &operator+=(ModInt rhs) {
    x += rhs.x;
    if (x >= MOD) x -= MOD;
    return *this;
  }
  ModInt &operator-=(ModInt rhs) {
    x -= rhs.x;
    if (x < 0) x += MOD;
    return *this;
  }
  ModInt &operator*=(ModInt rhs) {
    x = (long long)x * rhs.x % MOD;
    return *this;
  }
  ModInt &operator/=(ModInt rhs) { return *this *= rhs.inv(); }
  ModInt operator+(ModInt rhs) const { return ModInt(*this) += rhs; }
  ModInt operator-(ModInt rhs) const { return ModInt(*this) -= rhs; }
  ModInt operator*(ModInt rhs) const { return ModInt(*this) *= rhs; }
  ModInt operator/(ModInt rhs) const { return ModInt(*this) /= rhs; }
  ModInt inv() const {
    int a = x, b = MOD, u = 1, v = 0;
    while (b != 0) {
      int t = a / b;
      a -= t * b;
      u -= t * v;
      swap(a, b);
      swap(u, v);
    }
    return u;
  }
  ModInt pow(long long e) {
    ModInt r = 1, p = *this;
    while (e) {
      if (e & 1) r *= p;
      p *= p;
      e >>= 1;
    }
    return r;
  }
  bool operator==(ModInt rhs) { return x == rhs.x; }
  bool operator!=(ModInt rhs) { return x != rhs.x; }
  bool operator<(ModInt rhs) { return x < rhs.x; }
  bool operator<=(ModInt rhs) { return x <= rhs.x; }
  bool operator>(ModInt rhs) { return x > rhs.x; }
  bool operator>=(ModInt rhs) { return x >= rhs.x; }
  friend ostream &operator<<(ostream &os, ModInt i) { return os << i.x; }
};
const int MOD = 1e9 + 7;
using mint = ModInt<MOD>;
const int N = 310;
mint inv[N], fac[N], ifac[N], C[N][N];
void predo() {
  inv[1] = 1;
  for (int i = (2); i <= int(N - 1); i++)
    inv[i] = inv[MOD % i] * (MOD - MOD / i);
  fac[0] = 1;
  for (int i = (1); i <= int(N - 1); i++) fac[i] = fac[i - 1] * i;
  ifac[0] = 1;
  for (int i = (1); i <= int(N - 1); i++) ifac[i] = ifac[i - 1] * inv[i];
  C[0][0] = 1;
  for (int i = (1); i <= int(N - 1); i++) {
    C[i][0] = 1;
    for (int j = (1); j <= int(i); j++) C[i][j] = C[i - 1][j] + C[i - 1][j - 1];
  }
}
int n, a[N], b[N], c[N];
int f(int x) {
  static map<vector<int>, int> mp;
  vector<int> v;
  for (int i = 2; i * i <= x; i++) {
    int cnt = 0;
    while (x % i == 0) {
      x /= i;
      cnt ^= 1;
    }
    if (cnt) v.push_back(i);
  }
  if (x > 1) v.push_back(x);
  if (mp.count(v)) return mp[v];
  auto id = ((int)(mp).size()) + 1;
  return mp[v] = id;
}
mint _dp[2][N];
auto dp = _dp[0], dq = _dp[1];
int main() {
  predo();
  R(n);
  for (int i = 0; i < int(n); i++) R(a[i]);
  for (int i = 0; i < int(n); i++) b[i] = f(a[i]);
  sort(b, b + n);
  int m = 0;
  for (int i = 0, j; i < n; i = j) {
    for (j = i; j < n && b[i] == b[j]; j++)
      ;
    c[m++] = j - i;
  }
  sort(c, c + m);
  dp[0] = 1;
  int s = 1;
  for (int i = 0; i < int(m); i++) {
    memset(dq, 0, N * sizeof(mint));
    for (int j = (0); j <= int(s); j++) {
      mint me = dp[j];
      if (me.x == 0) continue;
      int uk = min(c[i], s);
      for (int k = (1); k <= int(uk); k++) {
        int lx = k - j;
        int ux = min(k, s - j);
        for (int x = (lx); x <= int(ux); x++) {
          int y = k - x;
          dq[j - y + (c[i] - k)] +=
              C[s - j][x] * C[j][y] * C[c[i] - 1][k - 1] * me;
        }
      }
    }
    s += c[i];
    swap(dp, dq);
  }
  mint ans = dp[0];
  for (int i = 0; i < int(m); i++) ans *= fac[c[i]];
  W(ans);
  return 0;
}

```