/* Count # DAGs after shrinking cycles, and accounting for weights. */
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
string to_string(string s) { return '"' + s + '"'; }
string to_string(const char* s) { return to_string(string(s)); }
string to_string(bool b) { return to_string(int(b)); }
string to_string(vector<bool>::reference b) { return to_string(int(b)); }
string to_string(char b) { return "'" + string(1, b) + "'"; }
template <typename A, typename B>
string to_string(pair<A, B> p) { return "(" + to_string(p.first) + ", " + to_string(p.second) + ")"; }
template <typename A>
string to_string(A v) {
	string res = "{";
	for (const auto& x : v) res += (res == "{" ? "" : ", ") + to_string(x);
	return res + "}";
}
void debug() { cerr << endl; }
template <typename Head, typename... Tail>
void debug(Head H, Tail... T) {
	cerr << " " << to_string(H);
	debug(T...);
}
#define db(...) cerr << "[" << #__VA_ARGS__ << "]:", debug(__VA_ARGS__)
#else
#define db(...) 42
#endif
typedef long long ll;
typedef long double ld;
const int MOD = 1000000007;
struct Mint {
	int val;
	Mint() { val = 0; }
	Mint(ll x) {
		val = (-MOD <= x && x < MOD) ? x : x % MOD;
		if (val < 0) val += MOD;
	}
	template <typename U>
	explicit operator U() const { return (U)val; }
	friend bool operator==(const Mint& a, const Mint& b) { return a.val == b.val; }
	friend bool operator!=(const Mint& a, const Mint& b) { return !(a == b); }
	friend bool operator<(const Mint& a, const Mint& b) { return a.val < b.val; }
	Mint& operator+=(const Mint& m) { if ((val += m.val) >= MOD) val -= MOD; return *this; }
	Mint& operator-=(const Mint& m) { if ((val -= m.val) < 0) val += MOD; return *this; }
	Mint& operator*=(const Mint& m) { val = (ll)val * m.val % MOD; return *this; }
	friend Mint modex(Mint a, ll p) {
		assert(p >= 0);
		Mint ans = 1;
		for (; p; p >>= 1, a *= a) if (p & 1) ans *= a;
		return ans;
	}
	Mint& operator/=(const Mint& m) { return *this *= modex(m, MOD - 2); }
	Mint& operator++() { return *this += 1; }
	Mint& operator--() { return *this -= 1; }
	Mint operator++(int) { Mint result(*this); *this += 1; return result; }
	Mint operator--(int) { Mint result(*this); *this -= 1; return result; }
	Mint operator-() const { return Mint(-val); }
	friend Mint operator+(Mint a, const Mint& b) { return a += b; }
	friend Mint operator-(Mint a, const Mint& b) { return a -= b; }
	friend Mint operator*(Mint a, const Mint& b) { return a *= b; }
	friend Mint operator/(Mint a, const Mint& b) { return a /= b; }
	friend ostream& operator<<(ostream& os, const Mint& x) { return os << x.val; }
	friend string to_string(const Mint& b) { return to_string(b.val); }
};
vector<int> getCycle(int x, const vector<int>& perm) {
	int init = x;
	vector<int> ret;
	do {
		ret.push_back(x);
		x = perm[x];
	} while (x != init);
	return ret;
}
vector<Mint> fac(1, 1), invfac(1, 1);
Mint binom(int n, int k) {
	if (k < 0 || k > n) return 0;
	while (fac.size() <= n) {
		fac.push_back(fac.back() * fac.size());
		invfac.push_back(1 / fac.back());
	}
	return fac[n] * invfac[k] * invfac[n - k];
}
int main() {
	binom(40, 0);
	int n;
	scanf("%d", &n);
	vector<int> P(n);
	for (int i = 0; i < n; ++i) {
		scanf("%d", &P[i]);
		--P[i];
	}
	vector<bool> used(n);
	vector<int> A;
	for (int i = 0; i < n; ++i) {
		if (used[i]) continue;
		auto cyc = getCycle(i, P);
		for (int x : cyc) used[x] = true;
		A.push_back(cyc.size());
	}
	sort(A.begin(), A.end());
	auto lens = A;
	lens.resize(unique(lens.begin(), lens.end()) - lens.begin());
	int m = lens.size();
	vector<int> lenCnt(m);
	for (int x : A) {
		lenCnt[lower_bound(lens.begin(), lens.end(), x) - lens.begin()]++;
	}
	db(lens, lenCnt);
	vector<int> wlen(m + 1, 1);
	for (int i = 1; i <= m; ++i) {
		wlen[i] = wlen[i - 1] * (lenCnt[i - 1] + 1);
	}
	db(wlen[m]);
	int ns = wlen[m];
	vector<vector<int>> subs(ns);
	for (int s = 0; s < ns; ++s) {
		function<void(int, int)> dfs = [&](int p, int sub) {
			if (p == m) {
				if (sub) subs[s].push_back(sub);
				return;
			}
			int lim = s / wlen[p] % (lenCnt[p] + 1);
			for (int i = 0; i <= lim; ++i) {
				dfs(p + 1, sub + wlen[p] * i);
			}
		};
		dfs(0, 0);
	}
	vector<Mint> dp(ns);
	dp[0] = 1;
	for (int s = 1; s < ns; ++s) {
		for (int ss : subs[s]) {
			db(s, ss);
			Mint add = -1;
			int sumB = 0, sumC = 0;
			for (int i = 0; i < m; ++i) {
				int ai = s / wlen[i] % (lenCnt[i] + 1);
				int bi = ss / wlen[i] % (lenCnt[i] + 1);
				assert(bi <= ai);
				sumB += lens[i] * bi;
				sumC += lens[i] * (ai - bi);
				add *= binom(ai, bi);
				if (bi & 1) add = -add;
			}
			assert(sumB > 0);
			// Edges
			Mint emult = 0;
			for (int i = 0; i <= sumC; ++i) {
				emult += binom(sumC, i) * fac[i] * fac[n - i - 1];
			}
			add *= modex(emult, sumB);
			dp[s] += dp[s - ss] * add;
		}
	}
	printf("%d\n", dp[ns - 1].val);
}

