[0] 記号の定義

• \(f(x)=\displaystyle \prod_{i=0}^{N-1} (Bx+Di+C)\) とします．
• \(G=\gcd(f(0), f(1), \ldots,f(N))\) とします．
• 正整数 \(n\) および素数 \(p\) に対して，\(n\) を \(p\) で割り切れる回数を \(v_p(n)\) とします．
• 整数 \(a,b(a\neq 0)\) に対し，\(b\) が \(a\) の倍数のとき \(a\mid b\)，そうでないとき \(a \nmid b\) とします．

---

[1] \(\gcd(B, C, D)=1\) への帰着

\(g=\gcd(B, C, D)\) とします．\(g\neq 1\) の場合は \(f(x)=g^N \displaystyle \prod_{i=0}^{N-1} \left(\frac{B}{g}x+\frac{D}{g}i+\frac{C}{g}\right)\) となるため，\(B\to \dfrac{B}{g}, C\to \dfrac{C}{g}, D\to \dfrac{D}{g}\) とした問題の答えを \(g^N\) 倍した値が答えになります．

以降の議論では，\(g=1\) を仮定します．

---

[2] 各素数 \(p\) に対する \(v_p(G)\)

各素数 \(p\) に対して \(v_p(G)\) を考えます．結論を先に述べます．

1. \(p\mid B\) かつ \(p\mid D\) のとき：\(v_p(G)=0\)
2. \(p\nmid B\) かつ \(p\mid D\) のとき：\(v_p(G)=0\)
3. \(p\nmid B\) かつ \(p\nmid D\) のとき：\(v_p(G)=v_p(N!)\)
4. \(p\mid B\) かつ \(p\nmid D\) のとき：計算によって \(v_p(G)\) を高速に求められる

---

[2-1] \(p\mid B\) かつ \(p\mid D\) のとき

仮定より \(p\nmid C\) です．したがって，任意の \(x,i\) に対して \(p\nmid (Bx+Di+C)\) なので \(v_p(G)=0\) です．

---

[2-2] \(p\nmid B\) かつ \(p\mid D\) のとき

\(n=0,1\) の少なくとも一方に対して \(p\nmid (Bn+C)\) です．そのような \(n\) に対して \(p\nmid f(n)\) なので \(v_p(G)=0\) です．

---

[2-3] \(p\nmid B\) かつ \(p\nmid D\) のとき

\(v_p(G)=v_p(N!)\) となることを証明します．

[2-3-1] \(v_p(G)\leq v_p(N!)\) の証明

関数 \(f(x)\) に対して，以下の作用素を定義します．

• 恒等作用素 \(I\)：\(If(x)=f(x)\)
• シフト作用素 \(E\)：\(Ef(x)=f(x+1)\)
• 差分作用素 \(\Delta\)：\(\Delta f(x)=f(x+1)-f(x)(=(E-I)f(x))\)

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

二項定理より \(\displaystyle \Delta^N=(E-I)^N=\sum_{k=0}^N \binom{N}{k}E^k(-I)^{N-k}\) です．

\(E^k f(x)=f(x+k), (-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)\) が成り立ちます．

[2-3-1-2] \(\Delta^N f(x)=B^N N!\) の証明

正整数 \(m\) および実数 \(r\) に対して，\(\Delta (rx^m)=rmx^{m-1} + (m-2 次以下の多項式)\) になることが二項定理により確かめられます．\(f(x)\) の \(x^N\) の係数は \(B^N\) なので，\(\Delta^N f(x)=B^N N!\) です．

[2-3-1-1] と [2-3-1-2] より \(\displaystyle\sum_{k=0}^N \binom{N}{k}(-1)^{N-k}f(x+k)=B^N N!\) が成り立ちます．\(x=0\) を代入すると \(\displaystyle\sum_{k=0}^N \binom{N}{k}(-1)^{N-k}f(k)=B^NN!\) です．

\(f(0),f(1),\ldots,f(N)\) の線形結合が \(B^N N!\) であることより，\(\gcd(f(0),f(1),\ldots, f(N))\) は \(B^N N!\) の約数です．特に，\(p\nmid B\) を満たす素数 \(p\) について，\(v_p(G)\leq v_p(B^N N!)=v_p(N!)\) です．

[2-3-2] \(v_p(G)\geq v_p(N!)\) の証明

\(y_i=Di+(Bn+C)\) とします．

ここで方程式 \(y_i\equiv 0 \pmod {p^k}\) を考えると，\(D\) は \(p\) の倍数でないので，ある整数 \(B',C'\) が存在して \(i\equiv\dfrac{-(Bn+C)}{D} \equiv B^{\prime}n+C^{\prime}\pmod {p^k}\) です．\(p^k\nmid B'\) なので，この合同式を満たす \(i\) は，連続する \(p^k\) 個の整数のなかにちょうど \(1\) つ含まれます．よって，\(y_1,y_2,\ldots,y_N\) の中に \(p^k\) で割り切れるものは少なくとも \(\Bigl\lfloor\dfrac{N}{p^k}\Bigr\rfloor\) 個存在します．よって，ルジャンドルの定理より

\[\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!)\]

です．

以上より，\(v_p(G)\leq v_p(N!)\) と \(v_p(G)\geq v_p(N!)\) が示されたので，\(v_p(G)=v_p(N!)\) です．

---

[2-4] \(p\mid B\) かつ \(p\nmid D\) のとき

\(M=v_p(B)\) とします．

整数 \(e, n\) に対し，関数 \(h(e, n)\) を，\(Bn+C, Bn+D+C, \dots, Bn+(N-1)D+C\) の中にある \(p^e\) の倍数の個数とします．\(\displaystyle v_p(f(n))=\sum_{k=1}^{\infty} h(k,n)\) です．

ここで，\(\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] \(1\leq k\leq M\) のときの \(h(k,n)\)

\(h(k,n)\) の値は \(n\) によらず常に一定で，その値は \((C, 2D+C, \ldots,(N-1)D+C)\) の中にある \(p^k\) の倍数の個数です．

[2-4-2] \(k\gt M\) のときの \(h(k,n)\)

\(B'=\dfrac{B}{p^M}\) とします．\(Di+C\) が \(p^M\) の倍数となるような \(i\) は公差 \(p^M\) の等差数列になるので，\(i=p^Mi'+j\) と置くと，ある整数 \(C', D'\) を用いて \(Bx+Di+C=p^MB'+D(p^M i'+j)+C=p^M(B'x+Di'+C')\) となります．\(B',D\) はともに \(p\) の倍数ではないため，ここから先は [2-3] と同じ議論ができます．\(Di+C\) が \(p^M\) の倍数となるような \(i (0\leq i\leq N-1)\) の個数を \(N'\) とすると，\(\displaystyle \min\left\lbrace\sum_{k=M+1}^{\infty} h(k,n)\mid 0\leq n\leq N\right\rbrace=v_p(N'!)\) です．

---

[3] 答えの計算

ほとんどの素数に対して \(v_p(G)=v_p(N!)\) が成り立つため，答えの初期値を \(N!\) としておき，\(v_p(G)\neq v_p(N!)\) であるような \(p\) に対して，差分を計算することにします．\(p\mid B\) かつ \(p\mid D\) のときや，\(p\nmid B\) かつ \(p\mid D\) のときは簡単です．\(p\mid B\) かつ \(p\nmid D\) のときは，[2-4] で書いた計算方法で求められます．一次不定方程式を解くことにより，素数 \(p\) に対して \(O(v_p(B)\log N)\) 時間で計算できます．

階乗や素因数分解を高速に求められるように前計算しておくことで，テストケースあたり \(O(\log^2(B+D+N))\) 時間で答えを求められます．

[4] 実装例

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):
    # return v_p(n!)
    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))