Problem: You are given an integer sequence a of length N. How many permutations p of the integers 1 through N satisfy the following condition?

* For each 1 ≤ i ≤ N, at least one of the following holds: p_i = a_i and p_{p_i} = a_i.



Find the count modulo 10^9 + 7.

Constraints

* 1 ≤ N ≤ 10^5
* a_i is an integer.
* 1 ≤ a_i ≤ N

Input

The input is given from Standard Input in the following format:


N
a_1 a_2 ... a_N


Output

Print the number of the permutations p that satisfy the condition, modulo 10^9 + 7.

Examples

Input

3
1 2 3


Output

4


Input

2
1 1


Output

1


Input

3
2 1 1


Output

2


Input

3
1 1 1


Output

0


Input

13
2 1 4 3 6 7 5 9 10 8 8 9 11


Output

6

Solution (in CPP):

Solution:
#include <bits/stdc++.h>
using namespace std;

const int RLEN=1<<18|1;
inline char nc() {
	static char ibuf[RLEN],*ib,*ob;
	(ib==ob) && (ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
	return (ib==ob) ? -1 : *ib++;
}
inline int rd() {
	char ch=nc(); int i=0,f=1;
	while(!isdigit(ch)) {if(ch=='-')f=-1; ch=nc();}
	while(isdigit(ch)) {i=(i<<1)+(i<<3)+ch-'0'; ch=nc();}
	return i*f;
}

const int N=1e5+50, mod=1e9+7;

inline int add(int x,int y) {return (x+y>=mod) ? (x+y-mod) : (x+y);}
inline int dec(int x,int y) {return (x-y<0) ? (x-y+mod) : (x-y);}
inline int mul(int x,int y) {return (long long)x*y%mod;}
inline int power(int a,int b,int rs=1) {for(;b;b>>=1,a=mul(a,a)) if(b&1) rs=mul(rs,a); return rs;}
inline int cinv(int a) {return power(a,mod-2);}

int n,ans=1,fa[N],incir[N],vis[N];
int cnt[N];
vector <int> edge[N];

struct combin {
	int fac[N],ifac[N];
	combin() {
		fac[0]=1;
		for(int i=1;i<N;i++) fac[i]=mul(fac[i-1],i);
		ifac[N-1]=cinv(fac[N-1]);
		for(int i=N-2;~i;i--) ifac[i]=mul(ifac[i+1],i+1);
	}
	inline int C(int a,int b) {return (a<b) ? 0 : (mul(fac[a],mul(ifac[b],ifac[a-b])));}
} C;

inline void dfs(int x) {
	vis[x]=1; 
	for(auto v:edge[x]) {
		if(!vis[v]) fa[v]=x, dfs(v);
		else {
			for(int u=x;;u=fa[u]) {
				if(!incir[u]) incir[u]=1;
				else {puts("0"); exit(0);}
				if(u==v) break;
			}
		}
	}
}

inline int calc_dep(int x) {
	if(edge[x].size()>=2) return -1e9;
	if(!edge[x].size()) return 1;
	return calc_dep(edge[x][0])+1;
}

vector <int> cir;
vector <int> ft;
int pos[N],dep[N],f[N];

inline int calc(int x) {
	cir.clear(); ft.clear();
	int now=x;
	while(1) {
		cir.push_back(now); pos[now]=cir.size();
		vis[now]=1;
		if(edge[now].size()>=3) return 0; 
		for(auto v:edge[now]) {
			if(!incir[v]) {
				int t=calc_dep(v);
				if(t>0) ft.push_back(now), dep[now]=t;
				else return 0;
			}
		}
		for(auto v:edge[now]) if(incir[v]) {now=v; break;}
		if(now==x) break;
	}
	if(!ft.size()) return ++cnt[cir.size()],1;
	else {
		int sum=1;
		for(int i=0;i<ft.size();++i) {
			int dis=(ft.size()==1) ? cir.size() : (ft[i]==ft.back() ? cir.size()-pos[ft[i]]+pos[ft[0]] : pos[ft[i+1]]-pos[ft[i]]);
			if(dis<dep[ft[i]]) return 0;
			else if(dis>dep[ft[i]]) sum=mul(sum,2); 
		} return sum;
	}
}

int main() {
	n=rd();
	for(int i=1;i<=n;i++) 
		edge[rd()].push_back(i);
	for(int i=1;i<=n;i++)
		if(!vis[i] && edge[i].size()>=2) dfs(i);
	for(int i=1;i<=n;i++) 
		if(!vis[i]) dfs(i);
	memset(vis,0,sizeof(vis));
	for(int i=1;i<=n;i++)
		if(incir[i] && !vis[i] && ans) ans=mul(ans,calc(i));
	f[0]=f[2]=1;
	for(int i=4;i<=n;i+=2) f[i]=add(f[i-2],mul(mul(i-2,i-3),f[i-4]));
	for(int i=1;i<=n;i++) {
		int sum=0;
		for(int j=0;j<=cnt[i];j+=2) {
			int tp=mul(C.C(cnt[i],j),f[j]);
			if(i!=1) {
				if(i&1) tp=mul(tp,power(2,cnt[i]-j));
				tp=mul(tp,power(i,j/2));
			} sum=add(sum,tp);
		} ans=mul(ans,sum);
	}
	cout<<ans<<'\n'; 
}