首先考虑相乘等于平方相当于两个数除去平方因子后相同。所以我们把所有数都除去平方因子，然后问的相当于是序列有多少排列满足相邻两项不等。为了更方便地 DP，我们对修改后的序列排序。
我们发现这个序列是可以被分成块的，满足每一块内数字相同。于是我们考虑每次插入好多个数而不止一个数。设第 \(i\) 块有 \(s_i\) 个相同的数。设 \(f(i,j)\) 表示目前插入了前 \(i\) 块，并且有 \(j\) 个位置相邻的数相同。
对于如何转移，我们考虑先枚举把第 \(i\) 块拆分成 \(k\) 个小块（这种方案数设其为 \(g(s_i,k)\)，比较好弄，之后再说）。首先其内部会产生 \(s_i-k\) 对相邻的相同的数。然后我们考虑枚举有多少个小块插入到了满足左右两边数值一样（即在插入前位置相邻且数值相同），设为 \(p\)，则对位置相邻且数值相同的对数的贡献为 \(-p\)。
设 \(q\) 为第 \(i\) 块前插入了多少个数（即块的大小总和），于是一共有 \(q+1\) 个空位可以插。
\[f(i,j+s_i-k-p)\gets \sum_{k}\binom{j}{p}\binom{q+1-j}{k-p}g(s_i,k)f(i-1,j) \]我们看看 \(g\) 的转移。注意到块内的数是有顺序而言的。
\[g(i,j)\gets (i-1+j)\times g(i-1,j)+j\times g(i-1,j-1) \]初值 \(f(0,0)=g(0,0)=1\)。
时间复杂度 \(\mathcal O(\sum s_in^2)=\mathcal O(n^3)\)。
Code：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 305, mod = 1e9 + 7;
int n;
int a[N];
int fac[N], inv[N];
int f[N][N], g[N][N];
int s[N], cnt;
int add(int a, int b) { return a + b >= mod ? a + b - mod : a + b; }
int mul(int a, int b) { return 1ll * a * b % mod; }
int qpow(int x, int y) {
int res = 1;
while (y) {
if (y & 1) res = mul(res, x);
x = mul(x, x);
y >>= 1;
}
return res;
}
void init(int maxn) {
fac[0] = 1;
for (int i = 1; i <= maxn; ++i) fac[i] = mul(fac[i - 1], i);
inv[maxn] = qpow(fac[maxn], mod - 2);
for (int i = maxn - 1; ~i; --i) inv[i] = mul(inv[i + 1], i + 1);
}
int C(int n, int m) {
if (n < 0 || m < 0 || n < m) return 0;
return mul(fac[n], mul(inv[n - m], inv[m]));
}
int div(int x) {
int t = sqrt(x);
for (int i = t; i; --i) {
if (x % (i * i) == 0) {
x /= i * i;
break;
}
}
return x;
}
int main() {
init(300);
scanf("%d", &n);
for (int i = 1; i <= n; ++i) scanf("%d", &a[i]), a[i] = div(a[i]);
sort(a + 1, a + n + 1);
for (int i = 1; i <= n; ++i) {
if (a[i] != a[i - 1]) s[++cnt] = 1;
else ++s[cnt];
}
g[0][0] = 1;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) g[i][j] = add(mul(i - 1 + j, g[i - 1][j]), mul(j, g[i - 1][j - 1]));
}
f[0][0] = 1;
for (int i = 0, sum = 0; i < cnt; ++i, sum += s[i]) {
for (int j = 0; j <= sum; ++j)
for (int k = 1; k <= s[i + 1]; ++k) {
for (int p = 0; p <= j; ++p)
f[i + 1][j + s[i + 1] - k - p] = add(f[i + 1][j + s[i + 1] - k - p], mul(mul(C(j, p), C(sum + 1 - j, k - p)), mul(g[s[i + 1]][k], f[i][j])));
}
}
printf("%d", f[cnt][0]);
return 0;
}