ヒント → https://atcoder.jp/contests/arc192/editorial/12151

\(p_i\) を \(i\) 番目に小さい素数とします。

また、\(d_{y}(x)\) を、\(x\) が \(y\) で割り切れる最大の回数とします。

各 \(i\) について、\(d_{p_i}\left(f\left(\dfrac{x}{y}\right)\right)=\left|d_{p_i}(x)-d_{p_i}(y)\right|\) となることがわかるので、素因数ごとに独立に処理できることがわかります。素因数ごとの答えをかけ合わせれば答えが求まるので、素因数ごとに処理することを考えます。

素因数 \(p_i\) に対する答えを考えます。\(s_i=(d_{p_i}(S_1),d_{p_i}(S_2),\dots,d_{p_i}(S_N))\) とすると、\(s_i\) としてあり得るものは次の条件をどちらも満たすものです。

• \(\min(s_i)=0\)
• \(1\leq j\lt N\) なる各 \(j\) について \(\left|s_{i,j}-s_{i,j+1}\right|=d_{p_i}\left(A_i\right)\)

このような \(s_i\) すべてについての \(p_i^{\text{sum}(s_i)}\) の総和を求める問題を解けばよくて、これは現時点での \(s_{i,j}\) と \(\min(s_i)=0\) かのフラグをもってDPすることで実現できます。遷移は実装例を参照してください。

計算量を見積もりましょう。\(\sum_{i=1}^{\infty}\sum_{j=1}^{N-1}d_{p_i}(A_j)\) は高々 \(O(N\log\max(A))\) なので、DPテーブルの大きさは \(O(N^2\log\max(A))\) となります。遷移は \(O(1)\) 時間なので、素因数分解するのも含めると計算量は \(O(N\sqrt{\max(A)}+N^2\log\max(A))\) となり、十分高速です。素因数分解に \(O(N\sqrt{\max(A)})\) ではなく\(O(N\max(A))\) かけても通るはずです。

実装例 (PyPy3)

MOD = 998244353
N = int(input())
A = list(map(int, input().split()))
ans = 1

for i in range(2, max(A) + 1):
    a = [0] * (N - 1)
    for j in range(N - 1):
        while A[j] % i == 0:
            a[j] += 1
            A[j] //= i

    sa = sum(a)

    if sa == 0:
        continue

    pows = [1] * (sa + 1)

    for j in range(1, sa + 1):
        pows[j] = (pows[j - 1] * i) % MOD

    dp = [[pows[j], 0] for j in range(sa + 1)]
    dp[0] = [0, 1]

    for j in range(N - 1):
        ndp = [[0, 0] for _ in range(sa + 1)]
        for k in range(sa + 1):
            if k - a[j] >= 0:
                nxt = 0
                if k - a[j] == 0:
                    nxt = 1

                ndp[k - a[j]][nxt] += dp[k][0] * pows[k - a[j]]
                ndp[k - a[j]][nxt] %= MOD

                ndp[k - a[j]][1] += dp[k][1] * pows[k - a[j]]
                ndp[k - a[j]][1] %= MOD

            if a[j] != 0 and k + a[j] <= sa:
                ndp[k + a[j]][0] += dp[k][0] * pows[k + a[j]]
                ndp[k + a[j]][0] %= MOD

                ndp[k + a[j]][1] += dp[k][1] * pows[k + a[j]]
                ndp[k + a[j]][1] %= MOD

        dp = ndp

    nsum = 0
    for j in range(sa + 1):
        nsum = (nsum + dp[j][1]) % MOD

    ans = (ans * nsum) % MOD

print(ans)