#include <bits/stdc++.h>
const int INF = 1e9 + 1;
const double PI = acos(-1.0);
const double EPS = 1e-9;
using namespace std;
template <typename F, typename S>
ostream &operator<<(ostream &os, const pair<F, S> &p) {
  return os << "(" << p.first << ", " << p.second << ")";
}
template <typename T>
ostream &operator<<(ostream &os, const vector<T> &v) {
  os << "{";
  typename vector<T>::const_iterator it;
  for (it = v.begin(); it != v.end(); it++) {
    if (it != v.begin()) os << ", ";
    os << *it;
  }
  return os << "}";
}
template <typename T>
ostream &operator<<(ostream &os, const set<T> &v) {
  os << "[";
  typename set<T>::const_iterator it;
  for (it = v.begin(); it != v.end(); it++) {
    if (it != v.begin()) os << ", ";
    os << *it;
  }
  return os << "]";
}
template <typename F, typename S>
ostream &operator<<(ostream &os, const map<F, S> &v) {
  os << "[";
  typename map<F, S>::const_iterator it;
  for (it = v.begin(); it != v.end(); it++) {
    if (it != v.begin()) os << ", ";
    os << it->first << " = " << it->second;
  }
  return os << "]";
}
template <class T>
void deb_array(T *arr, int length) {
  for (int i = 0; i < length; i++) {
    cout << arr[i] << ' ';
  }
  cout << '\n';
}
const int maxn = 305;
const long long mod = 1e9 + 7;
int n;
long long a[maxn], fac[maxn];
long long comb[maxn][maxn], dp[maxn][maxn];
vector<int> group;
map<int, int> mp;
void calc() {
  for (int i = 0; i < maxn; i++) {
    for (int j = 0; j <= i; j++) {
      if (i == j || j == 0)
        comb[i][j] = 1;
      else
        comb[i][j] = (comb[i - 1][j - 1] + comb[i - 1][j]) % mod;
    }
  }
  fac[0] = 1;
  for (int i = 1; i < maxn; i++) {
    fac[i] = fac[i - 1] * i;
    fac[i] %= mod;
  }
}
int main() {
  ios_base::sync_with_stdio(0);
  cin >> n;
  for (int i = 0; i < n; i++) {
    cin >> a[i];
    int b = a[i], num;
    for (int j = 2; j * j <= b; j++) {
      num = 0;
      while (b % j == 0) {
        num++;
        b /= j;
      }
      if (num & 1) {
        b *= j;
      }
    }
    mp[b]++;
  }
  for (auto it = mp.begin(); it != mp.end(); it++) {
    group.push_back(it->second);
  }
  calc();
  dp[0][group[0] - 1] = fac[group[0]];
  int sum = group[0] + 1;
  for (int x = 1; x < group.size(); x++) {
    for (int y = 0; y <= sum; y++) {
      for (int ij = 1; ij <= group[x]; ij++) {
        for (int j = 0; j <= min(ij, y); j++) {
          long long cur = dp[x - 1][y];
          if (cur == 0) continue;
          cur *= fac[group[x]];
          cur %= mod;
          cur *= comb[y][j];
          cur %= mod;
          cur *= comb[sum - y][ij - j];
          cur %= mod;
          cur *= comb[group[x] - 1][ij - 1];
          dp[x][y - j + group[x] - ij] += cur;
          dp[x][y - j + group[x] - ij] %= mod;
        }
      }
    }
    sum += group[x];
  }
  cout << dp[group.size() - 1][0] << '\n';
  return 0;
}

### Prompt

Please create a solution in Cpp 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>
const int INF = 1e9 + 1;
const double PI = acos(-1.0);
const double EPS = 1e-9;
using namespace std;
template <typename F, typename S>
ostream &operator<<(ostream &os, const pair<F, S> &p) {
  return os << "(" << p.first << ", " << p.second << ")";
}
template <typename T>
ostream &operator<<(ostream &os, const vector<T> &v) {
  os << "{";
  typename vector<T>::const_iterator it;
  for (it = v.begin(); it != v.end(); it++) {
    if (it != v.begin()) os << ", ";
    os << *it;
  }
  return os << "}";
}
template <typename T>
ostream &operator<<(ostream &os, const set<T> &v) {
  os << "[";
  typename set<T>::const_iterator it;
  for (it = v.begin(); it != v.end(); it++) {
    if (it != v.begin()) os << ", ";
    os << *it;
  }
  return os << "]";
}
template <typename F, typename S>
ostream &operator<<(ostream &os, const map<F, S> &v) {
  os << "[";
  typename map<F, S>::const_iterator it;
  for (it = v.begin(); it != v.end(); it++) {
    if (it != v.begin()) os << ", ";
    os << it->first << " = " << it->second;
  }
  return os << "]";
}
template <class T>
void deb_array(T *arr, int length) {
  for (int i = 0; i < length; i++) {
    cout << arr[i] << ' ';
  }
  cout << '\n';
}
const int maxn = 305;
const long long mod = 1e9 + 7;
int n;
long long a[maxn], fac[maxn];
long long comb[maxn][maxn], dp[maxn][maxn];
vector<int> group;
map<int, int> mp;
void calc() {
  for (int i = 0; i < maxn; i++) {
    for (int j = 0; j <= i; j++) {
      if (i == j || j == 0)
        comb[i][j] = 1;
      else
        comb[i][j] = (comb[i - 1][j - 1] + comb[i - 1][j]) % mod;
    }
  }
  fac[0] = 1;
  for (int i = 1; i < maxn; i++) {
    fac[i] = fac[i - 1] * i;
    fac[i] %= mod;
  }
}
int main() {
  ios_base::sync_with_stdio(0);
  cin >> n;
  for (int i = 0; i < n; i++) {
    cin >> a[i];
    int b = a[i], num;
    for (int j = 2; j * j <= b; j++) {
      num = 0;
      while (b % j == 0) {
        num++;
        b /= j;
      }
      if (num & 1) {
        b *= j;
      }
    }
    mp[b]++;
  }
  for (auto it = mp.begin(); it != mp.end(); it++) {
    group.push_back(it->second);
  }
  calc();
  dp[0][group[0] - 1] = fac[group[0]];
  int sum = group[0] + 1;
  for (int x = 1; x < group.size(); x++) {
    for (int y = 0; y <= sum; y++) {
      for (int ij = 1; ij <= group[x]; ij++) {
        for (int j = 0; j <= min(ij, y); j++) {
          long long cur = dp[x - 1][y];
          if (cur == 0) continue;
          cur *= fac[group[x]];
          cur %= mod;
          cur *= comb[y][j];
          cur %= mod;
          cur *= comb[sum - y][ij - j];
          cur %= mod;
          cur *= comb[group[x] - 1][ij - 1];
          dp[x][y - j + group[x] - ij] += cur;
          dp[x][y - j + group[x] - ij] %= mod;
        }
      }
    }
    sum += group[x];
  }
  cout << dp[group.size() - 1][0] << '\n';
  return 0;
}

```