洛谷

(乱搞艹爆正解系列)

    对不起，由于博主太弱了，并不会正解的多项式exp(甚至多项式exp我都不会2333)。

    只能来说一说我是怎么乱搞的啦QWQ

    首先这个题最关键的性质是： 一个在原置换 g 中长度为 l 的环，在 f 中会被拆成 gcd(l,n) 个 长度为 l/gcd(l,n) 的环  [可以类比同余系中的一些变换]

    所以知道了这个，我们怎么反着搞回去呢？？？

    先 O(N) 找出所有循环的长度，然后记 num[i] 为 f 中长度为i的环有多少个。这个时候，我们再知道了 哪些长度一样的环被合并到了一起，就可以计算出答案了。

    因为长度不一样的环之间互不影响，所以我们可以把每个长度的环的答案乘起来，就是最终答案了(如果f中没有长度为i的环，那么这个长度的环的答案是1)。

    于是我们就可以先从1到n扫一遍，把每个 i 存到 i/gcd(i,n) 的 vector 中去，表示 gcd(i,n) 个 i/gcd(i,n) 可以合并成一个 长度为i的环。

    然后对于每个长度，我们只需要扫它的倍数dp即可，复杂度O(玄学)。 

    (至于中间怎么推dp的式子就不说了，就是一些组合数学的基本功练习2333)

  

    为什么是 O(玄学)呢？

    单纯的扫一遍倍数已经是 O（N * log (N)）了，但是因为很多地方可以不用扫(比如没出现的环长，或者一种环长的环的长度和很小，不需要用到很大的dp数组)，仔细算一下，扫的过程其实是 O(图中环的个数) 的，已经是低于线性的了。

    但是可能一个vector里面会有很多元素啊，这样dp的时候会不会炸啊？？？

    就比如 i=1 的时候，我们可能要扫 1~n中所有数的dp数组，再乘一个 vector 的大小，难道不是秒变 O(N^2)了吗？

    但可以发现此时vector中的元素只有 n 的约数(请自行yy为什么)，1e5以下的约数最多的到100吗(不太清楚)？？

   所以就算用心卡我，也很难卡掉吧QWQ

   (达成目标 ： 在代码量是标程 1/3的情况下 速度 比标程快了8倍  2333)

#include<bits/stdc++.h>
#define ll long long
using namespace std;
#define pb push_back
const int maxn=100005,ha=998244353;
inline int add(int x,int y){ x+=y; return x>=ha?x-ha:x;}
inline void ADD(int &x,int y){ x+=y; if(x>=ha) x-=ha;}
inline int ksm(int x,int y){
	int an=1;
	for(;y;y>>=1,x=x*(ll)x%ha) if(y&1) an=an*(ll)x%ha;
	return an;
}

vector<int> g[maxn],w[maxn];
int n,a[maxn],num[maxn],ans=1;
int jc[maxn],ni[maxn],f[maxn];
bool v[maxn];

int gcd(int x,int y){ return y?gcd(y,x%y):x;}

inline int P(int x,int y){ return jc[x]*(ll)ni[x-y]%ha;}

inline int work(int x){
	int an=0;
	while(!v[x]) an++,v[x]=1,x=a[x];
	return an;
}

inline void solve(){
	jc[0]=1;
	for(int i=1;i<=n;i++) jc[i]=jc[i-1]*(ll)i%ha;
	ni[n]=ksm(jc[n],ha-2);
	for(int i=n;i;i--) ni[i-1]=ni[i]*(ll)i%ha;
	
	for(int i=1;i<=n;i++) if(!v[i]) num[work(i)]++;
	
	for(int i=n;i;i--) g[i/gcd(i,n)].pb(i);
	
	for(int i=1;i<=n;i++)
	    for(int j=0;j<g[i].size();j++) w[i].pb(ksm(i,g[i][j]/i-1));
	
	
	for(int i=1;i<=n;i++) if(num[i]){
		int T=num[i]*i,now; 
		f[0]=1;
		
		for(int j=i;j<=T;j+=i){
			f[j]=0;
			for(int l=g[i].size()-1;l>=0;l--){
				now=g[i][l];
				if(now>j) break;
				
				ADD(f[j],f[j-now]*(ll)P(j/i-1,now/i-1)%ha*(ll)w[i][l]%ha);
			}
		}
		
		ans=ans*(ll)f[T]%ha;
//		cout<<f[T]<<endl;
	}
}

int main(){
//	freopen("data.in","r",stdin);
//	freopen("data.out","w",stdout);
	
	scanf("%d",&n);
	for(int i=1;i<=n;i++) scanf("%d",a+i);
	
	solve();
	
	printf("%d
",ans);
	return 0;
}