[0] Definitions

• Let \(f(x) = \displaystyle \prod_{i=0}^{N-1} (Bx + Di + C)\).
• Let \(G = \gcd(f(0), f(1), \ldots, f(N))\).
• For a positive integer \(n\) and a prime \(p\), let \(v_p(n)\) denote the number of times \(n\) is divisible by \(p\).
• For integers \(a\) and \(b\) with \(a \neq 0\), write \(a \mid b\) if \(b\) is a multiple of \(a\), and \(a \nmid b\) otherwise.

---

[1] Reduction to \(\gcd(B, C, D) = 1\)

Let \(g = \gcd(B, C, D)\). If \(g \neq 1\), then \(f(x) = g^N \displaystyle \prod_{i=0}^{N-1} \left(\frac{B}{g}x + \frac{D}{g}i + \frac{C}{g}\right)\), so the answer is \(g^N\) times the answer to the problem with \(B \to \dfrac{B}{g}, C \to \dfrac{C}{g}, D \to \dfrac{D}{g}\).

In the following discussion, we assume \(g = 1\).

---

[2] \(v_p(G)\) for each prime \(p\)

We consider \(v_p(G)\) for each prime \(p\). We state the conclusion first.

1. If \(p \mid B\) and \(p \mid D\): \(v_p(G) = 0\)
2. If \(p \nmid B\) and \(p \mid D\): \(v_p(G) = 0\)
3. If \(p \nmid B\) and \(p \nmid D\): \(v_p(G) = v_p(N!)\)
4. If \(p \mid B\) and \(p \nmid D\): \(v_p(G)\) can be computed efficiently by calculation

---

[2-1] When \(p \mid B\) and \(p \mid D\)

By assumption, \(p \nmid C\). Thus, \(p \nmid (Bx + Di + C)\) for any \(x, i\), so \(v_p(G) = 0\).

---

[2-2] When \(p \nmid B\) and \(p \mid D\)

We have \(p \nmid (Bn + C)\) for at least one of \(n = 0, 1\). For such \(n\), \(p \nmid f(n)\), so \(v_p(G) = 0\).

---

[2-3] When \(p \nmid B\) and \(p \nmid D\)

We prove that \(v_p(G) = v_p(N!)\).

[2-3-1] Proof of \(v_p(G) \leq v_p(N!)\)

We define the following operators on functions \(f(x)\):

• Identity operator \(I\): \(If(x) = f(x)\)
• Shift operator \(E\): \(Ef(x) = f(x+1)\)
• Difference operator \(\Delta\): \(\Delta f(x) = f(x+1) - f(x)\ (= (E-I)f(x))\)

[2-3-1-1] Proof of \(\displaystyle \Delta^N f(x) = \sum_{k=0}^N \binom{N}{k}(-1)^{N-k}f(x+k)\)

By the binomial theorem, \(\displaystyle \Delta^N = (E-I)^N = \sum_{k=0}^N \binom{N}{k} E^k (-I)^{N-k}\).

Since \(E^k f(x) = f(x+k)\) and \((-I)^{N-k} f(x) = (-1)^{N-k} f(x)\),

\(\displaystyle \Delta^N f(x) = \sum_{k=0}^N \binom{N}{k} E^k (-I)^{N-k} f(x) = \sum_{k=0}^N \binom{N}{k}(-1)^{N-k}f(x+k)\) holds.

[2-3-1-2] Proof of \(\Delta^N f(x) = B^N N!\)

By the binomial theorem, one can verify that for a positive integer \(m\) and real number \(r\), \(\Delta(rx^m) = rmx^{m-1} + (\text{polynomial of degree} \leq m-2)\). Since the coefficient of \(x^N\) in \(f(x)\) is \(B^N\), we have \(\Delta^N f(x) = B^N N!\).

From [2-3-1-1] and [2-3-1-2], \(\displaystyle \sum_{k=0}^N \binom{N}{k}(-1)^{N-k}f(x+k) = B^N N!\) holds. Substituting \(x = 0\) gives \(\displaystyle \sum_{k=0}^N \binom{N}{k}(-1)^{N-k}f(k) = B^N N!\).

Since a linear combination of \(f(0), f(1), \ldots, f(N)\) equals \(B^N N!\), \(\gcd(f(0), f(1), \ldots, f(N))\) divides \(B^N N!\). Particularly, for a prime \(p\) with \(p \nmid B\), we have \(v_p(G) \leq v_p(B^N N!) = v_p(N!)\).

[2-3-2] Proof of \(v_p(G) \geq v_p(N!)\)

Let \(y_i = Di + (Bn + C)\).

Considering the equation \(y_i \equiv 0 \pmod{p^k}\): since \(D\) is not a multiple of \(p\), there exist integers \(B'\) and \(C'\) such that \(i \equiv \dfrac{-(Bn+C)}{D} \equiv B'n + C' \pmod{p^k}\). Since \(p^k \nmid B'\), among any \(p^k\) consecutive integers, exactly one satisfies this congruence. Thus, among \(y_1, y_2, \ldots, y_N\), at least \(\Bigl\lfloor \dfrac{N}{p^k} \Bigr\rfloor\) are divisible by \(p^k\). Hence, by Legendre’s formula,

\[\displaystyle v_p(f(n)) = v_p(y_1) + v_p(y_2) + \ldots + v_p(y_N) \geq \sum_{k=1}^\infty \Bigl\lfloor \dfrac{N}{p^k} \Bigr\rfloor = v_p(N!).\]

Both \(v_p(G) \leq v_p(N!)\) and \(v_p(G) \geq v_p(N!)\) have been shown, so we conclude \(v_p(G) = v_p(N!)\).

---

[2-4] When \(p \mid B\) and \(p \nmid D\)

Let \(M = v_p(B)\).

For integers \(e\) and \(n\), define \(h(e, n)\) as the number of multiples of \(p^e\) among \(Bn+C, Bn+D+C, \ldots, Bn+(N-1)D+C\). Then \(\displaystyle v_p(f(n)) = \sum_{k=1}^{\infty} h(k, n)\).

We split this as \(\displaystyle \sum_{k=1}^{\infty} h(k,n) = \sum_{k=1}^{M} h(k,n) + \sum_{k=M+1}^{\infty} h(k,n)\).

[2-4-1] \(h(k, n)\) for \(1 \leq k \leq M\)

The value of \(h(k, n)\) is constant regardless of \(n\), and equals the number of multiples of \(p^k\) among \((C, 2D+C, \ldots, (N-1)D+C)\).

[2-4-2] \(h(k, n)\) for \(k > M\)

Let \(B' = \dfrac{B}{p^M}\). The indices \(i\) for which \(Di + C\) is a multiple of \(p^M\) form an arithmetic sequence with common difference \(p^M\), so setting \(i = p^M i' + j\), there exist integers \(C'\) and \(D'\) such that \(Bx + Di + C = p^M B' + D(p^M i' + j) + C = p^M(B'x + Di' + C')\). Since neither \(B'\) nor \(D\) is a multiple of \(p\), the same argument as in [2-3] applies from here. Letting \(N'\) be the number of \(i\) with \(0 \leq i \leq N-1\) for which \(Di + C\) is a multiple of \(p^M\), we have \(\displaystyle \min\left\lbrace \sum_{k=M+1}^{\infty} h(k,n) \mid 0 \leq n \leq N \right\rbrace = v_p(N'!)\).

---

[3] Computing the Answer

Since \(v_p(G) = v_p(N!)\) holds for most primes, we initialize the answer as \(N!\) and compute the difference for primes \(p\) where \(v_p(G) \neq v_p(N!)\). The cases where \(p \mid B\) and \(p \mid D\), or \(p \nmid B\) and \(p \mid D\), are straightforward. When \(p \mid B\) and \(p \nmid D\), we use the computation method described in [2-4]. By solving a linear Diophantine equation, the computation takes \(O(v_p(B) \log N)\) time per prime \(p\).

By precomputing factorials and prime factorizations, we can compute the answer in \(O(\log^2(B+D+N))\) time per test case.

[4] Implementation Example

from math import gcd

mod = 998244353
M = 10**6 + 10
fac = [1] * (M + 1)
inv = [1] * (M + 1)
for i in range(2, M + 1):
    fac[i] = fac[i - 1] * i % mod
    inv[i] = -(mod // i) * inv[mod % i] % mod
sieve = [0] * M
for i in range(2, M):
    if sieve[i] == 0:
        for j in range(i, M, i):
            sieve[j] = i


def factorization(n):
    while n != 1:
        p = sieve[n]
        yield p
        while n % p == 0:
            n //= p


def calc(n, p):
    e = 0
    tmp = n
    while tmp:
        tmp //= p
        e += tmp
    return e


def solve(N, B, C, D):
    g = gcd(B, gcd(C, D))
    if g != 1:
        return pow(g, N, mod) * solve(N, B // g, C // g, D // g) % mod
    ans = fac[N]
    for p in factorization(D):
        ans *= pow(inv[p], calc(N, p), mod)
        ans %= mod
    for p in factorization(B):
        if D % p == 0:
            continue
        ans *= pow(inv[p], calc(N, p), mod)
        ans %= mod
        M = 0
        tmp = B
        while tmp % p == 0:
            tmp //= p
            M += 1
        v_p = 0
        pe = 1
        for e in range(1, M + 1):
            pe *= p
            min_i = (-C * pow(D, -1, pe)) % pe + 1
            res = (N + pe - min_i) // pe
            v_p += res
            if e == M:
                v_p += calc(res, p)
        ans *= pow(p, v_p, mod)
        ans %= mod
    return ans


T = int(input())
for i in range(T):
    N, B, C, D = map(int, input().split())
    print(solve(N, B, C, D))