Expected Components
Given a cyclic array $$$a$$$ of size $$$n$$$, where $$$a_i$$$ is the value of $$$a$$$ in the $$$i$$$-th position, there may be repeated values. Let us define that a permutation of $$$a$$$ is equal to another permutation of $$$a$$$ if and only if their values are the same for each position $$$i$$$ or we can transform them to each other by performing some cyclic rotation. Let us define for a cyclic array $$$b$$$ its number of components as the number of connected components in a graph, where the vertices are the positions of $$$b$$$ and we add an edge between each pair of adjacent positions of $$$b$$$ with equal values (note that in a cyclic array the first and last position are also adjacents).Find the expected value of components of a permutation of $$$a$$$ if we select it equiprobably over the set of all the different permutations of $$$a$$$.
In the first test case there is only $$$1$$$ different permutation of $$$a$$$:  $$$[1, 1, 1, 1]$$$ has $$$1$$$ component.  Therefore the expected value of components is $$$\frac{1}{1} = 1$$$ In the second test case there are $$$4$$$ ways to permute the cyclic array $$$a$$$, but there is only $$$1$$$ different permutation of $$$a$$$:  $$$[1, 1, 1, 2]$$$, $$$[1, 1, 2, 1]$$$, $$$[1, 2, 1, 1]$$$ and $$$[2, 1, 1, 1]$$$ are the same permutation and have $$$2$$$ components.  Therefore the expected value of components is $$$\frac{2}{1} = 2$$$ In the third test case there are $$$6$$$ ways to permute the cyclic array $$$a$$$, but there are only $$$2$$$ different permutations of $$$a$$$:  $$$[1, 1, 2, 2]$$$, $$$[2, 1, 1, 2]$$$, $$$[2, 2, 1, 1]$$$ and $$$[1, 2, 2, 1]$$$ are the same permutation and have $$$2$$$ components.  $$$[1, 2, 1, 2]$$$ and $$$[2, 1, 2, 1]$$$ are the same permutation and have $$$4$$$ components.  Therefore the expected value of components is $$$\frac{2+4}{2} = \frac{6}{2} = 3$$$ In the fourth test case there are $$$120$$$ ways to permute the cyclic array $$$a$$$, but there are only $$$24$$$ different permutations of $$$a$$$:  Any permutation of $$$a$$$ has $$$5$$$ components.  Therefore the expected value of components is $$$\frac{24\cdot 5}{24} = \frac{120}{24} = 5$$$ 
The input consists of multiple test cases. The first line contains a single integer $$$t$$$ ($$$1 \le t \le 10^5$$$) — the number of test cases. Description of the test cases follows.The first line of each test case contains a single integer $$$n$$$ ($$$1 \le n \le 10^6$$$) — the size of the cyclic array $$$a$$$.The second line of each test case contains $$$n$$$ integers, the $$$i$$$-th of them is the value $$$a_i$$$ ($$$1 \le a_i \le n$$$).It is guaranteed that the sum of $$$n$$$ over all test cases does not exceed $$$10^6$$$.It is guaranteed that the total number of different permutations of $$$a$$$ is not divisible by $$$M$$$
For each test case print a single integer — the expected value of components of a permutation of $$$a$$$ if we select it equiprobably over the set of all the different permutations of $$$a$$$ modulo $$$998\,244\,353$$$. Formally, let $$$M = 998\,244\,353$$$. It can be shown that the answer can be expressed as an irreducible fraction $$$\frac{p}{q}$$$, where $$$p$$$ and $$$q$$$ are integers and $$$q \not \equiv 0 \pmod{M}$$$. Output the integer equal to $$$p \cdot q^{-1} \bmod M$$$. In other words, output such an integer $$$x$$$ that $$$0 \le x < M$$$ and $$$x \cdot q \equiv p \pmod{M}$$$.
import sys
I = lambda : [*map(int, sys.stdin.readline().split())]
import math
M = 998244353
facs = [1]
ifacs = [1] * 1000001
for i in range(1, 1000001):
    facs.append(facs[-1] * i % M)
ifacs[-1] = pow(facs[-1], M - 2, M)
for i in range(999999, 0, -1):
    ifacs[i] = ifacs[i + 1] * (i + 1) % M
(t,) = I()
for _ in range(t):
    (n,) = I()
    a = I()
    ccount = [0] * n
    for guy in a:
        ccount[guy - 1] += 1
    count = []
    for guy in ccount:
        if guy > 0:
            count.append(guy)
    x = count[0]
    for guy in count:
        x = math.gcd(x, guy)
    if x == n:
        print(1)
        continue
    nperm = 0
    ncomps = 0
    added = {}
    for i in range(1, x + 1):
        if x % i == 0:
            added[i] = 0
    for i in range(1, x + 1):
        if x % i == 0:
            perms = facs[n // i]
            same = 0
            for guy in count:
                y = guy // i
                same += y * (y - 1)
                perms = perms * ifacs[y] % M
            diff = n // i * (n // i - 1) - same
            totcomps = diff * i * pow(n // i - 1, M - 2, M) * perms % M
            need = i - added[i]
            nperm += perms * need
            ncomps += totcomps * need
            for guy in added:
                if guy % i == 0:
                    added[guy] += need
    print(ncomps * pow(nperm, M - 2, M) % M)
(standard_input, packages, output_together) = (1, 1, 1)
(dfs, hashing, read_from_file) = (0, 0, 0)
de = 1
if 1:
    if standard_input:
        import io, os, sys
        input = lambda : sys.stdin.readline().strip()
        inf = float('inf')

        def I():
            return input()

        def II():
            return int(input())

        def MII():
            return map(int, input().split())

        def LI():
            return list(input().split())

        def LII():
            return list(map(int, input().split()))

        def LFI():
            return list(map(float, input().split()))

        def GMI():
            return map(lambda x: int(x) - 1, input().split())

        def LGMI():
            return list(map(lambda x: int(x) - 1, input().split()))
    if packages:
        from io import BytesIO, IOBase
        import math
        import random
        import os
        import bisect
        import typing
        from collections import Counter, defaultdict, deque
        from copy import deepcopy
        from functools import cmp_to_key, lru_cache, reduce
        from heapq import merge, heapify, heappop, heappush, heappushpop, nlargest, nsmallest
        from itertools import accumulate, combinations, permutations, count, product
        from operator import add, iand, ior, itemgetter, mul, xor
        from string import ascii_lowercase, ascii_uppercase, ascii_letters
        from typing import *
        BUFSIZE = 4096
    if output_together:

        class FastIO(IOBase):
            newlines = 0

            def __init__(self, file):
                self._fd = file.fileno()
                self.buffer = BytesIO()
                self.writable = 'x' in file.mode or 'r' not in file.mode
                self.write = self.buffer.write if self.writable else None

            def read(self):
                while True:
                    b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
                    if not b:
                        break
                    ptr = self.buffer.tell()
                    (self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
                self.newlines = 0
                return self.buffer.read()

            def readline(self):
                while self.newlines == 0:
                    b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
                    self.newlines = b.count(b'\n') + (not b)
                    ptr = self.buffer.tell()
                    (self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
                self.newlines -= 1
                return self.buffer.readline()

            def flush(self):
                if self.writable:
                    os.write(self._fd, self.buffer.getvalue())
                    (self.buffer.truncate(0), self.buffer.seek(0))

        class IOWrapper(IOBase):

            def __init__(self, file):
                self.buffer = FastIO(file)
                self.flush = self.buffer.flush
                self.writable = self.buffer.writable
                self.write = lambda s: self.buffer.write(s.encode('ascii'))
                self.read = lambda : self.buffer.read().decode('ascii')
                self.readline = lambda : self.buffer.readline().decode('ascii')
        sys.stdout = IOWrapper(sys.stdout)
    if dfs:
        from types import GeneratorType

        def bootstrap(f, stack=[]):

            def wrappedfunc(*args, **kwargs):
                if stack:
                    return f(*args, **kwargs)
                else:
                    to = f(*args, **kwargs)
                    while True:
                        if type(to) is GeneratorType:
                            stack.append(to)
                            to = next(to)
                        else:
                            stack.pop()
                            if not stack:
                                break
                            to = stack[-1].send(to)
                    return to
            return wrappedfunc
    if hashing:
        RANDOM = random.getrandbits(20)

        class Wrapper(int):

            def __init__(self, x):
                int.__init__(x)

            def __hash__(self):
                return super(Wrapper, self).__hash__() ^ RANDOM
    if read_from_file:
        file = open('input.txt', 'r').readline().strip()[1:-1]
        fin = open(file, 'r')
        input = lambda : fin.readline().strip()
        output_file = open('output.txt', 'w')

        def fprint(*args, **kwargs):
            print(*args, **kwargs, file=output_file)
    if de:

        def debug(*args, **kwargs):
            print('\x1b[92m', end='')
            print(*args, **kwargs)
            print('\x1b[0m', end='')

class Factorial:

    def __init__(self, N, mod) -> None:
        N += 1
        self.mod = mod
        self.f = [1 for _ in range(N)]
        self.g = [1 for _ in range(N)]
        for i in range(1, N):
            self.f[i] = self.f[i - 1] * i % self.mod
        self.g[-1] = pow(self.f[-1], mod - 2, mod)
        for i in range(N - 2, -1, -1):
            self.g[i] = self.g[i + 1] * (i + 1) % self.mod

    def fac(self, n):
        return self.f[n]

    def fac_inv(self, n):
        return self.g[n]

    def combi(self, n, m):
        if n < m or m < 0 or n < 0:
            return 0
        return self.f[n] * self.g[m] % self.mod * self.g[n - m] % self.mod

    def combi_inv(self, n, m):
        return self.g[n] * self.f[m] % self.mod * self.f[n - m] % self.mod

    def permu(self, n, m):
        if n < m or m < 0 or n < 0:
            return 0
        return self.f[n] * self.g[n - m] % self.mod

    def catalan(self, n):
        return (self.combi(2 * n, n) - self.combi(2 * n, n - 1)) % self.mod

    def inv(self, n):
        return self.f[n - 1] * self.g[n] % self.mod
mod = 998244353
fact = Factorial(10 ** 6, mod)

def main():
    n = II()
    nums = LII()
    cnt = list(Counter(nums).values())
    g = math.gcd(*cnt)
    if g == n:
        print(1)
        return
    freq = Counter(cnt)
    counts = Counter()
    ans = Counter()
    factors = []
    for i in range(1, g + 1):
        if g % i == 0:
            factors.append(i)
            tot = fact.fac(n // i)
            expectation = 0
            for x in freq:
                tot = tot * pow(fact.fac_inv(x // i), freq[x], mod) % mod
                (a, b) = ((n - x) // i, x // i)
                expectation += freq[x] * (a + b) * fact.combi(a + b - 2, a - 1) % mod * fact.combi_inv(a + b, a) % mod
            counts[i] = tot
            ans[i] = expectation * tot % mod
    k = len(factors)
    for i in range(k - 1, -1, -1):
        for j in range(i + 1, k):
            if factors[j] % factors[i] == 0:
                counts[factors[i]] -= counts[factors[j]] * (n // factors[j])
                counts[factors[i]] %= mod
                ans[factors[i]] -= ans[factors[j]] * n // factors[i] % mod
                ans[factors[i]] %= mod
        counts[factors[i]] = counts[factors[i]] * fact.inv(n) % mod * factors[i] % mod
        ans[factors[i]] = ans[factors[i]] * fact.inv(n) % mod * factors[i] % mod
    res = 0
    for x in counts:
        res += ans[x] * x % mod
    print(res * pow(sum(counts.values()), -1, mod) % mod)
    return
t = II()
for _ in range(t):
    main()
import math
import sys
r = range
m = lambda x: pow(x, M - 2, M)
M = 998244353
f = [1]
F = [1] * 10 ** 6
for i in r(1, 1000001):
    f.append(f[-1] * i % M)
F.append(m(f[-1]))
for i in r(10 ** 6, 1, -1):
    F[i - 1] = F[i] * i % M
I = sys.stdin.readline
for _ in [0] * int(I()):
    n = int(I())
    a = [*map(int, I().split())]
    C = [0] * n
    (P, c) = (0, 0)
    A = {}
    for g in a:
        C[g - 1] += 1
    C = [g for g in C if g > 0]
    x = C[0]
    for g in C:
        x = math.gcd(x, g)
    if x == n:
        print(1)
        continue
    for i in r(1, n):
        if x % i < 1:
            A[i] = 0
    for i in r(1, n):
        if x % i < 1:
            X = n // i
            p = f[X]
            s = 0
            e = i - A[i]
            for g in C:
                y = g // i
                s += y * (y - 1)
                p = p * F[y] % M
            P += p * e
            c += e * (X * (X - 1) - s) * i * m(X - 1) * p % M
            for g in A:
                if g % i < 1:
                    A[g] += e
    print(c * m(P) % M)