#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<vector>
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  return p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++;
}
inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

const int N=100005;
const int P=1e9+7;

struct edge{
  int u,v,next;
}G[N<<1]; int head[N],inum=1;
inline void add(int u,int v,int p){
  G[p].u=u; G[p].v=v; G[p].next=head[u]; head[u]=p;
}

int n,a[N];
int deg[N];

namespace Calc{
  const int INV2=(P+1)/2;
  ll fac[N],inv2[N],pow[N],inv[N];
  inline void Pre(int n){
    fac[0]=1; for (int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;
    inv[1]=1; for (int i=2;i<=n;i++) inv[i]=inv[P%i]*(P-P/i)%P;
    inv[0]=1; for (int i=1;i<=n;i++) (inv[i]*=inv[i-1])%=P;
    inv2[0]=1; for (int i=1;i<=n;i++) inv2[i]=inv2[i-1]*INV2%P;
    pow[0]=1; for (int i=1;i<=n;i++) pow[i]=pow[i-1]*2%P;
  }
  inline ll C(int n,int m){
    return fac[n]*inv[m]%P*inv[n-m]%P;
  }
  inline ll Pow(ll a,int b){
    ll ret=1;
    for (;b;a=a*a%P,b>>=1)
      if (b&1)
    (ret*=a)%=P;
    return ret;
  }
  int cnt[N];
  inline ll Solve(){
    Pre(n);
    ll ans=1;
    for (int i=1;i<=n;i++){
      if (!cnt[i]) continue;
      ll ret=0,tmp;
      for (int a=0;2*a<=cnt[i];a++){
    tmp=C(cnt[i],2*a)*(fac[2*a]*inv2[a]%P)%P*inv[a]%P*Pow(i,a)%P;
    if ((i&1) && (i>1)) (tmp*=pow[cnt[i]-2*a])%=P;
    ret=(ret+tmp)%P;
      }
      (ans*=ret)%=P;
    }
    return ans;
  }
}

int cnt;
int depth[N],fat[N];
int path[N],len,_l[N];
int inp[N];
#define V G[p].v
inline void dfs(int u,int fa){
  cnt++;
  for (int p=head[u];p;p=G[p].next)
    if (p!=(fa^1)){
      if (!depth[V])
    fat[V]=u,depth[V]=depth[u]+1,dfs(V,p);
      else if (depth[V]<=depth[u]){
    int t=u;
    while (t!=V) path[++len]=t,t=fat[t]; path[++len]=t;
      }
    }
}

int maxd;

inline void dfs2(int u,int fa){
  cnt++;
  depth[u]=depth[fa]+1; maxd=max(maxd,depth[u]);
  for (int p=head[u];p;p=G[p].next)
    if (V!=fa && !inp[V])
      dfs2(V,u);
}

int l[N],d[N];

int main(){
  ll Ans=1;
  read(n);
  for (int i=1;i<=n;i++) read(a[i]),deg[a[i]]++,add(i,a[i],++inum),add(a[i],i,++inum);
  for (int i=1;i<=n;i++) if (deg[i]>2) return printf("0\n"),0;
  for (int i=1;i<=n;i++){
    if (depth[i]) continue;
    cnt=len=0; depth[i]=1; dfs(i,0);
    if (len==2 && path[1]==path[2]) len--;
    if (cnt==len || cnt==1){
      Calc::cnt[::cnt]++; continue;
    }
    for (int j=1;j<=len;j++) inp[path[j]]=1;
    for (int j=1;j<=len;j++){
      cnt=maxd=0; dfs2(path[j],0);
      if (cnt>maxd) return printf("0\n"),0;
      _l[j]=maxd;
    }
    if (len>1 && a[path[1]]!=path[2]){
      reverse(path+1,path+len+1);
      reverse(_l+1,_l+len+1);
    }
    ll tmp=1; *l=*d=0;
    for (int j=1;j<=len;j++) if (_l[j]>1) l[++*l]=_l[j]-1,d[++*d]=j;
    l[++*l]=l[1]; d[++*d]=d[1]+len;
    for (int j=2;j<=*l && tmp;j++)
      if (d[j]-d[j-1]>=l[j]){
    if (l[j]<d[j]-d[j-1])
      (tmp*=2)%=P;
      }else
    tmp=0;
    (Ans*=tmp)%=P;
  }
  (Ans*=Calc::Solve())%=P;
  printf("%lld\n",Ans);
  return 0;
}

Measured memory_bytes: 12032