神仙题。
排列计数，一种常见的做法是 (i) 向 (p_i) 连边。
然而这里这个就逼迫我们只能从 (i) 向 (a_i) 连边。
不过没关系，考虑从 (i) 向 (p_i) 连边的图（为方便叫 (G_1)）和从 (i) 向 (a_i) 连边的图（为方便叫 (G_2)）的区别。
首先 (G_1) 中每个点入度和出度都是 (1)，所以是一堆环构成的。
考虑一个环：（下面建议画图，懒的建议看 litble 学姐的博客，自己不敢直接把图拿过来）
- 如果上面所有点都满足 (p_i=a_i)，那么这个环在 (G_2) 中也出现了，长得一样。
- 如果上面所有点都满足 (p_{p_i}=a_i)，且这个环长度为奇数，那么有一个相同长度的环在 (G_2) 中出现了，但是长得不一样。
- 如果上面所有点都满足 (p_{p_i}=a_i)，且这个环长度为偶数，那么有两个长度都为这个环长度的一半的环在 (G_2) 中出现了。
- 如果上面的点又有满足 (p_i=a_i) 的，又有满足 (p_{p_i}=a_i) 的，那么会有一个点数相同的基环内向树在 (G_2) 中出现。（这种情况比较复杂，等会再说）
现在知道了 (G_2) ，问 (G_1) 的方案数。
环和基环树独立。先看环。
令长度为 (i) 的环有 (cnt_i) 个。每个环要么是单独在 (G_1) 中出现，要么是与另一个环拼成一个大环在 (G_1) 中出现。
枚举与别的环一起拼的环的个数 (2j)，那么把下面这些全步乘起来：
- (j e 0) 时，先选出这些环，(inom{cnt_i}{2j})；
- (j e 0) 时，然后想象成是个二分图，枚举左边的环是哪些，再把右边的环分给左边的环。注意实际上没有顺序，所以要再除掉一堆 (2)。(frac{inom{2j}{j}j!}{2^j})；
- (j e 0) 时，每对环有 (i) 种拼法，(i^j)；
- (i) 为奇数且 (i e 1) 时，没有与别的环拼起来的环可以以两种形态在 (G_1) 中出现，(2^{cnt_i-2j})。
再看基环树。基环树之间也相互独立。
挂在基环树上的一堆链要压到环上。图的话，建议继续看学姐的博客。
抓比较重要的几点来说：1、每条链不能重叠，所以压链的话不能超过上一个有链的点。2、如果不考虑 1 中的限制，每条链有两种压法，上面的第一个点在 (G_1) 中直接连向 (u) 和两步连向 (u)。
所以说，每条链的压法只有 (0,1,2)，且互相独立。直接乘起来。
一些无解情况也可以很简单判了：在环上的点入度 (le 2)，不在环上的点入度 (le 1)。联想一下基环树的方案数算法应该很好理解。
时间复杂度 (O(nlog n))，如果有闲心也可以做到 (O(n))。
#include<bits/stdc++.h>
using namespace std;
const int maxn=100010,mod=1000000007,inv2=500000004;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline int read(){
int x=0,f=0;char ch=getchar();
while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar();
while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
return f?-x:x;
}
int n,a[maxn],deg[maxn],ans=1,stk[maxn],tp,q[maxn],h,r,len[maxn],seq[maxn],tot,cnt[maxn];
int fac[maxn],inv[maxn],invfac[maxn];
bool vis[maxn],ins[maxn],cyc[maxn];
void dfs(int u){
if(vis[u]){
if(ins[u]){
ROF(i,tp,1){
cyc[stk[i]]=true;
if(stk[i]==u) break;
}
}
return;
}
vis[u]=ins[u]=true;
stk[++tp]=u;
dfs(a[u]);
ins[u]=false;
}
bool dfs2(int u){
if(vis[u]) return true;
vis[u]=true;
seq[++tot]=len[u];
return dfs2(a[u]) && !len[u];
}
inline int C(int n,int m){
return 1ll*fac[n]*invfac[m]%mod*invfac[n-m]%mod;
}
inline int qpow(int a,int b){
int ans=1;
for(;b;b>>=1,a=1ll*a*a%mod) if(b&1) ans=1ll*ans*a%mod;
return ans;
}
int main(){
n=read();
FOR(i,1,n) a[i]=read(),deg[a[i]]++;
FOR(i,1,n) dfs(i);
FOR(i,1,n) if(cyc[i] && deg[i]>=3|| !cyc[i] && deg[i]>=2) return puts("0"),0;
h=1;r=0;
FOR(i,1,n) if(!deg[i]) q[++r]=i;
while(h<=r){
int u=q[h++];
len[a[u]]=len[u]+1;
if(!cyc[a[u]]) q[++r]=a[u];
}
MEM(vis,0);
FOR(i,1,n) if(cyc[i] && !vis[i]){
tot=0;
if(dfs2(i)) cnt[tot]++;
else{
int pre=0;
FOR(j,1,tot) if(seq[j]){
if(pre){
int at=j-seq[j];
if(at<pre) return puts("0"),0;
if(at>pre && tot>=2) ans=2*ans%mod;
}
pre=j;
}
FOR(j,1,tot) if(seq[j]){
int at=j-seq[j]+tot;
if(at<pre) return puts("0"),0;
if(at>pre && tot>=2) ans=2*ans%mod;
break;
}
}
}
fac[0]=fac[1]=inv[1]=invfac[0]=invfac[1]=1;
FOR(i,2,n){
fac[i]=1ll*fac[i-1]*i%mod;
inv[i]=mod-1ll*(mod/i)*inv[mod%i]%mod;
invfac[i]=1ll*invfac[i-1]*inv[i]%mod;
}
FOR(i,1,n) if(cnt[i]){
int s=0;
FOR(j,0,cnt[i]/2){
int x=1ll*C(cnt[i],2*j)*C(2*j,j)%mod*fac[j]%mod;
if(j) x=1ll*x*qpow(inv2,j)%mod*qpow(i,j)%mod;
if(i%2==1 && i!=1) x=1ll*x*qpow(2,cnt[i]-2*j)%mod;
s=(s+x)%mod;
}
ans=1ll*ans*s%mod;
}
printf("%d
",ans);
}