When analyzing subgroups of the symmetric group generated by a set of permutations, a crucial observation is the uniform distribution property across reachable configurations. Consider a subgroup $H \leq S_n$ generated by permutations $\pi_1, \ldots, \pi_k$. For any fixed index $i$, if the $i$-th position takes values $v_1, \ldots, v_m$ across elements of $H$, each value appears equally often. This extends to tuples of positions: for indices $i_1, \ldots, i_t$, the distribution of $(h(i_1), \ldots, h(i_t))$ is uniform over all attainable tuples in $H$.

To establish this, model the group action as a graph. For tuples of length $m$, vertices represent $m$-tuples of distinct indices. Directed edges correspond to generator applications: from $(a_1, \ldots, a_m)$ to $(\pi(a_1), \ldots, \pi(a_m))$ labeled with $\pi$. Since $\pi^{-1} \in H$, edges are effectively bidirectional. Fix a root tuple $r$ and consider a spanning tree. Each vertex $v$ corresponds to a unique path $r \to v$, yielding a group element $g_v$. The mapping $g \mapsto g \circ g_v^{-1}$ establishes a bijection between paths to $v$ and paths to $r$, proving uniformity.

This property enables efficient computation of average inversion counts. For the problem of computing expected inversions over the generated subgroup, construct a graph on ordered pairs $(i,j)$ with $i \neq j$. For each generator $\pi$, add undirected edges between $(i,j)$ and $(\pi(i), \pi(j))$. Each connected component $C$ contains some number $a$ of pairs with $i < j$ and $b$ with $i > j$. By uniformity, the probability of an inversion for a random element of $H$ restricted to pairs in $C$ is $\frac{b}{a+b}$, contributing $\frac{ab}{a+b}$ to the total expectation.

Direct computation uses $O(n^3)$ time via union-find on $O(n^2)$ pair nodes. However, the generator set can be reduced randomly. Sampling $O(\log n)$ random subsets of the original generators, composing them arbitrarily, produces a small generating set with high probability. This yields an $O(n^2 \log n)$ algorithm:

const int MAXN = 505;
const int MOD = 998244353;

int dim, genCnt;
int generators[MAXN][MAXN];
int perm[MAXN], temp[MAXN];
int parent[MAXN * MAXN], rev[MAXN * MAXN];
int lessCnt[MAXN * MAXN], greatCnt[MAXN * MAXN];

int encode(int x, int y) { return (x - 1) * dim + y; }

int find(int x) {
    return parent[x] == x ? x : parent[x] = find(parent[x]);
}

int main() {
    dim = read(); genCnt = read();
    for (int g = 1; g <= genCnt; ++g)
        for (int i = 1; i <= dim; ++i)
            generators[g][i] = read();
    
    for (int i = 1; i <= dim * dim; ++i) parent[i] = i;
    
    for (int trial = 0; trial < LOG; ++trial) {
        for (int i = 1; i <= dim; ++i) perm[i] = i;
        for (int g = 1; g <= genCnt; ++g) {
            if (rng() & 1) {
                for (int i = 1; i <= dim; ++i) temp[perm[i]] = generators[g][i];
                for (int i = 1; i <= dim; ++i) perm[i] = temp[i];
            }
        }
        for (int i = 1; i <= dim; ++i) {
            for (int j = 1; j <= dim; ++j) {
                if (i == j) continue;
                parent[find(encode(i, j))] = find(encode(perm[i], perm[j]));
            }
        }
    }
    
    for (int i = 1; i <= dim; ++i)
        for (int j = 1; j <= dim; ++j)
            if (i != j) {
                if (i < j) ++lessCnt[find(encode(i, j))];
                else ++greatCnt[find(encode(i, j))];
            }
    
    rev[1] = 1;
    for (int i = 2; i <= dim * dim; ++i)
        rev[i] = (long long)rev[MOD % i] * (MOD - MOD / i) % MOD;
    
    int ans = 0;
    for (int i = 1; i <= dim * dim; ++i)
        if (parent[i] == i && lessCnt[i] && greatCnt[i])
            ans = (ans + (long long)lessCnt[i] * greatCnt[i] % MOD * rev[lessCnt[i] + greatCnt[i]]) % MOD;
    
    printf("%d\n", ans);
    return 0;
}

For computing the subgroup order (the size of $H$), apply a recursive dimension reduction. The number of distinct values at position $n$ equals the size of the orbit of $n$, computable via the graph method above. To count stabilizers (elements fixing $n$), construct new generators from cycles involving $n$: for each edge $(u,v)$ in the orbit graph, the permutation mapping $n \to u \to v \to n$ stabilizes $n$ when appropriately conjugated.

Specifically, build a BFS tree from node $n$ in the action graph on single points. For each tree edge to node $u$ via generator $g$, and for each outgoing edge from $u$ to $v$ via generator $h$, the element $\text{path}[u] \cdot g \cdot h \cdot \text{path}[v]^{-1}$ fixes $n$, where $\text{path}[x]$ denotes the tree path from $n$ to $x$. This generates the stabilizer subgroup $H_n$. Randomly reducing these $O(n \cdot \text{genCnt})$ elements to $O(n)$ generators allows recursion on $n-1$ dimnesions.

The recursive algorithm operates in $O(n^5)$ time:

struct Permutation {
    int img[MAXN];
    void identity() { for (int i = 1; i <= curSize; ++i) img[i] = i; }
    friend Permutation compose(const Permutation& a, const Permutation& b) {
        Permutation c;
        for (int i = 1; i <= curSize; ++i) c.img[i] = b.img[a.img[i]];
        return c;
    }
    Permutation invert() {
        Permutation inv;
        for (int i = 1; i <= curSize; ++i) inv.img[img[i]] = i;
        return inv;
    }
};

Permutation gens[MAXN], reduced[MAXN];
Permutation treePath[MAXN << 1];
vector<pair<int,int>> adj[MAXN];
bool visited[MAXN];
int curSize, genCount;
long long result[20];
int resLen;
long long orbitSize;

void dfs(int u, Permutation cur) {
    treePath[u] = cur;
    visited[u] = true;
    ++orbitSize;
    for (auto [v, gid] : adj[u]) {
        if (visited[v]) continue;
        dfs(v, compose(cur, gens[gid]));
    }
}

void reduceGenerators() {
    for (int i = 1; i <= curSize; ++i) adj[i].clear();
    for (int g = 1; g <= genCount; ++g)
        for (int i = 1; i <= curSize; ++i)
            adj[i].push_back({gens[g].img[i], g});
    
    Permutation id;
    id.identity();
    dfs(curSize, id);
    
    for (int i = 1; i <= curSize; ++i) {
        treePath[i + curSize] = treePath[i].invert();
        adj[i].clear();
    }
    
    for (int t = 1; t <= BOUND; ++t) {
        reduced[t].identity();
        for (int g = 1; g <= genCount; ++g)
            for (int i = 1; i <= curSize; ++i)
                if (visited[i] && (rng() & 1))
                    reduced[t] = compose(reduced[t], compose(compose(treePath[i], gens[g]), treePath[gens[g].img[i] + curSize]));
    }
    
    for (int i = 1; i <= curSize; ++i) visited[i] = false;
    genCount = BOUND;
    for (int i = 1; i <= genCount; ++i) gens[i] = reduced[i];
    
    --curSize;
    long long carry = 0;
    for (int i = 0; i <= resLen; ++i) {
        result[i] = result[i] * orbitSize + carry;
        carry = result[i] / BASE;
        result[i] %= BASE;
    }
    if (carry) result[++resLen] = carry;
    orbitSize = 0;
}

int main() {
    curSize = read(); genCount = read();
    for (int i = 1; i <= genCount; ++i)
        for (int j = 1; j <= curSize; ++j)
            gens[i].img[j] = read();
    
    result[resLen = 0] = 1;
    while (curSize > 0) reduceGenerators();
    
    printf("%lld", result[resLen]);
    for (int i = resLen - 1; i >= 0; --i)
        printf("%016lld", result[i]);
    puts("");
    return 0;
}

The rank of a finite group—the minimum size of a generating set—can reach $\Theta(n)$, necessitating retention of $O(n)$ generators for order computation. However, for connectivity problems like inversion counting, $O(\log n)$ random generators suffice to connect the graph with high probability, explaining the efficacy of the randomized reduction.