Commutativity
Source: aizu

Commutativity
Problem Statement

You are given a function $F: \{1, 2, ..., N\} \rightarrow \{1, 2, ..., N\}$. In other words, $F$ is a function that takes an integer $x$ between $1$ and $N$ inclusive and returns an integer $F(x)$ between $1$ and $N$ inclusive. Your task is to count the number of functions $G:\{1, 2, ..., N\} \rightarrow \{1, 2, ..., N\}$ such that $F$ and $G$ commute under composition: that is, $F(G(x))=G(F(x))$ holds for any $x \in \{1, 2, ..., N\}$.

As this number could be large, print the answer modulo $998,244,353$.

Input

The input consists of a single test case of the following format.

$N$
$F_1$ $F_2$ ... $F_N$


The first line consists of an integer $N$ between $1$ and $5,000$, inclusive. The second line consists of $N$ positive integers $F_1, F_2, ..., F_N$. For each $i$ ($1 \leq i \leq N$), $F_i$ represents the value of $F(i)$. It is guaranteed that $1 \leq F_i \leq N$.

Output

Print the number of possible functions $G: \{1, 2, ..., N\} \rightarrow \{1, 2, ..., N\}$ modulo $998,244,353$.

Sample Input 1
5
4 5 3 2 1

Sample Output 1
5

Sample Input 2
8
2 3 1 3 6 7 5 5

Sample Output 2
64

Sample Input 3
10
3 1 4 1 5 9 2 6 5 3

Sample Output 3
64

Sample Input 4
15
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Sample Output 4
567381138


C++ solution:
// AOJ #3639 Commutativity
// 2025.3.16

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = 998244353;

int N;
vector<int> f;
vector<vector<int>> inv;

vector<int> col;
vector<bool> cyc;
vector<int> cl;
vector<vector<int>> cycList;

vector<vector<int>> ch;

void dfsCycle(int v){
    vector<int> st;
    int cur = v;
    while(true){
        if(col[cur]==0){
            col[cur]=1;
            st.push_back(cur);
            cur = f[cur];
        } else if(col[cur]==1){
            int idx=-1;
            for (int i=0; i<st.size(); i++){
                if(st[i]==cur){ idx = i; break; }
            }
            int L = st.size()-idx;
            vector<int> cycV;
            for (int i=idx; i<st.size(); i++){
                int u = st[i];
                cyc[u]=true; cl[u] = L;
                cycV.push_back(u);
            }
            cycList.push_back(cycV);
            for (int u : st) col[u]=2;
            break;
        } else {
            for (int u : st) col[u]=2;
            break;
        }
    }
}

vector<vector<int>> dp;
vector<bool> dpDone;

void dfsDP(int v){
    for (int w : ch[v]) dfsDP(w);
    for (int x=1; x<=N; x++){
        ll sum = 0;
        for (int u : inv[x]){
            ll prod = 1;
            for (int w : ch[v]) prod = (prod * dp[w][u]) % MOD;
            sum = (sum + prod) % MOD;
        }
        dp[v][x] = (int)sum;
    }
    dpDone[v] = true;
}

int Fval(int v, int x){
    ll prod = 1;
    for (int w : ch[v]) prod = (prod * dp[w][x]) % MOD;
    return (int) prod;
}

int main(){
    ios::sync_with_stdio(0);
    cin.tie(0);

    cin >> N;
    f.resize(N+1);
    for (int i=1; i<=N; i++) cin >> f[i];

    inv.assign(N+1, vector<int>());
    for (int i=1; i<=N; i++) inv[f[i]].push_back(i);

    col.assign(N+1, 0);
    cyc.assign(N+1, false);
    cl.assign(N+1, 0);
    for (int i=1; i<=N; i++){
        if(col[i]==0) dfsCycle(i);
    }

    ch.assign(N+1, vector<int>());
    for (int i=1; i<=N; i++){
        if(!cyc[i]){
            int p = f[i];
            ch[p].push_back(i);
        }
    }

    dp.assign(N+1, vector<int>(N+1,0));
    dpDone.assign(N+1, false);
    for (int i=1; i<=N; i++){
        if(!cyc[i] && !dpDone[i]) dfsDP(i);
    }

    ll ans = 1;
    for (auto &cycV : cycList){
        int start = cycV[0];
        vector<int> dC;
        int cur = start;
        do {
            dC.push_back(cur);
            cur = f[cur];
        } while(cur != start);
        int L = dC.size();

        ll sumC = 0;
        for (int a = 1; a<=N; a++){
            if(!cyc[a]) continue;
            if(L % cl[a] != 0) continue;
            ll prod = 1;
            int t = a;
            for (int i=0; i<L; i++){
                prod = (prod * Fval(dC[i], t)) % MOD;
                t = f[t];
            }
            sumC = (sumC + prod) % MOD;
        }
        ans = (ans * sumC) % MOD;
    }
    cout << ans << endl;
    return 0;
}

Go solution:
package main

import (
	"bufio"
	"fmt"
	"os"
	"strconv"
	"strings"
)

const MOD = 998244353

func main() {
	scanner := bufio.NewScanner(os.Stdin)
	scanner.Scan()
	N, _ := strconv.Atoi(scanner.Text())

	scanner.Scan()
	line := scanner.Text()
	parts := strings.Fields(line)
	f := make([]int, N+1)
	for i := 1; i <= N; i++ {
		f[i], _ = strconv.Atoi(parts[i-1])
	}

	inv := make([][]int, N+1)
	for i := 0; i <= N; i++ {
		inv[i] = make([]int, 0)
	}
	for i := 1; i <= N; i++ {
		inv[f[i]] = append(inv[f[i]], i)
	}

	col := make([]int, N+1)
	cyc := make([]bool, N+1)
	cl := make([]int, N+1)
	var cycList [][]int

	var dfsCycle func(int)
	dfsCycle = func(v int) {
		stack := make([]int, 0)
		cur := v
		for {
			if col[cur] == 0 {
				col[cur] = 1
				stack = append(stack, cur)
				cur = f[cur]
			} else if col[cur] == 1 {
				idx := -1
				for i := 0; i < len(stack); i++ {
					if stack[i] == cur {
						idx = i
						break
					}
				}
				L := len(stack) - idx
				cycV := make([]int, 0)
				for i := idx; i < len(stack); i++ {
					u := stack[i]
					cyc[u] = true
					cl[u] = L
					cycV = append(cycV, u)
				}
				cycList = append(cycList, cycV)
				for _, u := range stack {
					col[u] = 2
				}
				break
			} else {
				for _, u := range stack {
					col[u] = 2
				}
				break
			}
		}
	}

	for i := 1; i <= N; i++ {
		if col[i] == 0 {
			dfsCycle(i)
		}
	}

	ch := make([][]int, N+1)
	for i := 0; i <= N; i++ {
		ch[i] = make([]int, 0)
	}
	for i := 1; i <= N; i++ {
		if !cyc[i] {
			p := f[i]
			ch[p] = append(ch[p], i)
		}
	}

	dp := make([][]int, N+1)
	for i := 0; i <= N; i++ {
		dp[i] = make([]int, N+1)
	}
	dpDone := make([]bool, N+1)

	var dfsDP func(int)
	dfsDP = func(v int) {
		for _, w := range ch[v] {
			dfsDP(w)
		}
		for x := 1; x <= N; x++ {
			sum := 0
			for _, u := range inv[x] {
				prod := 1
				for _, w := range ch[v] {
					prod = (prod * dp[w][u]) % MOD
				}
				sum = (sum + prod) % MOD
			}
			dp[v][x] = sum
		}
		dpDone[v] = true
	}

	for i := 1; i <= N; i++ {
		if !cyc[i] && !dpDone[i] {
			dfsDP(i)
		}
	}

	Fval := func(v, x int) int {
		prod := 1
		for _, w := range ch[v] {
			prod = (prod * dp[w][x]) % MOD
		}
		return prod
	}

	ans := 1
	for _, cycV := range cycList {
		start := cycV[0]
		dC := make([]int, 0)
		cur := start
		for {
			dC = append(dC, cur)
			cur = f[cur]
			if cur == start {
				break
			}
		}
		L := len(dC)

		sumC := 0
		for a := 1; a <= N; a++ {
			if !cyc[a] {
				continue
			}
			if L%cl[a] != 0 {
				continue
			}
			prod := 1
			t := a
			for i := 0; i < L; i++ {
				prod = (prod * Fval(dC[i], t)) % MOD
				t = f[t]
			}
			sumC = (sumC + prod) % MOD
		}
		ans = (ans * sumC) % MOD
	}
	fmt.Println(ans)
}

Java solution:
import java.io.*;
import java.util.*;

public class Main {
    static final int MOD = 998244353;

    public static void main(String[] args) throws IOException {
        BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
        int N = Integer.parseInt(br.readLine());
        int[] f = new int[N + 1];
        StringTokenizer tokenizer = new StringTokenizer(br.readLine());
        for (int i = 1; i <= N; i++) {
            f[i] = Integer.parseInt(tokenizer.nextToken());
        }

        List<List<Integer>> inv = new ArrayList<>();
        for (int i = 0; i <= N; i++) {
            inv.add(new ArrayList<>());
        }
        for (int i = 1; i <= N; i++) {
            inv.get(f[i]).add(i);
        }

        int[] col = new int[N + 1];
        boolean[] cyc = new boolean[N + 1];
        int[] cl = new int[N + 1];
        List<List<Integer>> cycList = new ArrayList<>();

        for (int i = 1; i <= N; i++) {
            if (col[i] == 0) {
                List<Integer> stack = new ArrayList<>();
                int cur = i;
                while (true) {
                    if (col[cur] == 0) {
                        col[cur] = 1;
                        stack.add(cur);
                        cur = f[cur];
                    } else if (col[cur] == 1) {
                        int idx = -1;
                        for (int j = 0; j < stack.size(); j++) {
                            if (stack.get(j) == cur) {
                                idx = j;
                                break;
                            }
                        }
                        int L = stack.size() - idx;
                        List<Integer> cycle = new ArrayList<>();
                        for (int j = idx; j < stack.size(); j++) {
                            int u = stack.get(j);
                            cyc[u] = true;
                            cl[u] = L;
                            cycle.add(u);
                        }
                        cycList.add(cycle);
                        for (int u : stack) {
                            col[u] = 2;
                        }
                        break;
                    } else {
                        for (int u : stack) {
                            col[u] = 2;
                        }
                        break;
                    }
                }
            }
        }

        List<List<Integer>> ch = new ArrayList<>();
        for (int i = 0; i <= N; i++) {
            ch.add(new ArrayList<>());
        }
        for (int i = 1; i <= N; i++) {
            if (!cyc[i]) {
                int p = f[i];
                ch.get(p).add(i);
            }
        }

        int[][] dp = new int[N + 1][N + 1];
        boolean[] dpDone = new boolean[N + 1];

        for (int i = 1; i <= N; i++) {
            if (!cyc[i] && !dpDone[i]) {
                dfsDP(i, dp, dpDone, ch, inv, N);
            }
        }

        long ans = 1;
        for (List<Integer> cycle : cycList) {
            int start = cycle.get(0);
            List<Integer> cycleNodes = new ArrayList<>();
            int cur = start;
            do {
                cycleNodes.add(cur);
                cur = f[cur];
            } while (cur != start);
            int L = cycleNodes.size();

            long sumC = 0;
            for (int a = 1; a <= N; a++) {
                if (!cyc[a]) continue;
                if (L % cl[a] != 0) continue;
                long prod = 1;
                int t = a;
                for (int i = 0; i < L; i++) {
                    int node = cycleNodes.get(i);
                    prod = (prod * Fval(node, t, dp, ch)) % MOD;
                    t = f[t];
                }
                sumC = (sumC + prod) % MOD;
            }
            ans = (ans * sumC) % MOD;
        }
        System.out.println(ans);
    }

    static void dfsDP(int v, int[][] dp, boolean[] dpDone, List<List<Integer>> ch, List<List<Integer>> inv, int N) {
        for (int w : ch.get(v)) {
            dfsDP(w, dp, dpDone, ch, inv, N);
        }
        for (int x = 1; x <= N; x++) {
            long sum = 0;
            for (int u : inv.get(x)) {
                long prod = 1;
                for (int w : ch.get(v)) {
                    prod = (prod * dp[w][u]) % MOD;
                }
                sum = (sum + prod) % MOD;
            }
            dp[v][x] = (int) sum;
        }
        dpDone[v] = true;
    }

    static int Fval(int v, int x, int[][] dp, List<List<Integer>> ch) {
        long prod = 1;
        for (int w : ch.get(v)) {
            prod = (prod * dp[w][x]) % MOD;
        }
        return (int) prod;
    }
}

JavaScript solution:
const MOD = 998244353;

function main() {
    const fs = require('fs');
    const data = fs.readFileSync(0, 'utf-8').split(/\s+/);
    let ptr = 0;
    const N = parseInt(data[ptr++]);
    const F = new Array(N + 1).fill(0);
    for (let i = 1; i <= N; i++) {
        F[i] = parseInt(data[ptr++]);
    }

    // Build inverse map
    const inv = new Array(N + 1).fill().map(() => []);
    for (let i = 1; i <= N; i++) {
        inv[F[i]].push(i);
    }

    // Detect cycles
    const col = new Array(N + 1).fill(0); // 0: unvisited, 1: visiting, 2: visited
    const cyc = new Array(N + 1).fill(false);
    const cl = new Array(N + 1).fill(0); // cycle length for each node
    const cycList = [];

    function dfsCycle(v) {
        const st = [];
        let cur = v;
        while (true) {
            if (col[cur] === 0) {
                col[cur] = 1;
                st.push(cur);
                cur = F[cur];
            } else if (col[cur] === 1) {
                let idx = -1;
                for (let i = 0; i < st.length; i++) {
                    if (st[i] === cur) {
                        idx = i;
                        break;
                    }
                }
                const L = st.length - idx;
                const cycV = [];
                for (let i = idx; i < st.length; i++) {
                    const u = st[i];
                    cyc[u] = true;
                    cl[u] = L;
                    cycV.push(u);
                }
                cycList.push(cycV);
                for (const u of st) col[u] = 2;
                break;
            } else {
                for (const u of st) col[u] = 2;
                break;
            }
        }
    }

    for (let i = 1; i <= N; i++) {
        if (col[i] === 0) {
            dfsCycle(i);
        }
    }

    // Build children (tree parts)
    const ch = new Array(N + 1).fill().map(() => []);
    for (let i = 1; i <= N; i++) {
        if (!cyc[i]) {
            const p = F[i];
            ch[p].push(i);
        }
    }

    // DP for tree parts
    const dp = new Array(N + 1).fill().map(() => new Array(N + 1).fill(0));
    const dpDone = new Array(N + 1).fill(false);

    function dfsDP(v) {
        for (const w of ch[v]) {
            dfsDP(w);
        }
        for (let x = 1; x <= N; x++) {
            let sum = 0;
            for (const u of inv[x]) {
                let prod = 1;
                for (const w of ch[v]) {
                    prod = (prod * dp[w][u]) % MOD;
                }
                sum = (sum + prod) % MOD;
            }
            dp[v][x] = sum;
        }
        dpDone[v] = true;
    }

    for (let i = 1; i <= N; i++) {
        if (!cyc[i] && !dpDone[i]) {
            dfsDP(i);
        }
    }

    function Fval(v, x) {
        let prod = 1;
        for (const w of ch[v]) {
            prod = (prod * dp[w][x]) % MOD;
        }
        return prod;
    }

    let ans = 1;
    for (const cycV of cycList) {
        const start = cycV[0];
        const dC = [];
        let cur = start;
        do {
            dC.push(cur);
            cur = F[cur];
        } while (cur !== start);
        const L = dC.length;

        let sumC = 0;
        for (let a = 1; a <= N; a++) {
            if (!cyc[a]) continue;
            if (L % cl[a] !== 0) continue;
            let prod = 1;
            let t = a;
            for (let i = 0; i < L; i++) {
                prod = (prod * Fval(dC[i], t)) % MOD;
                t = F[t];
            }
            sumC = (sumC + prod) % MOD;
        }
        ans = (ans * sumC) % MOD;
    }
    console.log(ans);
}

main();

Python solution:
import sys
from collections import defaultdict

MOD = 998244353

def main():
    input = sys.stdin.read().split()
    ptr = 0
    N = int(input[ptr])
    ptr += 1
    F = list(map(int, input[ptr:ptr+N]))
    ptr += N
    F = [0] + F  # 1-based indexing

    # Build inverse mapping
    inv = defaultdict(list)
    for i in range(1, N+1):
        inv[F[i]].append(i)

    # Detect cycles
    col = [0] * (N + 1)
    cyc = [False] * (N + 1)
    cl = [0] * (N + 1)
    cycList = []

    for i in range(1, N+1):
        if col[i] == 0:
            stack = []
            cur = i
            while True:
                if col[cur] == 0:
                    col[cur] = 1
                    stack.append(cur)
                    cur = F[cur]
                elif col[cur] == 1:
                    idx = -1
                    for j in range(len(stack)):
                        if stack[j] == cur:
                            idx = j
                            break
                    L = len(stack) - idx
                    cycV = []
                    for j in range(idx, len(stack)):
                        u = stack[j]
                        cyc[u] = True
                        cl[u] = L
                        cycV.append(u)
                    cycList.append(cycV)
                    for u in stack:
                        col[u] = 2
                    break
                else:
                    for u in stack:
                        col[u] = 2
                    break

    # Build children for tree parts
    ch = [[] for _ in range(N+1)]
    for i in range(1, N+1):
        if not cyc[i]:
            p = F[i]
            ch[p].append(i)

    # DP tables
    dp = [[0] * (N + 1) for _ in range(N + 1)]
    dpDone = [False] * (N + 1)

    def dfsDP(v):
        for w in ch[v]:
            dfsDP(w)
        for x in range(1, N + 1):
            total = 0
            for u in inv[x]:
                product = 1
                for w in ch[v]:
                    product = (product * dp[w][u]) % MOD
                total = (total + product) % MOD
            dp[v][x] = total
        dpDone[v] = True

    for i in range(1, N + 1):
        if not cyc[i] and not dpDone[i]:
            dfsDP(i)

    def Fval(v, x):
        product = 1
        for w in ch[v]:
            product = (product * dp[w][x]) % MOD
        return product

    ans = 1
    for cycV in cycList:
        start = cycV[0]
        dC = []
        cur = start
        while True:
            dC.append(cur)
            cur = F[cur]
            if cur == start:
                break
        L = len(dC)

        sumC = 0
        for a in range(1, N + 1):
            if not cyc[a]:
                continue
            if L % cl[a] != 0:
                continue
            product = 1
            t = a
            for i in range(L):
                product = (product * Fval(dC[i], t)) % MOD
                t = F[t]
            sumC = (sumC + product) % MOD
        ans = (ans * sumC) % MOD

    print(ans)

if __name__ == '__main__':
    main()

Ruby solution:
require 'set'

MOD = 998244353

def main
  n = gets.to_i
  f = gets.split.map(&:to_i)
  f.unshift(0) # Make it 1-based

  # Build inverse mapping
  inv = Array.new(n + 1) { [] }
  (1..n).each { |i| inv[f[i]] << i }

  # Initialize cycle detection variables
  color = Array.new(n + 1, 0)
  in_cycle = Array.new(n + 1, false)
  cycle_length = Array.new(n + 1, 0)
  cycles = []

  # Detect cycles
  (1..n).each do |i|
    next if color[i] != 0

    path = []
    current = i
    while true
      if color[current] == 0
        color[current] = 1
        path << current
        current = f[current]
      elsif color[current] == 1
        # Found a cycle
        cycle_start = path.index(current)
        cycle_nodes = path[cycle_start..-1]
        len = cycle_nodes.size
        cycle_nodes.each do |node|
          in_cycle[node] = true
          cycle_length[node] = len
        end
        cycles << cycle_nodes
        path.each { |node| color[node] = 2 }
        break
      else
        path.each { |node| color[node] = 2 }
        break
      end
    end
  end

  # Build tree structure (children)
  children = Array.new(n + 1) { [] }
  (1..n).each do |i|
    next if in_cycle[i]
    children[f[i]] << i
  end

  # Initialize DP tables
  dp = Array.new(n + 1) { Array.new(n + 1, 0) }
  dp_done = Array.new(n + 1, false)

  # Process nodes in reverse order for DP
  (1..n).reverse_each do |i|
    next if in_cycle[i] || dp_done[i]

    stack = [[i, false]]
    until stack.empty?
      node, processed = stack.pop
      if processed
        (1..n).each do |x|
          sum = 0
          inv[x].each do |u|
            product = 1
            children[node].each { |child| product = (product * dp[child][u]) % MOD }
            sum = (sum + product) % MOD
          end
          dp[node][x] = sum
        end
        dp_done[node] = true
      else
        stack << [node, true]
        children[node].reverse_each do |child|
          stack << [child, false] unless dp_done[child]
        end
      end
    end
  end

  result = 1
  cycles.each do |cycle_nodes|
    cycle = []
    current = cycle_nodes.first
    loop do
      cycle << current
      current = f[current]
      break if current == cycle_nodes.first
    end
    len = cycle.size

    cycle_sum = 0
    (1..n).each do |a|
      next unless in_cycle[a]
      next unless len % cycle_length[a] == 0

      product = 1
      t = a
      cycle.each do |node|
        temp_product = 1
        children[node].each { |child| temp_product = (temp_product * dp[child][t]) % MOD }
        product = (product * temp_product) % MOD
        t = f[t]
      end
      cycle_sum = (cycle_sum + product) % MOD
    end
    result = (result * cycle_sum) % MOD
  end

  puts result
end

main