Problem:
 Score : 100 points    Problem Statement You are given two integer arrays a_0, a_1, ..., a_{N - 1} and b_0, b_1, ..., b_{M - 1} . Calculate the array c_0, c_1, ..., c_{(N - 1) + (M - 1)} , defined by c_i = \sum_{j = 0}^i a_j b_{i - j} \bmod 998244353 .      Constraints  1 \leq N, M \leq 524288  0 \leq a_i, b_i < 998244353  All values in Input are integer.         Input Input is given from Standard Input in the following format:  N  M  a_0  a_1 ... a_{N-1}  b_0  b_1 ... b_{M-1}       Output Print the answer in the following format:  c_0  c_1 ... c_{(N - 1) + (M - 1)}         Sample Input 1 4 5
1 2 3 4
5 6 7 8 9      Sample Output 1 5 16 34 60 70 70 59 36       Sample Input 2 1 1
10000000
10000000      Sample Output 2 871938225

Solution (Python):
class NTT:
    def __init__(self, convolution_rank, binary_level):
        self.cr = convolution_rank
        self.bl = binary_level
        self.modulo_list = [1]*convolution_rank
        self.primal_root_list = [1]*convolution_rank
        self.bin_list = [1]*(binary_level+1)
        self.primal_base_matrix = [[0 for _ in range(binary_level+1)] for __ in range(convolution_rank)]
        self.inverse_base_matrix = [[0 for _ in range(binary_level+1)] for __ in range(convolution_rank)]
        for i in range(binary_level):
            self.bin_list[i+1] = self.bin_list[i] * 2
    
    def doubling(self, n, m, modulo=0):
        y = 1
        tmp = m
        bas = n
        while tmp:
            if tmp % 2:
                y *= bas
                if modulo:
                    y %= modulo
            bas *= bas
            if modulo:
                bas %= modulo
            tmp >>= 1
        return y
    
    def powlimit(self, n):
        y = 1
        cnt = 0
        while y < n:
            y *= 2
            cnt += 1
        return y, cnt
    
    def IsPrime(self, num):
        p = 2
        if num <= 1:
            return False
        while p * p <= num:
            if num % p == 0:
                return False
            p += 1
        return True
    
    def extgcd(self, a, b, c):
        x, y, u, v, k, l = 1, 0, 0, 1, a, b
        while l:
            x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
            k, l = l, k % l
        if c % k:
            return "No Solution"
        return x * c, y * c
    
    def inved(self, a, modulo):
        x, y = self.extgcd(a, modulo, 1)
        return (x+modulo)%modulo
    
    def root_manual(self):
        for i in range(self.cr):
            r = 1
            flg = True
            while flg:
                fflg = True
                for j in range(self.bl):
                    if self.doubling(r, self.bin_list[j], self.modulo_list[i]) == 1:
                        fflg = False
                if self.doubling(r, self.bin_list[-1], self.modulo_list[i]) != 1:
                    fflg = False
                if fflg:
                    flg = False
                else:
                    r += 1
                    if r >= self.modulo_list[i]:
                        break
            self.primal_root_list[i] = r
    
    def make_prime_root(self):
        cnt = 0
        j = 1
        last = self.bin_list[-1]
        while cnt < self.cr:
            if self.IsPrime(j*last+1):
                flg = True
                r = 1
                while flg:
                    fflg = True
                    for i in range(self.bl):
                        if self.doubling(r, self.bin_list[i], j*last+1) == 1:
                            fflg = False
                    if self.doubling(r, last, j*last+1) != 1:
                        fflg = False
                    if fflg:
                        flg = False
                    else:
                        r += 1
                        if r >= j*last+1:
                            break
                if flg==False:
                    self.modulo_list[cnt] = j*last+1
                    self.primal_root_list[cnt] = r
                    cnt += 1
            j += 2
    
    def make_basis(self):
        for i in range(self.cr):
            for j in range(self.bl):
                tmp = self.doubling(2, self.bl-j)
                self.primal_base_matrix[i][j] = self.doubling(self.primal_root_list[i], tmp, self.modulo_list[i])
                self.inverse_base_matrix[i][j] = self.inved(self.primal_base_matrix[i][j], self.modulo_list[i])
    
    def NTT(self, f, n, idx, depth, inverse=False, surface=True):
        mod = self.modulo_list[idx]
        if n == 1:
            return f
        if n == 2:
            xf = [(f[0] + f[1])%mod, (f[0] - f[1])%mod]
            if inverse * surface:
                xf[0] *= self.inved(2, mod)
                xf[0] %= mod
                xf[1] *= self.inved(2, mod)
                xf[1] %= mod
            return xf
        hn = n // 2
        fe = [f[2*i+0] for i in range(hn)]
        fo = [f[2*i+1] for i in range(hn)]
        fe = self.NTT(fe, hn, idx, depth-1, inverse, False)
        fo = self.NTT(fo, hn, idx, depth-1, inverse, False)
        xf = [0 for _ in range(n)]
        grow = 1
        if inverse:
            seed = self.primal_base_matrix[idx][depth]
        else:
            seed = self.inverse_base_matrix[idx][depth]
        for i in range(hn):
            right = (fo[i] * grow) % mod
            xf[i+0*hn] = (fe[i] + right) % mod
            xf[i+1*hn] = (fe[i] - right) % mod
            grow *= seed
            grow %= mod
        if inverse * surface:
            invn = self.inved(n, mod)
            for i in range(n):
                xf[i] *= invn
                xf[i] %= mod
        return xf
    
    def CRT(self, num, a_list, m_list):
        r = a_list[0]
        bas = m_list[0]
        x, y = 0, 0
        for i in range(1, num):
            x, y = self.extgcd(bas, -m_list[i], -a_list[i]-r)
            r += bas * x
            bas *= b[i]
        return r % bas

N, M = map(int, input().split())
A = list(map(int, input().split()))
B = list(map(int, input().split()))
ntt = NTT(1, 23)
MOD = 998244353
ntt.modulo_list[0] = MOD
ntt.root_manual()
ntt.make_basis()
pl, cnt = ntt.powlimit(N+M)
A = A[:] + [0 for _ in range(pl-N)]
B = B[:] + [0 for _ in range(pl-M)]
X = ntt.NTT(A, pl, 0, cnt)
Y = ntt.NTT(B, pl, 0, cnt)
Z = [X[i]*Y[i]%MOD for i in range(pl)]
Z = ntt.NTT(Z, pl, 0, cnt, True)
for i in range(N+M-1):
    print(Z[i], end = ' '*(i!=N+M-2))
print()

Solution (Python):
class FastModuloTransform():
    def __init__(self, n=2**20, mod=998244353, root=3):
        self.n = n
        self.mod = mod
        self.omega = pow(root, (mod - 1) // n, mod)
        self.rev = pow(self.omega, mod - 2, mod)

    def _bit_reverse_(self, A):
        n = self.n
        m = self.n // 2
        i = 0
        for j in range(0, m, 2):
            if j < i:
                A[i], A[j] = A[j], A[i]
                A[i + m + 1], A[j + m + 1] = A[j + m + 1], A[i + m + 1]
            A[i + 1], A[j + m] = A[j + m], A[i + 1]
            k = m // 2
            i ^= k
            while k > i:
                k //= 2
                i ^= k
        return A

    def _fmt_(self, A, base):
        halfpow = pow(base, self.n // 2, self.mod)
        n = self.n
        m = 1
        while n > 1:
            n //= 2
            w = pow(base, n, self.mod)
            k = 1
            for j in range(m):
                for i in range(j, self.n, 2 * m):
                    u = A[i]
                    v = (A[i + m] * k) % self.mod
                    A[i] = (u + v) % self.mod
                    A[i + m] = (u + v * halfpow) % self.mod
                k = (k * w) % self.mod
            m *= 2
        return A

    def fmt(self, f):
        A = f + [0] * (self.n - len(f))
        self._bit_reverse_(A)
        F = self._fmt_(A, self.omega)
        return F

    def ifmt(self, F):
        A = F + [0] * (self.n - len(F))
        self._bit_reverse_(A)
        nrev = pow(self.n, self.mod - 2, self.mod)
        f = [(e * nrev) % self.mod for e in self._fmt_(A, self.rev)]
        return f

    def convolute(self, f, g):
        F = self.fmt(f)
        G = self.fmt(g)
        H = [(s * t) % self.mod for s, t in zip(F, G)]
        h = self.ifmt(H)
        return h

    def autocorrelation(self, f):
        F = self.fmt(f)
        H = [(s * s) % self.mod for s in F]
        h = self.ifmt(H)
        return h

fmt = FastModuloTransform()

N, M = map(int, input().split())
A = list(map(int, input().split()))
B = list(map(int, input().split()))

res = fmt.convolute(A, B)[:N + M - 1]

print(*res)

Solution (Python):

def gcd(x, y):
    k, l = x, y
    if y < 0:
        k, l = -k, -l
    while l:
        k, l = l, k % l
    return k
    
def extgcd(a, b, c):
    if b == 0:
        return c // a, 0
    if b < 0:
        a, b, c = -a, -b, -c
    x, y, u, v, k, l = 1, 0, 0, 1, a, b
    while l != 0:
        x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
        k, l = l, k % l
    return c * x, c * y


def IsPrime(n):
    if n == 1:
        return False
    p = 2
    while p * p <= n:
        if n % p == 0:
            return False
        p += 1
    return True


def fracsum(f1, f2, sign = 1):
    fs = [f1[0]*f2[1]+sign*f1[1]*f2[0], f1[1]*f2[1]]
    g = gcd(fs[0], fs[1])
    fs[0] //= g
    fs[1] //= g
    return fs
def fracprod(f1, f2):
    fs = [f1[0]*f2[0], f1[1]*f2[1]]
    g = gcd(fs[0], fs[1])
    fs[0] //= g
    fs[1] //= g
    return fs


def factorization(n):
    D = {}
    a = n
    p = 2
    while a != 1:
        cnt = 0
        while a % p == 0:
            cnt += 1
            a //= p
        if cnt != 0:
            D[p] = cnt
        p += 1
        if p * p > n and a != 1:
            D[a] = 1
            break
    return D
def divisors_list(D):
    F = [1]
    for i in D:
        mL = F[:]
        dp = i
        for j in range(D[i]):
            for k in mL:
                F.append(dp*k)
            dp *= i
    F.sort()
    return F

def doubling(n, m, modulo=0):
    y = 1
    base = n
    tmp = m
    while tmp:
        if tmp % 2 == 1:
            y *= base
            if modulo > 0:
                y %= modulo
        base *= base
        if modulo > 0:
            base %= modulo
        tmp //= 2
    return y

def inved(a, modulo):
    x, y, u, v, k, l = 1, 0, 0, 1, a, modulo
    while l != 0:
        x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
        k, l = l, k % l
    return x % modulo
def bic(n): 
    cnt = 0
    while n % 2 == 0:
        cnt += 1
        n //= 2
    return cnt
def powlimit(n):
    y = 1
    cnt = 0
    while y < n:
        y *= 2
        cnt += 1
    return y, cnt
def CRT(n, a, b): 
    r = a[0]
    bas = b[0]
    x, y = 0, 0
    for i in range(1, n):
        x, y = extgcd(bas, -b[i], a[i] - r)
        r += bas * x
        bas *= b[i]
    return r % bas


def mod_maker(cnt, num):
    bin_list = [1 for _ in range(cnt+1)]
    for i in range(cnt, 0, -1):
        bin_list[i-1] = 2*bin_list[i]
    L = []
    pm = []
    bas = doubling(2, cnt)
    ccnt = 0
    j = 1
    while ccnt < num:
        if IsPrime(j*bas+1):
            x = j*bas+1
            L.append(j*bas+1)
            ccnt += 1
            st = 1
            flg = 0
            while not flg:
                flg = (doubling(st, bin_list[0], x) == 1)
                for k in range(cnt):
                    flg *= (doubling(st, bin_list[k+1], x) != 1)
                if flg:
                    pm.append(st)
                    break
                st += 1
        j += 2
    return L, pm

def NTT(f, n, cn, num, inverse=False):
    if n == 1:
        return f
    if n == 2:
        return [(f[0]+f[1])%mod[num], (f[0]-f[1])%mod[num]]
    dn = n // 2
    fe = [f[2*i+0] for i in range(dn)]
    fo = [f[2*i+1] for i in range(dn)]
    fe = NTT(fe, dn, cn-1, num, inverse)
    fo = NTT(fo, dn, cn-1, num, inverse)
    seed = inverse * prim[num][cn] + (1 - inverse) * invp[num][cn]
    grow = 1
    xf = [0 for _ in range(n)]
    for i in range(dn):
        right = fo[i] * grow % mod[num]
        xf[i+0*dn] = (fe[i] + right) % mod[num]
        xf[i+1*dn] = (fe[i] - right) % mod[num]
        grow *= seed
        grow %= mod[num]
    return xf


N, M = map(int, input().split())
A = list(map(int, input().split()))
B = list(map(int, input().split()))

pl, bn = powlimit(2*max(N, M))
A = A[:] + [0]*(pl-N)
B = B[:] + [0]*(pl-M)
convolution_time = 1 
binary_count = 23
mod, primod = [998244353], [31]
prim = [[doubling(primod[j], doubling(2, binary_count-i), mod[j]) for i in range(binary_count)] for j in range(convolution_time)]
invp = [[inved(prim[j][i], mod[j]) for i in range(binary_count)] for j in range(convolution_time)]

X1 = NTT(A, pl, bn, 0)
X2 = NTT(B, pl, bn, 0)
XX = [X1[i]*X2[i]%mod[0] for i in range(pl)]
Y = NTT(XX, pl, bn, 0, 1)
ipl = inved(pl, mod[0])
for i in range(N+M-1):
    Y[i] *= ipl
    Y[i] %= mod[0]
    print(Y[i], end = ' '*(i!=N+M-2))
print()

Solution (Python):
import sys
import numpy as np
import numba
from numba import njit, b1, i4, i8, f8

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 998_244_353
MOD_root = 3
max_ntt = 20


@njit((i8, i8), cache=True)
def mpow(a, n):
    p = 1
    while n:
        if n & 1:
            p = p * a % MOD
        a = a * a % MOD
        n >>= 1
    return p


@njit((i8, i8), cache=True)
def ntt_precompute(N, primitive_root):
    roots = np.empty(1 << N, np.int64)
    iroots = np.empty(1 << N, np.int64)
    roots[0] = iroots[0] = 1
    for n in range(N):
        x = mpow(primitive_root, (MOD - 1) >> n + 2)
        y = mpow(x, MOD - 2)
        for i in range(1 << n):
            roots[(1 << n) | i] = roots[i] * x % MOD
            iroots[(1 << n) | i] = iroots[i] * y % MOD
    return roots, iroots


@njit((i8[:], i8[:], i8[:], numba.b1), cache=True)
def ntt(roots, iroots, A, inverse):
    N = len(A)
    if not inverse:
        m = N >> 1
        while m:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m] * roots[k]
                    A[i], A[i + m] = x + y, x - y
            A %= MOD
            m >>= 1
    else:
        m = 1
        while m < N:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m]
                    A[i], A[i + m] = x + y, (x - y) * iroots[k]
            A %= MOD
            m <<= 1
        for i in range(N):
            A[i] = A[i] * mpow(N, MOD - 2) % MOD
    return A


def ntt_convolve(roots, iroots, A, B):
    N = len(A) + len(B) - 1
    L = 1
    while 2 * L < N:
        L *= 2
    L *= 2
    A = np.concatenate((A, np.zeros(L - len(A), np.int64)))
    B = np.concatenate((B, np.zeros(L - len(B), np.int64)))
    ntt(roots, iroots, A, False)
    ntt(roots, iroots, B, False)
    A *= B
    A %= MOD
    ntt(roots, iroots, A, True)
    return A[:N]


roots, iroots = ntt_precompute(max_ntt, MOD_root)

readline()
A = np.array(readline().split(), np.int64)
B = np.array(readline().split(), np.int64)

C = ntt_convolve(roots, iroots, A, B)
print(' '.join(map(str, C.tolist())))

Solution (Python):





import numpy as np

fft_len = 2 ** 20   
MOD = 998244353     

def convolve10(As, Bs):
    "As, Bsの要素が10bitに入る(つまり2**10 = 1024以下)なら誤差はまず平気そう。"
    return np.round(np.real(np.fft.irfft(np.fft.rfft(As, fft_len) * np.fft.rfft(Bs, fft_len), fft_len))).astype(np.int64)

def convolve30(As, Bs):
    """
    10bitで3分割してToom-3
    参考: https://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication
    """
    mask = (1 << 10) - 1
    
    m0 = As & mask
    m1 = (As >> 10) & mask
    m2 = As >> 20

    p_tmp = m0 + m2 
    p0 = m0
    p1 = p_tmp + m1
    pm1 = p_tmp - m1
    pm2 = ((pm1 + m2) << 1) - m0
    pinf = m2

    n0 = Bs & mask
    n1 = (Bs >> 10) & mask
    n2 = Bs >> 20

    q_tmp = n0 + n2 
    q0 = n0
    q1 = q_tmp + n1
    qm1 = q_tmp - n1
    qm2 = ((qm1 + n2) << 1) - n0
    qinf = n2

    r0   = convolve10(  p0,   q0)
    r1   = convolve10(  p1,   q1)
    rm1  = convolve10( pm1,  qm1)
    rm2  = convolve10( pm2,  qm2)
    rinf = convolve10(pinf, qinf)

    r_0 = r0
    r_4 = rinf
    r_3 = (rm2 - r1) // 3
    r_1 = (r1 - rm1) >> 1
    r_2 = (rm1 - r0)
    r_3 = ((r_2 - r_3) >> 1) + (rinf << 1)
    r_2 += r_1 - r_4
    r_1 -= r_3
    
    ret = ((r_4 << 20) + (r_3 << 10) + r_2) % MOD
    ret = ((ret << 20) + (r_1 << 10) + r_0) % MOD

    return ret


if __name__ == "__main__":  
    import sys
    input = sys.stdin.buffer.readline

    N, M = map(int, input().split())
    fft_len = 1 << (N + M - 1).bit_length()
    As = np.array(list(map(int, input().split())), dtype = np.int64)
    Bs = np.array(list(map(int, input().split())), dtype = np.int64)
    Cs = convolve30(As, Bs)
    print(*list(Cs[:N + M - 1]))

Solution (Python):
import sys
readline = sys.stdin.readline

MOD = 998244353
def frac(limit):
    frac = [1]*limit
    for i in range(2,limit):
        frac[i] = i * frac[i-1]%MOD
    fraci = [None]*limit
    fraci[-1] = pow(frac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        fraci[i] = fraci[i+1] * (limit + i + 1) % MOD
    return frac, fraci
frac, fraci = frac(1341398)
class NTT:
    def __init__(self, MOD = 998244353, pr = 3, LS = 20):
        self.MOD = MOD
        self.N0 = 1<<LS
        omega = pow(pr, (MOD-1)//self.N0, MOD)
        omegainv = pow(omega, MOD-2, MOD)
        self.w = [0]*(self.N0//2)
        self.winv = [0]*(self.N0//2)
        self.w[0] = 1
        self.winv[0] = 1
        for i in range(1, self.N0//2):
            self.w[i] = (self.w[i-1]*omega)%MOD
            self.winv[i] = (self.winv[i-1]*omegainv)%MOD
        used = set()
        for i in range(self.N0//2):
            if i in used:
                continue
            j = 0
            for k in range(LS-1):
                j |= (i>>k&1) << (LS-2-k)
            used.add(j)
            self.w[i], self.w[j] = self.w[j], self.w[i]
            self.winv[i], self.winv[j] = self.winv[j], self.winv[i]
        
    def _fft(self, A):
        MOD = self.MOD
        M = len(A)
        bn = 1
        hbs = M>>1
        while hbs:
            for j in range(hbs):
                A[j], A[j+hbs] = A[j] + A[j+hbs], A[j] - A[j+hbs]
                if A[j] > MOD:
                    A[j] -= MOD
                if A[j+hbs] < 0:
                    A[j+hbs] += MOD
            for bi in range(1, bn):
                wi = self.w[bi]
                for j in range(bi*hbs*2, bi*hbs*2+hbs):
                    A[j], A[j+hbs] = (A[j] + wi*A[j+hbs])%MOD, (A[j] - wi*A[j+hbs])%MOD
            bn <<= 1
            hbs >>= 1

    def _ifft(self, A):
        MOD = self.MOD
        M = len(A)
        bn = M>>1
        hbs = 1
        while bn:
            for j in range(hbs):
                A[j], A[j+hbs] = A[j] + A[j+hbs], A[j] - A[j+hbs]
                if A[j] > MOD:
                    A[j] -= MOD
                if A[j+hbs] < 0:
                    A[j+hbs] += MOD
            for bi in range(1, bn):
                winvi = self.winv[bi]
                for j in range(bi*hbs*2, bi*hbs*2+hbs):
                    A[j], A[j+hbs] = (A[j] + A[j+hbs])%MOD, winvi*(A[j] - A[j+hbs])%MOD
            bn >>= 1
            hbs <<= 1
    
    def convolve(self, A, B):
        LA = len(A)
        LB = len(B)
        LC = LA+LB-1
        M = 1<<(LC-1).bit_length()
        A += [0]*(M-LA)
        B += [0]*(M-LB)
        self._fft(A)
        self._fft(B)
        C = [0]*M
        for i in range(M):
            C[i] = A[i]*B[i]%self.MOD
        self._ifft(C)
        minv = pow(M, self.MOD-2, self.MOD)
        for i in range(LC):
            C[i] = C[i]*minv%self.MOD
        return C[:LC]
        
        return C
    
    def inverse(self, A):
        LA = len(A)
        dep = (LA-1).bit_length()
        M = 1<<dep
        A += [0]*(M-LA)
        
        g = [pow(A[0], self.MOD-2, self.MOD)]
        for n in range(dep):
            dl = 1<<(n+1)
            f = A[:dl]
            fg = self.convolve(f, g[:])[:dl]
            fgg = self.convolve(fg, g[:])[:dl]
            ng = [None]*dl
            for i in range(dl//2):
                ng[i] = (2*g[i]-fgg[i])%self.MOD
            for i in range(dl//2, dl):
                ng[i] = self.MOD-fgg[i]
            g = ng[:]
        return g[:LA]
    
    def integral(self, A):
        A = [1] + [A[i]*fraci[i+2]%self.MOD for i in range(len(A))]

N, M = map(int, readline().split())
A = list(map(int, readline().split()))
B = list(map(int, readline().split()))
T = NTT()
print(*T.convolve(A, B))

Solution (Python):
import sys
input = lambda : sys.stdin.readline().rstrip()
sys.setrecursionlimit(max(1000, 10**9))
write = lambda x: sys.stdout.write(x+"\n")


n,m = list(map(int, input().split()))
a = list(map(int, input().split()))
b = list(map(int, input().split()))

import numpy as np
M = 998244353
a = np.array(a, dtype=np.int64)
b = np.array(b, dtype=np.int64)
def fft(a,b):
    l = 1
    while 2 * l < len(a) + len(b) - 1:
        l *= 2
    l *= 2
    c = np.fft.irfft((np.fft.rfft(a,l))*(np.fft.rfft(b,l)),l)
    c = np.rint(c).astype(np.int64)
    return c
def fft_large(a,b):
    d = 30000
    a1, a2 = np.divmod(a,d)
    b1, b2 = np.divmod(b,d)
    aa = fft(a1,b1) % M
    bb = fft(a2,b2) % M
    cc = (fft(a1+a2, b1+b2) - (aa+bb)) % M
    h = (((aa*d)%M)*d  + cc*d + bb) % M
    return h
def fft_large(a,b):
    d = 1<<10
    a1, a2 = np.divmod(a,d*d)
    a2, a3 = np.divmod(a2,d)
    b1, b2 = np.divmod(b,d*d)
    b2, b3 = np.divmod(b2,d)
    aa = fft(a1,b1) % M
    bb = fft(a2,b2) % M
    cc = fft(a3,b3) % M
    dd = (fft(a1+a2, b1+b2) - (aa+bb)) % M
    ee = (fft(a2+a3, b2+b3) - (bb+cc)) % M
    ff = (fft(a1+a3, b1+b3) - (aa+cc)) % M
    h = (((aa*d*d)%M)*d*d + ((dd*d*d)%M)*d + (bb+ff)*d*d + ee*d + cc) % M
    return h

c = fft_large(a,b)[:n+m-1]
write(" ".join(map(str, c)))

Solution (Python):
mod = 998244353
omega = pow(3,119,mod)
rev_omega = pow(omega,mod-2,mod)

N = 2*10**5
g1 = [1]*(N+1) 
g2 = [1]*(N+1) 
inv = [1]*(N+1) 

for i in range( 2, N + 1 ):
    g1[i]=( ( g1[i-1] * i ) % mod )
    inv[i]=( ( -inv[mod % i] * (mod//i) ) % mod )
    g2[i]=( (g2[i-1] * inv[i]) % mod )
inv[0]=0

def _ntt(f,L,reverse=False):
    F=[f[i] for i in range(L)]
    n = L.bit_length() - 1
    base = omega
    if reverse:
        base = rev_omega

    if not n:
        return F

    size = 2**n
    wj = pow(base,2**22,mod)
    res = [0]*2**n

    for i in range(n,0,-1):
        use_omega = pow(base,2**(22+i-n),mod)
        res = [0]*2**n
        size //= 2
        w = 1
        for j in range(0,L//2,size):
            for a in range(size):
                res[a+j] = (F[a+2*j] + w * F[a+size+2*j]) % mod
                t = (w * wj) % mod
                res[L//2+a+j] = (F[a+2*j] + t * F[a+size+2*j]) % mod
            w = (w * use_omega) % mod
        F = res

    return res

def ntt(f,L=0):
    l = len(f)
    if not L:
        L = 1<<((l-1).bit_length())
    while len(f)<L:
        f.append(0)
    f=f[:L]
    F = _ntt(f,L)
    return F

def intt(f,L=0):
    l = len(f)
    if not L:
        L = 1<<((l-1).bit_length())
    while len(f)<L:
        f.append(0)
    f=f[:L]
    F = _ntt(f,L,reverse=True)
    inv = pow(L,mod-2,mod)
    for i in range(L):
        F[i] *= inv
        F[i] %= mod
    return F

def convolve(f,g,limit):
    l = len(f)+len(g)-1
    L = 1<<((l-1).bit_length())

    F = ntt(f,L)
    G = ntt(g,L)

    H = [(F[i] * G[i]) % mod for i in range(L)]

    h = intt(H,L)

    return h[:limit]

N,M = map(int,input().split())
a = list(map(int,input().split()))
b = list(map(int,input().split()))

print(*convolve(a,b,N+M-1))

Solution (Python):
def convolve(a, b):
    def fft(f):
        d = n // 2
        v = w
        while d >= 1:
            u = 1
            for i in range(d):
                for j in range(i, n, 2*d):
                    f[j], f[j+d] = (f[j] + f[j+d]) % p, u * (f[j] - f[j+d]) % p
                u = u * v % p
            v = v * v % p
            d //= 2
 
    def ifft(f):
        d = 1
        while d < n:
            v = pow(invw, n // (2 * d), p)
            u = 1
            for i in range(d):
                for j in range(i, n, 2*d):
                    f[j+d] *= u
                    f[j], f[j+d] = (f[j] + f[j+d]) % p, (f[j] - f[j+d]) % p
                u = u * v % p
            d *= 2
 
    p, g = 998244353, 3
    n0, n1 = len(a), len(b)
    n = 1 << (max(n0, n1) - 1).bit_length() + 1
    a = a + [0] * (n-n0)
    b = b + [0] * (n-n1)
    w = pow(g, (p - 1) // n, p)
    invw = pow(w, p-2, p)
    fft(a), fft(b)
    for i in range(n):
        a[i] = a[i] * b[i] % p
    ifft(a)
    invn = pow(n, p - 2, p)
    return [a[i] * invn % p for i in range(n0 + n1 - 1)]

N, M = map(int, input().split())
A = [int(a) for a in input().split()]
B = [int(a) for a in input().split()]
print(*convolve(A, B))

Solution (Python):
import sys
input = sys.stdin.buffer.readline


def convolve(a, b):
    """配列a, bの畳み込みをNTT(number-theoretical transform)によって計算する。
    res[k] = (sum{i = 0..k} a[i] * b[k - i]) % MOD
    """
    MOD = 998244353
    ROOT = 5

    def ntt(a, inverse=False):
        
        for i in range(n):
            j = 0
            for k in range(log_sz):
                j |= (i >> k & 1) << (log_sz - 1 - k)
            if i < j:
                a[i], a[j] = a[j], a[i]

        h = pow(ROOT, (MOD - 1) * pow(n, MOD - 2, MOD) % MOD, MOD)
        if inverse:
            h = pow(h, MOD - 2, MOD)

        
        m = 1
        while m < n:
            zeta_pow = 1
            zeta = pow(h, n // (2 * m), MOD)
            for j in range(m):
                for k in range(0, n, 2 * m):
                    s = a[j + k]
                    t = a[j + k + m] * zeta_pow
                    a[j + k] = (s + t) % MOD
                    a[j + k + m] = (s - t) % MOD
                zeta_pow *= zeta
                zeta_pow %= MOD
            m <<= 1

        
        if inverse:
            n_inv = pow(n, MOD - 2, MOD)
            for i in range(n):
                a[i] *= n_inv
                a[i] %= MOD
        return a

    def intt(a):
        return ntt(a, inverse=True)

    n = 1 << (len(a) + len(b) - 1).bit_length()
    log_sz = n.bit_length() - 1
    a = a + [0] * (n - len(a))
    b = b + [0] * (n - len(b))

    
    ntt(a), ntt(b)
    for i, val in enumerate(b):
        a[i] *= val
        a[i] %= MOD
    intt(a)

    return a


n, m = map(int, input().split())
a = list(map(int, input().split()))
b = list(map(int, input().split()))


print(*convolve(a, b)[:n + m - 1])