#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;
#define endl '\n'
#define ll long long
#define pi pair<int, int>
#define f first
#define s second

const int mod = 1000000007;
const int mxn = 41, mxm = 1500;
int n;
ll a[mxn], b[mxn], c[mxn], p[mxn], f[mxn], iv[mxn], vis[mxn];
ll ff[mxn][mxn], g[mxn][mxn];
ll dp[mxm][mxn][mxn];

int main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	
	cin >> n;
	
	for(int i = 0; i < n; i++) cin >> a[i], a[i]--;
	
	for(int i = 0; i < n; i++) if(!vis[i]){
		int x = 0;
		for(int j = i; !vis[j]; j = a[j]) x++, vis[j] = 1;
		b[x]++;
	}
	
	p[n] = 1;
	for(int i = n; i; i--) p[i - 1] = (b[i] + 1) * p[i];
	
	f[0] = iv[0] = ff[0][0] = 1;
	for(int i = 1; i <= n; i++){
		iv[i] = ~-i ? mod - mod / i * iv[mod % i] % mod : 1;
		f[i] = i * f[i - 1] % mod;
		ff[i][0] = 1;
		for(int j = 1; j <= i; j++){
			ff[i][j] = (ff[i - 1][j] + ff[i - 1][j - 1]) % mod; 
		}
	}
	
	for(int i = 0; i <= n; i++){
		ll x = f[n] * iv[n - i] % mod;
		g[i][0] = 1;
		for(int j = 1; j <= n; j++) g[i][j] = x * g[i][j - 1] % mod; 
	}
	
	dp[0][0][0] = 1;
	for(int i = 0; i < p[0]; i++){ //cycle counts
		int x = i, y = 0;
		for(int j = 1; j <= n; j++) c[j] = x / p[j], x %= p[j], y += c[j] * j;
		for(int j = 0; j <= y; j++){ //rounds passed
			for(int k = 1; k <= n; k++) //cycle size added
			for(int l = 0; l <= b[k] - c[k]; l++){
				(dp[i + l * p[k]][j][k] += mod + (1 - 2 * (l & 1)) * 
					dp[i][j][k - 1] * g[j][k * l] % mod * ff[l + c[k]][c[k]]) %= mod;
			}
			if(j < y) (dp[i][y][0] += mod - dp[i][j][n]) %= mod;
		}
	}
	
	cout << dp[p[0] - 1][n][0] << endl;
	
	return 0;
}