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 <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <string>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>
#include <sstream>
#include <stack>
#include <iomanip>
using namespace std;
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
typedef double ld;
typedef vector<int> vi;
#define fi first
#define se second
#define fe first
#define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);}
#define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);}
#define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);}
#define es(x,e) (int e=fst[x];e;e=nxt[e])
#define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e])
#define SZ 666666
int n,a[SZ],ff[SZ],sz[SZ]; Edg
int gf(int x) {return ff[x]?ff[x]=gf(ff[x]):x;}
void uni(int a,int b)
{
	a=gf(a),b=gf(b);
	(a^b)?(ff[a]=b):0;
}
map<int,int> cp;
bool ic[SZ],vv[SZ];
int dfs(int x)
{
	int ans=0,chh=0; vv[x]=1;
	for esb(x,e,b) if(ic[b]||vv[b]);else
	{
		int t=dfs(b);++chh;
		if(t<0) return -1;
		ans=max(ans,t+1);
	}
	if(chh>=2) return -1;
	return ans;
}
int py(int a,int b)
{
	if(b>a) return 0;
	if(a==b) return 1;
	return 2;
}
const ll MOD=1e9+7;
ll f[SZ];
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;++i)
		scanf("%d",a+i),adde(i,a[i]),uni(i,a[i]);
	for(int i=1;i<=n;++i) ++sz[gf(i)];
	ll ans=1;
	for(int i=1;i<=n;++i) if(!ff[i])
	{
		int g=i;
		for(int j=1;j<=sz[i];++j) g=a[g];
		vector<int> cyc; cyc.pb(g);
		for(int h=a[g];h!=g;h=a[h]) cyc.pb(h);
		if(cyc.size()==sz[i])
		{++cp[sz[i]]; continue;}
		for(auto r:cyc) ic[r]=1;
		vector<int> tmp;
		for(auto r:cyc)
		{
			int t=dfs(r);
			if(t<0) {puts("0"); return 0;}
			tmp.pb(t);
		}
		int cl=0,pv=-1,bl=0;
		for(int i=0;i<tmp.size();++i)
		{
			++cl;
			if(!tmp[i]) continue;
			if(pv==-1) bl=cl,pv=tmp[i];
			else ans*=py(cl,tmp[i]),ans%=MOD;
			cl=0;
		}
		ans*=py(bl+cl,pv); ans%=MOD;
	}
	for(auto g:cp)
	{
		int x=g.fi,y=g.se;
		f[0]=1;
		for(int i=1;i<=y;++i)
			f[i]=(((i>=2)?(f[i-2]*(i-1)%MOD*x):0)+(1+(x%2==1&&x>1))*f[i-1])%MOD;
		ans*=f[y]; ans%=MOD;
	}
	cout<<ans<<"\n";
}