#include <bits/stdc++.h>

#define I inline
#define fi first
#define se second
#define LL long long
#define mp make_pair
#define reg int
#define pii pair<int,int>
#define fo(i, a, b) for(reg i = a; i <= b; i++)
#define fd(i, a, b) for(reg i = a; i >= b; i--)
#define bp __builtin_popcount
#define vi vector<int>
#define cr const reg&
using namespace std;
const int mod = 998244353, N = 3e6;

I int _max(cr x, cr y) {return x > y ? x : y;}
I int _min(cr x, cr y) {return x < y ? x : y;}
I int read() {
	reg x = 0, f = 1; char ch = getchar();
	while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
	while(ch >= '0' && ch <= '9') x = (x << 3) + (x << 1) + (ch ^ 48), ch = getchar();
	return x * f;
}
I void ptt(LL x) {if(x >= 10) ptt(x / 10); putchar(x % 10 + '0');}
I void put(LL x) {x < 0 ? putchar('-'), ptt(-x) : ptt(x);}
I void pr1(LL x) {put(x), putchar(' ');}
I void pr2(LL x) {put(x), puts("");}

I int pow_mod(reg a, reg k) {reg ans = 1; for(; k; k >>= 1, a = (LL)a * a % mod) if(k & 1) ans = (LL)ans * a % mod; return ans;}
namespace AFastExp
{
	typedef unsigned int uint;
	typedef long long unsigned int uint64;
	constexpr uint Max_size = 3e6;
	constexpr uint g = 3, Mod = 998244353;
	inline uint norm_2(const uint x){return x < Mod * 2 ? x : x - Mod * 2;}
	inline uint norm(const uint x){return x < Mod ? x : x - Mod;}
	struct Z
	{
		uint v;
		Z() { }
		Z(const uint _v) : v(_v) { }
	};
	
	inline Z operator+(const Z x1, const Z x2) { return norm(x1.v + x2.v); }
	inline Z operator-(const Z x1, const Z x2) { return norm(x1.v + Mod - x2.v); }
	inline Z operator-(const Z x) { return x.v ? Mod - x.v : 0; }
	inline Z operator*(const Z x1, const Z x2) { return static_cast<uint64>(x1.v) * x2.v % Mod; }
	inline Z &operator+=(Z &x1, const Z x2) { return x1 = x1 + x2; }
	inline Z &operator-=(Z &x1, const Z x2) { return x1 = x1 - x2; }
	inline Z &operator*=(Z &x1, const Z x2) { return x1 = x1 * x2; }
	inline Z Power(Z Base, int Exp)
	{
		Z res = 1;
		for (; Exp; Base *= Base, Exp >>= 1)if (Exp & 1)res *= Base;
		return res;
	}
	
	inline Z Rec(const Z x)
	{
		return Power(x, Mod - 2);
	}
	
	std::pair<Z, Z> sqrt(const Z x)
	{
		if (x.v == 0)
			return {0, 0};
		Z w = 1;
		while (Power(w * w - x, (Mod - 1) / 2).v == 1)
			++w.v;
		Z Base0 = w, Base1 = 1;
		w = w * w - x;
		Z res0 = 1, res1 = 0;
		int Exp = (Mod + 1) / 2;
		for (; Exp; std::tie(Base0, Base1) = std::make_pair(Base0 * Base0 + Base1 * Base1 * w, 2 * Base0 * Base1), Exp >>= 1)
			if (Exp & 1)
				std::tie(res0, res1) = std::make_pair(res0 * Base0 + res1 * Base1 * w, res0 * Base1 + Base0 * res1);
		Z res = res0.v < Mod - res0.v ? res0 : Mod - res0;
		return {res, Mod - res};
	}
	
	Z _Rec[Max_size];
	
	void init_Rec(const int n)
	{
		_Rec[1] = 1;
		for (int i = 2; i != n; ++i)
			_Rec[i] = _Rec[i - 1] * (i - 1);
		Z R = Rec(_Rec[n - 1] * (n - 1));
		for (int i = n - 1; i != 1; --i)
			_Rec[i] *= R, R *= i;
	}
	
	int size;
	uint w[Max_size], w_q[Max_size];
	
	inline uint mult_Shoup_2(const uint x, const uint y, const uint y_q)
	{
		uint q = static_cast<uint64>(x) * y_q >> 32;
		return x * y - q * Mod;
	}
	
	inline uint mult_Shoup(const uint x, const uint y, const uint y_q)
	{
		return norm(mult_Shoup_2(x, y, y_q));
	}
	
	inline uint mult_Shoup_q(const uint x, const uint y, const uint y_q)
	{
		uint q = static_cast<uint64>(x) * y_q >> 32;
		return q + (x * y - q * Mod >= Mod);
	}
	
	void init_w(const int n)
	{
		for (size = 2; size < n; size <<= 1)
			;
		uint pr = Power(g, (Mod - 1) / size).v;
		uint pr_q = (static_cast<uint64>(pr) << 32) / Mod;
		uint pr_r = (static_cast<uint64>(pr) << 32) % Mod;
		size >>= 1;
		w[size] = 1, w_q[size] = (static_cast<uint64>(w[size]) << 32) / Mod;
		for (int i = 1; i < size; ++i)
		{
			//w[size + i] = mult_Shoup(w[size + i - 1], pr, pr_q);
			uint x = w[size + i - 1];
			uint64 p = static_cast<uint64>(x) * pr_q;
			uint q = p >> 32;
			w[size + i] = norm(x * pr - q * Mod);
			w_q[size + i] = static_cast<uint>(p) + mult_Shoup_q(pr_r, w[size + i - 1], w_q[size + i - 1]);
		}
		for (int i = size - 1; i; --i)
			w[i] = w[i * 2], w_q[i] = w_q[i * 2];
		size <<= 1;
	}
	void DFT_fr_2(Z _A[], const int L)
	{
		if (L == 1)
			return;
		uint *A = reinterpret_cast<uint *>(_A);
	#define butterfly1(a, b)\
		do\
		{\
			uint _a = a, _b = b;\
			uint x = norm_2(_a + _b), y = norm_2(_a + Mod * 2 - _b);\
			a = x, b = y;\
		} while (0)
		if (L == 2)
		{
			butterfly1(A[0], A[1]);
			return;
		}
	#define butterfly(a, b, _w, _w_q)\
		do\
		{\
			uint _a = a, _b = b;\
			uint x = norm_2(_a + _b), y = mult_Shoup_2(_a + Mod * 2 - _b, _w, _w_q);\
			a = x, b = y;\
		} while (0)
		if (L == 4)
		{
			butterfly1(A[0], A[2]);
			butterfly(A[1], A[3], w[3], w_q[3]);
			butterfly1(A[0], A[1]);
			butterfly1(A[2], A[3]);
			return;
		}
		for (int d = L >> 1; d != 4; d >>= 1)
			for (int i = 0; i != L; i += d << 1)
				for (int j = 0; j != d; j += 4)
				{
					butterfly(A[i + j + 0], A[i + d + j + 0], w[d + j + 0], w_q[d + j + 0]);
					butterfly(A[i + j + 1], A[i + d + j + 1], w[d + j + 1], w_q[d + j + 1]);
					butterfly(A[i + j + 2], A[i + d + j + 2], w[d + j + 2], w_q[d + j + 2]);
					butterfly(A[i + j + 3], A[i + d + j + 3], w[d + j + 3], w_q[d + j + 3]);
				}
		for (int i = 0; i != L; i += 8)
		{
			butterfly1(A[i + 0], A[i + 4]);
			butterfly(A[i + 1], A[i + 5], w[5], w_q[5]);
			butterfly(A[i + 2], A[i + 6], w[6], w_q[6]);
			butterfly(A[i + 3], A[i + 7], w[7], w_q[7]);
		}
		for (int i = 0; i != L; i += 8)
		{
			butterfly1(A[i + 0], A[i + 2]);
			butterfly(A[i + 1], A[i + 3], w[3], w_q[3]);
			butterfly1(A[i + 4], A[i + 6]);
			butterfly(A[i + 5], A[i + 7], w[3], w_q[3]);
		}
		for (int i = 0; i != L; i += 8)
		{
			butterfly1(A[i + 0], A[i + 1]);
			butterfly1(A[i + 2], A[i + 3]);
			butterfly1(A[i + 4], A[i + 5]);
			butterfly1(A[i + 6], A[i + 7]);
		}
	#undef butterfly1
	#undef butterfly
	}
	
	void DFT_fr(Z _A[], const int L)
	{
		DFT_fr_2(_A, L);
		for (int i = 0; i != L; ++i)
			_A[i] = norm(_A[i].v);
	}
	void IDFT_fr(Z _A[], const int L)
	{
		if (L == 1)
			return;
		uint *A = reinterpret_cast<uint *>(_A);
	#define butterfly1(a, b)\
		do\
		{\
			uint _a = a, _b = b;\
			uint x = norm_2(_a), t = norm_2(_b);\
			a = x + t, b = x + Mod * 2 - t;\
		} while (0)
		if (L == 2)
		{
			butterfly1(A[0], A[1]);
			A[0] = norm(norm_2(A[0])), A[0] = A[0] & 1 ? A[0] + Mod : A[0], A[0] /= 2;
			A[1] = norm(norm_2(A[1])), A[1] = A[1] & 1 ? A[1] + Mod : A[1], A[1] /= 2;
			return;
		}
	#define butterfly(a, b, _w, _w_q)\
		do\
		{\
			uint _a = a, _b = b;\
			uint x = norm_2(_a), t = mult_Shoup_2(_b, _w, _w_q);\
			a = x + t, b = x + Mod * 2 - t;\
		} while (0)
		if (L == 4)
		{
			butterfly1(A[0], A[1]);
			butterfly1(A[2], A[3]);
			butterfly1(A[0], A[2]);
			butterfly(A[1], A[3], w[3], w_q[3]);
			std::swap(A[1], A[3]);
			for (int i = 0; i != L; ++i)
			{
				uint64 m = -A[i] & 3;
				A[i] = norm((A[i] + m * Mod) >> 2);
			}
			return;
		}
		for (int i = 0; i != L; i += 8)
		{
			butterfly1(A[i + 0], A[i + 1]);
			butterfly1(A[i + 2], A[i + 3]);
			butterfly1(A[i + 4], A[i + 5]);
			butterfly1(A[i + 6], A[i + 7]);
		}
		for (int i = 0; i != L; i += 8)
		{
			butterfly1(A[i + 0], A[i + 2]);
			butterfly(A[i + 1], A[i + 3], w[3], w_q[3]);
			butterfly1(A[i + 4], A[i + 6]);
			butterfly(A[i + 5], A[i + 7], w[3], w_q[3]);
		}
		for (int i = 0; i != L; i += 8)
		{
			butterfly1(A[i + 0], A[i + 4]);
			butterfly(A[i + 1], A[i + 5], w[5], w_q[5]);
			butterfly(A[i + 2], A[i + 6], w[6], w_q[6]);
			butterfly(A[i + 3], A[i + 7], w[7], w_q[7]);
		}
		for (int d = 8; d != L; d <<= 1)
			for (int i = 0; i != L; i += d << 1)
				for (int j = 0; j != d; j += 4)
				{
					butterfly(A[i + j + 0], A[i + d + j + 0], w[d + j + 0], w_q[d + j + 0]);
					butterfly(A[i + j + 1], A[i + d + j + 1], w[d + j + 1], w_q[d + j + 1]);
					butterfly(A[i + j + 2], A[i + d + j + 2], w[d + j + 2], w_q[d + j + 2]);
					butterfly(A[i + j + 3], A[i + d + j + 3], w[d + j + 3], w_q[d + j + 3]);
				}
	#undef butterfly1
	#undef butterfly
		std::reverse(A + 1, A + L);
		int k = __builtin_ctz(L);
		for (int i = 0; i != L; ++i)
		{
			uint64 m = -A[i] & (L - 1);
			A[i] = norm((A[i] + m * Mod) >> k);
		}
	}
	
	inline void reset(Z A[], const int n)
	{
		std::memset(reinterpret_cast<uint *>(A), 0, n * sizeof(Z));
	}
	
	const uint _M_size = Max_size * 5;
	Z _M[_M_size], *_M_pos = _M + _M_size;
	
	Z *Exp_2_res;
	int Exp_2_N, Exp_2_D;
	int Exp_2_K[5], Exp_2_Q[5];
	Z *Exp_2_dftc[5];
	Z *Exp_2_DC;
	uint64 Exp_2_tmp[Max_size];
	
	void Exp_2_rec(const int L, const int R, const int d)
	{
		Z *const res = Exp_2_res;
		if (d + 1 == Exp_2_D)
		{
			for (int i = L + !L; i != R && i != Exp_2_N; ++i)
			{
				uint64 r = 0;
				int j = 1;
				if (i - L >= 4)
				{
					uint64 r0 = 0, r1 = 0, r2 = 0, r3 = 0;
					for (; j + 3 <= i - L; j += 4)
					{
						r0 += static_cast<uint64>(res[i - j - 0].v) * Exp_2_DC[j + 0].v;
						r1 += static_cast<uint64>(res[i - j - 1].v) * Exp_2_DC[j + 1].v;
						r2 += static_cast<uint64>(res[i - j - 2].v) * Exp_2_DC[j + 2].v;
						r3 += static_cast<uint64>(res[i - j - 3].v) * Exp_2_DC[j + 3].v;
					}
					if (((i - L) & ~3) <= 16)
						r = (r0 + r1 + r2 + r3) % Mod;
					else if (((i - L) & ~3) >> 1 <= 16)
						r = (r0 + r1) % Mod + (r2 + r3) % Mod;
					else
						r = r0 % Mod + r1 % Mod + r2 % Mod + r3 % Mod;
				}
				for (; j <= i - L; ++j)
					r += static_cast<uint64>(res[i - j].v) * Exp_2_DC[j].v;
				res[i] = (res[i].v + r) % Mod;
				res[i] *= _Rec[i];
			}
			return;
		}
		int Q = Exp_2_Q[d], bK = Exp_2_K[d + 1], bS = 1 << bK;
		std::memset(Exp_2_tmp + L * 2, 0, sizeof(uint64) * (2 << Exp_2_K[d]));
		while (L + ((Q - 1) << bK) >= Exp_2_N)
			--Q;
		for (int q = 0; q != Q; ++q)
		{
			int l = L + (q << bK);
			if (l >= Exp_2_N)
				return;
			if (q)
			{
				Z *dftr = _M_pos -= bS * 2;
				for (int i = 0; i != bS * 2; ++i)
					dftr[i] = Exp_2_tmp[l * 2 + i] % Mod;
				IDFT_fr(dftr, bS * 2);
				for (int i = 0; i != bS; ++i)
					res[l + i] += dftr[bS + i];
				_M_pos += bS * 2;
			}
			Exp_2_rec(l, l + bS, d + 1);
			if (q != Q - 1)
			{
				Z *dftr = _M_pos -= bS * 2;
				std::copy_n(res + l, bS, dftr), reset(dftr + bS, bS);
				DFT_fr(dftr, bS * 2);
				for (int i = q + 1; i != Q; ++i)
				{
					Z *a = Exp_2_dftc[d] + ((i - 1 - q) << bK) * 2;
					uint64 *r = Exp_2_tmp + (L + (i << bK)) * 2;
					for (int j = 0; j != 2 << bK; j += 4)
					{
						r[j + 0] += static_cast<uint64>(dftr[j + 0].v) * a[j + 0].v;
						if (r[j + 0] >= 0ULL - (Mod - 1ULL) * (Mod - 1ULL))
							r[j + 0] %= Mod;
						r[j + 1] += static_cast<uint64>(dftr[j + 1].v) * a[j + 1].v;
						if (r[j + 1] >= 0ULL - (Mod - 1ULL) * (Mod - 1ULL))
							r[j + 1] %= Mod;
						r[j + 2] += static_cast<uint64>(dftr[j + 2].v) * a[j + 2].v;
						if (r[j + 2] >= 0ULL - (Mod - 1ULL) * (Mod - 1ULL))
							r[j + 2] %= Mod;
						r[j + 3] += static_cast<uint64>(dftr[j + 3].v) * a[j + 3].v;
						if (r[j + 3] >= 0ULL - (Mod - 1ULL) * (Mod - 1ULL))
							r[j + 3] %= Mod;
					}
				}
				_M_pos += bS * 2;
			}
		}
	}
	
	void Exp_2(const Z C[], Z res[], const int N, const int L)
	{
		if (L == 1)
			return;
		Exp_2_N = N;
		int (&K)[5] = Exp_2_K, (&Q)[5] = Exp_2_Q;
		Z *(&dftc)[5] = Exp_2_dftc;
		int D = 0;
		for (K[0] = __builtin_ctz(L); K[D]; ++D)
		{
			int k = static_cast<int>(std::log2(K[D]) + 2 + 0.5);
			if (1 << K[D] <= 64 || k >= K[D])
				k = K[D];
			K[D + 1] = K[D] - k;
			Q[D] = 1 << k;
		}
		Exp_2_D = D;
		for (int d = 0; d != D - 1; ++d)
		{
			dftc[d] = _M_pos -= 2 << K[d];
			for (int q = 0; q != Q[d] - 1; ++q)
			{
				for (int i = 0; i != 2 << K[d + 1]; ++i)
					dftc[d][(q << K[d + 1]) * 2 + i] = C[(q << K[d + 1]) + i] * ((q << K[d + 1]) + i);
				DFT_fr(dftc[d] + (q << K[d + 1]) * 2, 2 << K[d + 1]);
			}
		}
		Exp_2_DC = _M_pos -= 1 << K[D - 1];
		for (int i = 0; i != 1 << K[D - 1]; ++i)
			Exp_2_DC[i] = C[i] * i;
		Exp_2_res = res;
		reset(res, N), res[0] = 1;
		Exp_2_rec(0, L, 0);
		_M_pos += 1 << K[D - 1];
		for (int d = 0; d != D - 1; ++d)
			_M_pos += 2 << K[d];
	}
	
	int N, L;
	Z A[Max_size];
	
	void gtexp(int *a,int *res,int n)
	{
		N=n;
		for(int i=0;i<n;++i)A[i].v=a[i];
		for (L = 1; L < N; L <<= 1);
		init_w(L << 1);
		init_Rec(L << 1);
		
		Exp_2(A, A, N, L);
		for(int i=0;i<n;++i)res[i]=A[i].v;
	}
}//namespace AFastExp by duyi
using AFastExp::gtexp; 

int jc[N], inv[N];
int p[N], S[N], L[N], f[N], z[N]; bool v[N];

int main() {
//	freopen("a.in","r",stdin);
//	freopen("a.out","w",stdout);
	
		reg n = read();
		fo(i, 1, n) p[i] = read(), S[i] = L[i] = v[i] = 0;
		reg mn = n;
		fo(i, 1, n) if(!v[i]){
			reg s = 0;
			reg x = i; while(!v[x]) v[x] = 1, x = p[x], s++;
			S[s]++, mn = _min(mn, s);
		} reg ans;
		jc[0] = 1; fo(i, 1, n) jc[i] = (LL)jc[i - 1] * i % mod;
		inv[n] = pow_mod(jc[n], mod - 2); fd(i, n, 1) inv[i - 1] = (LL)inv[i] * i % mod;
		if(mn > 2) {puts("0"); return 0;}
		if(S[mn] == 1) ans = 1; else ans = pow_mod(S[mn], S[mn] - 2);
		if(mn == 2) ans = (LL)ans * S[mn] % mod * pow_mod(2, S[mn] - 1) % mod;
		fo(i, 1, n) {
			for(reg j = i + i; j <= n; j += i) L[j] = (L[j] + (LL)S[i] * i) % mod;
		} 
		fo(i, mn + 1, n) if(S[i]){
			f[0] = 1; f[1] = L[i];
			reg t = L[i];
			fo(j, 2, S[i]) f[j] = (LL)L[i] * pow_mod(j, j - 2) % mod * inv[j] % mod * pow_mod(i, j - 1) % mod * j % mod;
			gtexp(f, z, S[i] + 1);
			ans = (LL)ans * jc[S[i]] % mod * z[S[i]] % mod;
		} pr2(ans);
	return 0;
}





Measured memory_bytes: 17092