神仙题。
排列计数，一种常见的做法是 \(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\ne 0\) 时，先选出这些环，\(\binom{cnt_i}{2j}\)；
- \(j\ne 0\) 时，然后想象成是个二分图，枚举左边的环是哪些，再把右边的环分给左边的环。注意实际上没有顺序，所以要再除掉一堆 \(2\)。\(\frac{\binom{2j}{j}j!}{2^j}\)；
- \(j\ne 0\) 时，每对环有 \(i\) 种拼法，\(i^j\)；
- \(i\) 为奇数且 \(i\ne 1\) 时，没有与别的环拼起来的环可以以两种形态在 \(G_1\) 中出现，\(2^{cnt_i-2j}\)。
再看基环树。基环树之间也相互独立。
挂在基环树上的一堆链要压到环上。图的话，建议继续看学姐的博客。
抓比较重要的几点来说：1、每条链不能重叠，所以压链的话不能超过上一个有链的点。2、如果不考虑 1 中的限制，每条链有两种压法，上面的第一个点在 \(G_1\) 中直接连向 \(u\) 和两步连向 \(u\)。
所以说，每条链的压法只有 \(0,1,2\)，且互相独立。直接乘起来。
一些无解情况也可以很简单判了：在环上的点入度 \(\le 2\)，不在环上的点入度 \(\le 1\)。联想一下基环树的方案数算法应该很好理解。
时间复杂度 \(O(n\log 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\n",ans);
}