We introduce two solutions:

1. primitive \(\mathrm{O}(U \log^2 U)\) solution using Pólya enumeration theorem
2. advanced \(\mathrm{O}(U \log U)\) solution using Dirichlet generating functions

The latter solution using DGF involves magical techniques. It is a truly astonishing algorithm, so do read the editorial.

1. Solution using Pólya enumeration theorem

We already introduced the Pólya theorem in the editorial of ABC284-Ex, so we omit the details here.

First, let us solve the sequence case. Define \(A(n)\) as the number of length-\(n\) sequences whose product is \(i\).

By the constraints, the length is always not greater than \(\lfloor \log_2 U \rfloor\). Let \(V = \lfloor \log_2 U \rfloor\).
Then \(A(0), A(1), \dots, A(V)\) can be easily computed with the following DP in \(\mathrm{O}(U \log^2 U)\) time.

using mint = atcoder::modint998244353;
int freq[500003]; // Precomputed frequencies of b_i
mint A[21][500003];
for(int i = 1; i <= U; i++) A[0][i] = (i == 1 ? 1 : 0);
for(int i = 1; i <= V; i++)
  for(int s = 1; s <= U; s++)
    for(int t = 1; s * t <= U; t++)
      A[i + 1][s * t] += A[i][s] * freq[t];

Now consider a ring case using the obtained \(A(n)\).
Let us compute the number of necklaces of length \(L\) and product \(x\) using Pólya enumeration theorem. By taking the averaging of the numbers of fixed points, we obtain

\[ \begin{aligned} &\frac{1}{L} \sum_{i=0}^{L-1} A\left( \frac{L}{\gcd(i, L)}\right) [\sqrt[\gcd(i,L)]{x}] \\ &= \frac{1}{L} \sum_{d \vert L} \varphi(d) A \left( \frac{L}{d} \right) [\sqrt[d]{x}]. \end{aligned} \]

If \(\sqrt[d]{x}\) does not take an integer value, let the contribution of the term containing \(\sqrt[d]{x}\) be \(0\).

Just as the sequence case, \(L\) is capped by \(V\). Thus, the number of necklaces with product \(x\) can be found as the sum of the expression above for \(L=1,2,\dots,V\), which can be done in a total of \(\mathrm{O}(V \log V)\) time if totient functions are precomputed.

Therefore, the problem can be solved in \(\mathrm{O}(U \log^2 U)\) time, which is fast enough.

2. Solution using Dirichlet generating functions (DGF)

DGF is a generating function (GF) defined for \(a_1, a_2, \dots\) defined as follows:

\[A(s) = \frac{a_1}{1^s} + \frac{a_2}{2^s} + \dots = \sum_{n=1}^\infty \frac{a_n}{n^s}\]

DGFs are useful for problems like ours, where products are important.

When you put a jewel of beauty \(2\) and another of beauty \(3\), you obtain a jewel sequence with total product \(6\). This process corresponds to a product of DGFs:

\[\frac{1}{2^s} \times \frac{1}{3^s} = \frac{1}{6^s}\]

Now let us try to describe the problem using DGFs. Let \(A(s)\) be the DGF of the frequencies of jewel beauties.

By Pólya enumeration theorem, the GF of the necklaces of length \(L\) can be represented as:

\[C_L(s) = \frac{1}{L} \sum_{D \vert L} \varphi(d) A(ds)^{L/d}\]

So, the GF for the necklaces satisfies:

\[ \begin{aligned} C(s) &= \sum_L C_L(s) \\ &= \sum_L \frac{1}{L} \sum_{D \vert L} \varphi(d) A(ds)^{L/d} \\ &= \sum_{D \vert L}\frac{1}{L} \varphi(d) A(ds)^{L/d} \\ &= \sum_{D} \sum_{t \geq 1} \frac{1}{dt} \varphi(d) A(ds)^{t} \\ &= \sum_{D} \frac{\varphi(d)}{d}\sum_{t \geq 1} \frac{A(ds)^{t} }{t} \\ &= \sum_{D} \frac{\varphi(d)}{d} \log \left( \frac{1}{1 - A(ds)}\right). \end{aligned} \]

All that left is to evaluate this expression. Here, terms like \(\log \left( \frac{1}{1 - A(2s)}\right)\) and\(\log \left( \frac{1}{1 - A(3s)}\right)\) can be obtained by replacing terms like \(\frac{c}{n^s}\) in \(\log \left( \frac{1}{1 - A(s)}\right)\) with \(\frac{c}{(n^2)^s}\) or \(\frac{c}{(n^3)^s}\), so we can compute them all if we know \(\log \left( \frac{1}{1 - A(s)}\right)\). Also, it is obvious that it is sufficient to compute within \(1 \leq D \leq \lfloor \log_2 U \rfloor\).
Thus, it turns out that the problem can be solved if we can compute \(\log \left( \frac{1}{1 - A(s)}\right)\).

It can be easily confirmed that convolutions and inverse of DGFs up to the \(U\)-th term can be computed in \(\mathrm{O}(U \log U)\) (if you are not sure, think about it), so the only issue is how to compute \(\log\).
One possible approach follows from the fact that, in its Taylor series, the \(\left( \lfloor \log_2 U \rfloor + 1 \right)\)-th and further terms are all \(0\), This allows us to compute the first \(\lfloor \log_2 U \rfloor\) terms using convolutions, which gives us an\(\mathrm{O}(U \log^2 U)\) solution.
To optimize it more, let us use derivatives. Given a function \(G(s)\), consider how to find \(F(s) = \log G(s)\).

\[ \begin{aligned} &\log G = F \\ &\to e^F = G \\ &\to F' e^F = F' G = G'\\ &\to F' = \frac{G'}{G}, \end{aligned} \]

so convolutions and inverse operations on DGFs seem sufficient to obtain \(F\) from \(G\).

However we have an issue: the derivative of the DGF \(A(s)\) of \(a_1, a_2, \) is

\[A'(s) = \sum_{n =1}^\infty \frac{a_n \ln n}{n^s},\]

which contains a natural logarithm \(\ln n\). Of course, \(\ln n\) cannot be represented modulo \(998244353\), begin a large obstacle for computing the derivatives.

Although we are stuck, here comes the magic. According to apiad’s blog, this was discovered by _rqy and Elegia. The magic is \(\Omega(n)\), the function that returns the number of prime factors in \(n\), counted with multiplicity. Roughly speaking, it gives us how many multiplications of primes we need starting from \(1\). At a glance, it does not seem anything special, but in fact, we can replace \(\Omega(n)\) with \(\ln n\)!

We will explain formally. Define a special derivative operator \(\mathfrak{D}\) as:

\[\mathfrak{D} A(s) = \sum_{n =1}^\infty \frac{a_n \Omega(n)}{n^s}\]

Then surprisingly, we can prove that any DGFs \(F\) and \(G\) and an ordinary function \(H\) satisfies

\[\mathfrak{D} (FG) = F \mathfrak{D} G + G \mathfrak{D} F,\]

\[\mathfrak{D}(h \circ G) = (h' \circ G) \mathfrak{D} G.\]

In other words, this magic derivative can be treated in the same way as an ordinary derivatives!

Hence, \(F = \log(G)\) can be obtained as \(F\) satisfying

\[\mathfrak{D} F = \frac{ \mathfrak{D} G}{G},\]

which can be computed in \(\mathrm{O}(U \log U)\) time. As you can see, although we cannot define derivatives on DGFs (modulo \(998244353\)), we can tweak the definition of derivatives to get everything work. This concludes the explanation of the magical technique.
By implementing all these algorithms, the problem can be computed in \(\mathrm{O}(U \log U)\) time, which is very fast.

The following is sample code (Python, codon). Since the language updates have took place, I tried writing it in Codon. Even with a rough one-shot coding, the code runs in about 120 ms. It’s nice to see modint works fast enough in Python, isn’t it?
As of 2025, the mainstream Python runtime is PyPy, but it has a numerous fields that it is not good at (e.g. recursions, sorting tuples, data structures, modint, abstraction…). They have dropped shadow on the Python language in algo (algorithm contests in competitive programming) and burdened those who only write Python, as many contestants reading this editorial probably know. Codon claims to get rid of all these issues, and I strongly desire that it becomes a Tier-1 runtime unlike PyPy. (I want a template that absorbs subtle difference from Python, especially map behavior. Also I want ACL (AtCoder Library) for Codon…)

@tuple
class modint:
    x: u32

    def __new__():
        return modint(u32(0))

    def __new__(y: int) -> modint:
        z = y % 998244353
        if z < 0:
            z += 998244353
        return modint(u32(z))

    def __add__(self, rhs: modint) -> modint:
        res = self.x + rhs.x
        if res >= u32(998244353):
            res -= u32(998244353)
        return modint(res)

    def __sub__(self, rhs: modint) -> modint:
        res = self.x - rhs.x
        if res >= u32(998244353):
            res += u32(998244353)
        return modint(res)

    def __mul__(self, rhs: modint) -> modint:
        res = u64(self.x) * u64(rhs.x) % u64(998244353)
        return modint(u32(res))
    
    def __neg__(self) -> modint:
        return modint(u32(0) if self.x == u32(0) else u32(998244353) - self.x)

    def __str__(self) -> str:
        return str(self.x)

def dirichlet_conv(a, b):
    c = [modint(0)] * (U + 1)
    for i in range(1, U + 1):
        for j in range(1, U // i + 1):
            c[i * j] += a[i] * b[j]
    return c

def dirichlet_inv(a):
    assert a[1] == modint(1)
    b = [modint(0)] * (U + 1)
    b[1] = modint(1)
    for i in range(1, U + 1):
        for j in range(2, U // i + 1):
            b[i * j] -= b[i] * a[j]
    return b

N, U = [int(a) for a in input().split()]

# Omega(n)
lpf = [i for i in range(U + 1)]
for p in range(2, U + 1):
    if lpf[p] != p:
        continue
    for j in range(p * p, U + 1, p):
        lpf[j] = min(lpf[j], p)
Omega = [0] * (U + 1)
for i in range(2, U + 1):
    Omega[i] = Omega[i // lpf[i]] + 1

# totient
totient = [i for i in range(30)]
for i in range(1, 30):
    for j in range(2 * i, 30, i):
        totient[j] -= totient[i]

# modinv
inv = [modint(0)] + [modint(pow(i, 998244353 - 2, 998244353)) for i in range(1, 30)]

A = [modint(0)] * (U + 1)
for b in map(int, input().split()):
    A[b] += modint(1)

# 1 - A(s)
B = [modint(0), modint(1)] + [-x for x in A[2:]]
# 1 / (1 - A(s)) = inv(B)
G = dirichlet_inv(B)
# log(G) = F <-> F' = G' B
for i in range(1, U + 1):
    G[i] *= modint(Omega[i])
F = dirichlet_conv(G, B)
for i in range(2, U + 1):
    F[i] *= inv[Omega[i]]

ans = [modint(0)] * (U + 1)
for i in range(2, U + 1):
    j, d = i, 1
    while j <= U:
        ans[j] += F[i] * modint(totient[d]) * inv[d]
        j, d = j * i, d + 1
print(" ".join(map(str, ans[2:])))