Solution 1

Using the Euclidean algorithm, the problem can be reduced to the case \(C=0\).

while C != 0:
    int k = floor(A / C)
    A = A - k * C
    B = B - k * D
    swap(A, C)
    swap(B, D)

Let the values of \(A,B,C,D\) after this process be \(A',B',C'(=0),D'\). Then the answer is \(\displaystyle \sum_{i=1}^N \gcd(A'i+B',|D'|)\). Here, \(|D'|\leq 10^8\) (because \(|AD-BC|\) is always constant).

When \(D'=0\), the problem is straightforward.

When \(D'\neq 0\), for each positive divisor \(m\) of \(D'\), we can find the number of \(i\,(1\leq i\leq N)\) such that \(A'i+B'\) is a multiple of \(m\) by solving a linear Diophantine equation. Then, by Moebius inversion over divisors, we can find the number of \(m\) such that \(\gcd(Ai+B,|D'|)=m\).

The number of divisors of an integer up to \(10^8\) is at most \(768\), so even taking quadratic time for the Moebius inversion is fast enough.

Solution 2

Because \(C(Ai+B)-A(Ci+D)=(CB-AD)\), the greatest common divisor is a divisor of \(|CB-AD|\). When \(CB-AD=0\), the problem is straightforward.

For each divisor \(m\) of \(|CB-AD|\), the condition on \(i\) such that both \(Ai+B\) and \(Ci+D\) are multiples of \(m\) can be found by repeatedly solving linear Diophantine equations. The rest follows as in Solution 1.

Sample implementation of Solution 1

from atcoder.math import inv_mod
from math import gcd


def divisors(n):
    lower, upper = [], []
    i = 1
    while i * i <= n:
        if n % i == 0:
            lower.append(i)
            if i * i != n:
                upper.append(n // i)
        i += 1
    return lower + upper[::-1]


mod = 998244353


def solve():
    n, A, B, C, D = map(int, input().split())
    B += A
    D += C
    while C != 0:
        k = A // C
        A -= k * C
        B -= k * D
        A, B, C, D = C, D, A, B
    if D == 0:
        print((n * (n - 1) // 2 * A + B * n) % mod)
        return
    divs = divisors(abs(D))
    sz = len(divs)
    cnt = [0] * sz
    ans = 0
    for i in range(sz - 1, -1, -1):
        m = divs[i]
        g = gcd(A, m)
        if B % g == 0:
            p, q, r = A // g, B // g, m // g
            mn = (-q * inv_mod(p, r)) % r
            cnt[i] = max(0, (n - 1 - mn) // r + 1)
        for j in range(i + 1, sz):
            if divs[j] % m == 0:
                cnt[i] -= cnt[j]
        ans += m * cnt[i]
    print(ans % mod)


for _ in range(int(input())):
    solve()