#include<iostream>
#include<algorithm>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<string>
#include<stack>
#include<cstdio>
#include<cmath>
using namespace std;

typedef long long ll;
typedef long double ld;
typedef pair<int,int> P;
typedef pair<int,P> P1;

#define fr first
#define sc second
#define mp make_pair
#define pb push_back
#define rep(i,x) for(int i=0;i<x;i++)
#define rep1(i,x) for(int i=1;i<=x;i++)
#define rrep(i,x) for(int i=x-1;i>=0;i--)
#define rrep1(i,x) for(int i=x;i>0;i--)
#define sor(v) sort(v.begin(),v.end())
#define rev(s) reverse(s.begin(),s.end())
#define lb(vec,a) lower_bound(vec.begin(),vec.end(),a)
#define ub(vec,a) upper_bound(vec.begin(),vec.end(),a)
#define uniq(vec) vec.erase(unique(vec.begin(),vec.end()),vec.end())
#define mp1(a,b,c) P1(a,P(b,c))

const int INF=1000000000;
const int dir_4[4][2]={{1,0},{0,1},{-1,0},{0,-1}};
const int dir_8[8][2]={{1,0},{1,1},{0,1},{-1,1},{-1,0},{-1,-1},{0,-1},{1,-1}};

const ll M = 1000000007;

struct UF{
	int par[100010],r[100010];
	void init(){
		rep(i,100010){
			par[i] = i;
			r[i] = 0;
		}
	}
	int find(int x){
		if(par[x] == x)return x;
		return par[x] = find(par[x]);
	}
	bool same(int x,int y){
		return find(x) == find(y);
	}
	void unit(int x,int y){
		if(same(x,y))return;
		x = find(x);
		y = find(y);
		if(r[x] < r[y]){
			par[x] = y;
		}
		else {
			par[y] = x;
			if(r[x] == r[y]){
				r[x] ++;
			}
		}
	}
}uf;

ll POW[100010],FRA[100010];
ll POW_inv[100010],FRA_inv[100010];
ll modpow(ll x,ll k){
	if(k == 0)return 1;
	ll ret = modpow(x,k/2);
	ret *= ret; ret %= M;
	if(k%2 == 1){
		ret *= x;
		ret %= M;
	}
	return ret;
}
void init_math(){
	POW[0] = 1;
	rep1(i,100000){
		POW[i] = POW[i-1]*2;
		POW[i] %= M;
	}
	FRA[0] = 1;
	rep1(i,100000){
		FRA[i] = FRA[i-1]*i;
		FRA[i] %= M;
	}
	rep(i,100001){
		POW_inv[i] = modpow(POW[i],M-2);
		FRA_inv[i] = modpow(FRA[i],M-2);
	}
}

ll ret;
ll loop_cnt[100010];
int nex[100010];
vector<int> rG[100010];
bool loop[100010];
int hige[100010];
bool used[100010];
bool dfs_hige(int v){
	hige[v] = 0;
	for(int u: rG[v]){
		if(loop[u])continue;
		if(!dfs_hige(u))return false;
		if(hige[v] != 0)return false;
		hige[v] = hige[u];
	}
	hige[v] ++;
	return true;
}
void solve(vector<int> vs){
	/*for(int v: vs){
		cout << v << " ";
	}
	cout << endl;*/
	int v = vs[0];
	while(!used[v]){
		used[v] = true;
		v = nex[v];
	}
	while(!loop[v]){
		loop[v] = true;
		v = nex[v];
	}
	int cnt = 0;
	vector<int> vec,vec_;
	while(hige[v] == -1){
		if(!dfs_hige(v)){
			ret = 0;
			return;
		}
		cnt ++;
		vec.pb(hige[v]-1);
		vec_.pb(0);
		v = nex[v];
	}
	//puts("DEBUG");
	if(cnt == vs.size()){
		loop_cnt[cnt] ++;
		return;
	}
	//puts("DEBUG");
	/*cout << "vec(hige) check" << endl;
	rep(i,vec.size()){
		cout << vec[i] << " ";
	}
	cout << endl;*/
	rep(i,vec.size()){
		int j = vec[i];
		rep1(k,j){
			vec_[(i-k+cnt)%cnt] ++;
		}
	}
	rep(i,vec_.size()){
		if(vec_[i] > 1){
			ret = 0;
			return;
		}
	}
	//puts("DEBUG");
	ll ret_ = 1;
	rep(i,vec.size()){
		if(vec[i] > 0 && vec_[i] == 0){
			ret_ *= 2;
			ret_ %= M;
		}
	}
	ret *= ret_;
	ret %= M;
}

int main(){
	int n;
	static int a[100010];
	scanf("%d",&n);
	rep1(i,n)scanf("%d",&a[i]);
	
	init_math();
	uf.init();
	rep1(i,n){
		uf.unit(i,a[i]);
	}
	ret = 1;
	rep1(i,n){
		nex[i] = a[i];
		rG[a[i]].pb(i);
	}
	rep(i,100010){
		loop_cnt[i] = 0;
		hige[i] = -1;
		loop[i] = false;
		used[i] = false;
	}
	
	static vector<int> vs[100010];
	rep1(i,n){
		vs[uf.find(i)].pb(i);
	}
	rep1(i,n){
		if(uf.find(i) == i)solve(vs[i]);
	}
	
	/*cout << ret << endl;
	rep1(i,n){
		cout << loop_cnt[i] << " ";
	}
	cout << endl;*/
	
	for(int i = 1 ; i <= n ; i += 2){
		if(i == 3)i --;
		ll sum = 0;
		ll x = loop_cnt[i];
		for(int j = 0 ; 2*j <= x ; j ++){
			ll sum_ = FRA[x];
			sum_ *= POW_inv[j]; sum_ %= M;
			sum_ *= FRA_inv[j]; sum_ %= M;
			sum_ *= FRA_inv[x-2*j]; sum_ %= M;
			sum_ *= modpow(i,j); sum_ %= M;
			sum += sum_; if(sum >= M)sum -= M;
		}
		//cout << i << ":" << sum << endl;
		ret *= sum;
		ret %= M;
	}
	for(int i = 3 ; i <= n ; i += 2){
		ll sum = 0;
		ll x = loop_cnt[i];
		for(int j = 0 ; 2*j <= x ; j ++){
			ll sum_ = FRA[x];
			sum_ *= POW_inv[j]; sum_ %= M;
			sum_ *= FRA_inv[j]; sum_ %= M;
			sum_ *= FRA_inv[x-2*j]; sum_ %= M;
			sum_ *= POW[x-2*j]; sum_ %= M;
			sum_ *= modpow(i,j); sum_ %= M;
			sum += sum_; if(sum >= M)sum -= M;
			//cout << sum_ << endl;
		}
		//cout << i << ":" << sum << endl;
		ret *= sum;
		ret %= M;
	}
	cout << ret << endl;
}


Measured memory_bytes: 17908