#include <bits/stdc++.h>
using namespace std;
typedef long long lint;
typedef long double llf;
typedef pair<int, int> pi;
const int mod = 1e9 + 7;

int n, a[100005];
int indeg[100005];
int cnt[100005];
vector<int> gph[100005];
vector<int> v;

int dfs2(int x, int p){
	vector<int> v;
	for(auto &i : gph[x]){
		if(indeg[i] || i == p) continue;
		v.push_back(dfs2(i, x) + 1);
	}
	if(v.size() >= 2) return 1e8;
	if(v.size() == 0) return 0;
	if(v.size() == 1) return v[0];
}

lint ipow(lint x, lint p){
	lint ret = 1, piv = x % mod;
	while(p){
		if(p&1) ret *= piv;
		piv *= piv;
		ret %= mod;
		piv %= mod;
		p >>= 1;
	}
	return ret % mod;
}

int rv[100005], lb[100005];

lint solve(int x){
	queue<int> que;
	for(auto &i : v){
		if(indeg[i] == 0) que.push(i);
	}
	while(!que.empty()){
		int x = que.front();
		que.pop();
		indeg[a[x]]--;
		if(indeg[a[x]] == 0) que.push(a[x]);
	}
	vector<int> w;
	lint ret = 1;
	for(auto &i : v){
		if(indeg[i] == 1){
			w.push_back(i);
			for(int j=a[i]; j!=i; j=a[j]){
				w.push_back(j);
			}
			reverse(w.begin(), w.end());
			for(int i=0; i<w.size(); i++){
				rv[w[i]] = i;
				lb[w[i]] = dfs2(w[i], -1);
				if(lb[w[i]] > 1e7) return 0;
			}
			for(int i=0; i<w.size(); i++){
				if(lb[w[i]] > 0){
					int brk = 0;
					for(int j=i+1; j<=i+w.size(); j++){
						if(lb[w[j%w.size()]] > 0){
							brk = j - i;
							break;
						}
					}
					if(brk < lb[w[i]]) return 0;
					if(brk > lb[w[i]]) ret *= 2;
					ret %= mod;
				}
			}
			break;
		}
	}
	return ret;
}

int vis[100005];
void dfs(int x){
	v.push_back(x);
	vis[x] = 1;
	for(auto &i : gph[x]){
		if(!vis[i]) dfs(i);
	}
}

lint fact[100005], invf[100005], invp[100005];

lint match(int x, int k){
	return (fact[x] * invf[x-2*k] % mod) * (invf[k] * invp[k] % mod) % mod;
}


int main(){
	scanf("%d",&n);
	fact[0] = invf[0] = invp[0] = 1;
	for(int i=1; i<=n; i++){
		fact[i] = fact[i-1] * i % mod;
		invf[i] = ipow(fact[i], mod-2);
		invp[i] = invp[i-1] * (mod + 1) / 2;
		invp[i] %= mod;
	}
	for(int i=1; i<=n; i++){
		scanf("%d",&a[i]);
		indeg[a[i]]++;
		gph[i].push_back(a[i]);
		gph[a[i]].push_back(i);
	}
	lint ret = 1;
	for(int i=1; i<=n; i++){
		if(vis[i]) continue;
		v.clear();
		dfs(i);
		int notsimple = 0;
		for(auto &i : v){
			if(indeg[i] != 1){
				notsimple = 1;
			}
		}
		if(!notsimple){
			cnt[v.size()]++;
		}
		else{
			lint w = solve(i);
			ret *= w;
			ret %= mod;
		}
	}
	for(int i=1; i<=n; i++){
		if(cnt[i] >= 1){
			lint tmp = 0;
			for(int j=0; j<=cnt[i]/2; j++){
				lint tmp2 = match(cnt[i], j) * ipow(i, j)% mod;
				if(i > 2 && i % 2 == 1) tmp2 *= ipow(2, cnt[i] - 2 * j) % mod;
				tmp += tmp2;
				tmp %= mod;
			}
			ret *= tmp;
			ret %= mod;
		}
	}
	cout << ret;
}
