MOD = 10**9 + 7 def main(): import sys input = sys.stdin.read().split() idx = 0 N = int(input[idx]) idx += 1 P = list(map(int, input[idx:idx+N])) idx += N # Precompute factorials up to N max_n = N fact = [1] * (max_n + 1) for i in range(2, max_n + 1): fact[i] = fact[i-1] * i % MOD # Compute cycle decomposition visited = [False] * (N + 1) cycle_counts = {} for i in range(1, N+1): if not visited[i]: current = i cycle_length = 0 while not visited[current]: visited[current] = True current = P[current - 1] # P is 0-based in list, but elements are 1-based cycle_length += 1 if cycle_length in cycle_counts: cycle_counts[cycle_length] += 1 else: cycle_counts[cycle_length] = 1 # Check if there exists m with (m even and a_m >=1) or (m odd and a_m >=2) has_odd_centralizer = False for m in cycle_counts: a_m = cycle_counts[m] if (m % 2 == 0 and a_m >= 1) or (m % 2 == 1 and a_m >= 2): has_odd_centralizer = True break if has_odd_centralizer: # Compute denominator denominator = 1 for m in cycle_counts: a_m = cycle_counts[m] denominator = denominator * pow(m, a_m, MOD) % MOD denominator = denominator * fact[a_m] % MOD size = fact[N] * pow(denominator, MOD-2, MOD) % MOD print(size) else: if N == 1: print(1) else: # Compute denominator denominator = 1 for m in cycle_counts: a_m = cycle_counts[m] denominator = denominator * pow(m, a_m, MOD) % MOD denominator = denominator * fact[a_m] % MOD size = fact[N] * pow(denominator, MOD-2, MOD) % MOD inv_2 = 500000004 # Modular inverse of 2 mod 1e9+7 ans = size * inv_2 % MOD print(ans) if __name__ == "__main__": main()