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



Sample input:
4 5
1 2 3 4
5 6 7 8 9


Sample output:
5 16 34 60 70 70 59 36

Submission status: Time Limit Exceeded

Solution code:
# ----------- gcd ----------- #
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

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

# --- fraction calculator --- #
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

# ------ factorization ------ #
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): # count time divided by 2
    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): # x = a[i] mod b[i]
    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

# making lists of mod and prim. root
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

# ----------- NTT test zone ----------- #
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)+2)
A = A[:] + [0]*(pl-N)
B = B[:] + [0]*(pl-M)
convolution_time = 1 # convolution level (number of 'mod's)
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, j) for j in range(convolution_time)]
X2 = [NTT(B, pl, bn, j) for j in range(convolution_time)]
XX = [[X1[j][i]*X2[j][i]%mod[j] for i in range(pl)] for j in range(convolution_time)]
Y = [NTT(XX[j], pl, bn, j, 1) for j in range(convolution_time)]
Z = [0 for _ in range(pl)]
ipl = [inved(pl, mod[j]) for j in range(convolution_time)]
for j in range(convolution_time):
    for i in range(pl):
        Y[j][i] *= ipl[j]
        Y[j][i] %= mod[j]
for i in range(N+M-1):
    print(Y[0][i], end = ' '*(i!=N+M-2))
print()

Submission status: Time Limit Exceeded

Solution code:
# ----------- gcd ----------- #
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

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

# --- fraction calculator --- #
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

# ------ factorization ------ #
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): # count time divided by 2
    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): # x = a[i] mod b[i]
    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

# making lists of mod and prim. root
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

# ----------- NTT test zone ----------- #
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)+2)
A = A[:] + [0]*(pl-N)
B = B[:] + [0]*(pl-M)
convolution_time = 1 # convolution level (number of 'mod's)
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, j) for j in range(convolution_time)]
X2 = [NTT(B, pl, bn, j) for j in range(convolution_time)]
XX = [[X1[j][i]*X2[j][i]%mod[j] for i in range(pl)] for j in range(convolution_time)]
Y = [NTT(XX[j], pl, bn, j, 1) for j in range(convolution_time)]
Z = [0 for _ in range(pl)]
ipl = [inved(pl, mod[j]) for j in range(convolution_time)]
for j in range(convolution_time):
    for i in range(pl):
        Y[j][i] *= ipl[j]
        Y[j][i] %= mod[j]
for i in range(pl):
    Z[i] = CRT(convolution_time, [Y[j][i] for j in range(convolution_time)], mod)
    if i <= N+M-2:
        print(Z[i], end = ' '*(i!=N+M-2))
print()

Submission status: Time Limit Exceeded

Solution code:
# ----------- gcd ----------- #
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

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

# --- fraction calculator --- #
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

# ------ factorization ------ #
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): # count time divided by 2
    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): # x = a[i] mod b[i]
    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

# making lists of mod and prim. root
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

# ----------- NTT test zone ----------- #
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)+2)
A = A[:] + [0]*(pl-N)
B = B[:] + [0]*(pl-M)
convolution_time = 1 # convolution level (number of 'mod's)
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, j) for j in range(convolution_time)]
X2 = [NTT(B, pl, bn, j) for j in range(convolution_time)]
XX = [[X1[j][i]*X2[j][i]%mod[j] for i in range(pl)] for j in range(convolution_time)]
Y = [NTT(XX[j], pl, bn, j, 1) for j in range(convolution_time)]
Z = [0 for _ in range(pl)]
ipl = [inved(pl, mod[j]) for j in range(convolution_time)]
for j in range(convolution_time):
    for i in range(pl):
        Y[j][i] *= ipl[j]
        Y[j][i] %= mod[j]
for i in range(pl):
    Z[i] = CRT(convolution_time, [Y[j][i] for j in range(convolution_time)], mod)
    if i <= N+M-2:
        print(Z[i], end = ' '*(i!=N+M-2))
print()

Submission status: Time Limit Exceeded

Solution code:
k, K = 80, 1 << 80
nu = lambda L: int("".join([bin(a)[2:].zfill(k) for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:].zfill((len(A)+len(B)-1)*k) + "_"
li = lambda s: [int(a, 2) % 998244353 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(len(A)+len(B)-1)]]
input()
A = [int(a) for a in input().split()]
B = [int(a) for a in input().split()]
print(*li(st(nu(A) * nu(B))))

Submission status: Accepted

Solution code:
import numpy as np

MOD = 998244353

N, M = map(int, input().split())
As = list(map(int, input().split()))
Bs = list(map(int, input().split()))

def divideArr3(As, numD):
    arrA = np.array(As, dtype=np.int64)
    mask = (1<<numD) - 1
    arrA0 = arrA & mask
    arrA1 = (arrA>>numD) & mask
    arrA2 = arrA>>(2*numD)
    return arrA2, arrA1, arrA0

def convolveArr(arrA, arrB):
    L = len(arrA) + len(arrB) - 1
    lenFFT = 1 << ((L-1).bit_length())
    fA = np.fft.rfft(arrA, n=lenFFT)
    fB = np.fft.rfft(arrB, n=lenFFT)
    arrC = np.rint(np.fft.irfft(fA * fB, n=lenFFT)).astype(np.int64)
    return arrC[:L]

def convolveMOD(As, Bs, MOD):
    numD = 10
    div = 1<<numD
    arrA2, arrA1, arrA0 = divideArr3(As, numD)
    arrB2, arrB1, arrB0 = divideArr3(Bs, numD)

    arrC4 = convolveArr(arrA2, arrB2) % MOD
    arrC2 = convolveArr(arrA1, arrB1) % MOD
    arrC0 = convolveArr(arrA0, arrB0) % MOD
    arrC3 = (convolveArr(arrA2+arrA1, arrB2+arrB1) % MOD - arrC4 - arrC2) % MOD
    arrC1 = (convolveArr(arrA1+arrA0, arrB1+arrB0) % MOD - arrC2 - arrC0) % MOD
    arrC2 += (convolveArr(arrA2+arrA0, arrB2+arrB0) % MOD - arrC4 - arrC0) % MOD
    arrC2 %= MOD

    div2 = div*div%MOD
    div3 = div2*div%MOD
    div4 = div3*div%MOD
    arrAns = ((arrC4*div4)%MOD + (arrC3*div3)%MOD + (arrC2*div2)%MOD + (arrC1*div)%MOD + arrC0) % MOD
    return arrAns.tolist()

anss = convolveMOD(As, Bs, MOD)
print(' '.join(map(str, anss)))

Submission status: Accepted

Solution code:
import numpy as np
from scipy import signal

MOD = 998244353

N, M = map(int, input().split())
As = list(map(int, input().split()))
Bs = list(map(int, input().split()))

def divideArr3(As, numD):
    arrA = np.array(As, dtype=np.int64)
    mask = (1<<numD) - 1
    arrA0 = arrA & mask
    arrA1 = (arrA>>numD) & mask
    arrA2 = arrA>>(2*numD)
    return arrA2, arrA1, arrA0

def convolveArr(arrA, arrB):
    arrC = np.rint(signal.fftconvolve(arrA, arrB)).astype(np.int64)
    return arrC

def convolveMOD(As, Bs, MOD):
    numD = 10
    div = 1<<numD
    arrA2, arrA1, arrA0 = divideArr3(As, numD)
    arrB2, arrB1, arrB0 = divideArr3(Bs, numD)

    arrC4 = convolveArr(arrA2, arrB2) % MOD
    arrC2 = convolveArr(arrA1, arrB1) % MOD
    arrC0 = convolveArr(arrA0, arrB0) % MOD
    arrC3 = (convolveArr(arrA2+arrA1, arrB2+arrB1) % MOD - arrC4 - arrC2) % MOD
    arrC1 = (convolveArr(arrA1+arrA0, arrB1+arrB0) % MOD - arrC2 - arrC0) % MOD
    arrC2 += (convolveArr(arrA2+arrA0, arrB2+arrB0) % MOD - arrC4 - arrC0) % MOD
    arrC2 %= MOD

    div2 = div*div%MOD
    div3 = div2*div%MOD
    div4 = div3*div%MOD
    arrAns = ((arrC4*div4)%MOD + (arrC3*div3)%MOD + (arrC2*div2)%MOD + (arrC1*div)%MOD + arrC0) % MOD
    return arrAns.tolist()

anss = convolveMOD(As, Bs, MOD)
print(' '.join(map(str, anss)))

Submission status: Accepted

Solution code:
MOD = 998244353


def convolution(a, b):
    def butterfly(arr):
        n = len(arr)
        h = (n - 1).bit_length()
        sum_e = (911660635, 509520358, 369330050, 332049552, 983190778, 123842337, 238493703, 975955924, 603855026, 856644456, 131300601,
                 842657263, 730768835, 942482514, 806263778, 151565301, 510815449, 503497456, 743006876, 741047443, 56250497, 0, 0, 0, 0, 0, 0, 0, 0, 0)
        for ph in range(1, h + 1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            now = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = arr[i + offset]
                    r = arr[i + offset + p] * now
                    arr[i + offset] = (l + r) % MOD
                    arr[i + offset + p] = (l - r) % MOD
                now *= sum_e[(~s & -~s).bit_length() - 1]
                now %= MOD

    def butterfly_inv(arr):
        n = len(arr)
        h = (n - 1).bit_length()
        sum_ie = (86583718, 372528824, 373294451, 645684063, 112220581, 692852209, 155456985, 797128860, 90816748, 860285882, 927414960, 354738543,
                  109331171, 293255632, 535113200, 308540755, 121186627, 608385704, 438932459, 359477183, 824071951, 0, 0, 0, 0, 0, 0, 0, 0, 0)
        for ph in range(1, h + 1)[::-1]:
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            inow = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = arr[i + offset]
                    r = arr[i + offset + p]
                    arr[i + offset] = (l + r) % MOD
                    arr[i + offset + p] = (MOD + l - r) * inow % MOD
                inow *= sum_ie[(~s & -~s).bit_length() - 1]
                inow %= MOD

    n = len(a)
    m = len(b)
    if not n or not m:
        return []
    if min(n, m) <= 50:
        if n < m:
            n, m = m, n
            a, b = b, a
        res = [0] * (n + m - 1)
        for i in range(n):
            for j in range(m):
                res[i + j] += a[i] * b[j]
                res[i + j] %= MOD
        return res
    z = 1 << (n + m - 2).bit_length()
    a += [0] * (z - n)
    b += [0] * (z - m)
    butterfly(a)
    butterfly(b)
    for i in range(z):
        a[i] *= b[i]
        a[i] %= MOD
    butterfly_inv(a)
    a = a[:n + m - 1]
    iz = pow(z, MOD - 2, MOD)
    for i in range(n + m - 1):
        a[i] *= iz
        a[i] %= MOD
    return a


def atcoder_practice2_f():
    import sys
    input = sys.stdin.buffer.readline

    N, M = map(int, input().split())
    A = [int(x) for x in input().split()]
    B = [int(x) for x in input().split()]

    print(*convolution(A, B))


if __name__ == "__main__":
    atcoder_practice2_f()

Submission status: Accepted

Solution code:
MOD = 998244353

def convolution(a, b):
    def butterfly(arr):
        n = len(arr)
        h = (n - 1).bit_length()
        sum_e = (911660635, 509520358, 369330050, 332049552, 983190778, 123842337, 238493703, 975955924, 603855026, 856644456, 131300601,
                842657263, 730768835, 942482514, 806263778, 151565301, 510815449, 503497456, 743006876, 741047443, 56250497, 0, 0, 0, 0, 0, 0, 0, 0, 0)
        for ph in range(1, h + 1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            now = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = arr[i + offset]
                    r = arr[i + offset + p] * now
                    arr[i + offset] = (l + r) % MOD
                    arr[i + offset + p] = (l - r) % MOD
                now *= sum_e[(~s & -~s).bit_length() - 1]
                now %= MOD


    def butterfly_inv(arr):
        n = len(arr)
        h = (n - 1).bit_length()
        sum_ie = (86583718, 372528824, 373294451, 645684063, 112220581, 692852209, 155456985, 797128860, 90816748, 860285882, 927414960, 354738543,
                109331171, 293255632, 535113200, 308540755, 121186627, 608385704, 438932459, 359477183, 824071951, 0, 0, 0, 0, 0, 0, 0, 0, 0)
        for ph in range(1, h + 1)[::-1]:
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            inow = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = arr[i + offset]
                    r = arr[i + offset + p]
                    arr[i + offset] = (l + r) % MOD
                    arr[i + offset + p] = (MOD + l - r) * inow % MOD
                inow *= sum_ie[(~s & -~s).bit_length() - 1]
                inow %= MOD

    n = len(a)
    m = len(b)
    if not n or not m: return []
    if min(n, m) <= 50:
        if n < m:
            n, m = m, n
            a, b = b, a
        res = [0] * (n + m - 1)
        for i in range(n):
            for j in range(m):
                res[i + j] += a[i] * b[j]
                res[i + j] %= MOD
        return res
    z = 1 << (n + m - 2).bit_length()
    a += [0] * (z - n)
    b += [0] * (z - m)
    butterfly(a)
    butterfly(b)
    for i in range(z):
        a[i] *= b[i]
        a[i] %= MOD
    butterfly_inv(a)
    a = a[:n + m - 1]
    iz = pow(z, MOD - 2, MOD)
    for i in range(n + m - 1):
        a[i] *= iz
        a[i] %= MOD
    return a


def atcoder_practice2_f():
    import sys
    input = sys.stdin.buffer.readline

    N, M = map(int, input().split())
    A = [int(x) for x in input().split()]
    B = [int(x) for x in input().split()]

    print(*convolution(A, B))


if __name__ == "__main__":
    atcoder_practice2_f()

Submission status: Accepted

Solution code:

MOD = 998244353
sum_e = (911660635, 509520358, 369330050, 332049552, 983190778, 123842337, 238493703, 975955924, 603855026, 856644456, 131300601, 842657263, 730768835, 942482514, 806263778, 151565301, 510815449, 503497456, 743006876, 741047443, 56250497, 0, 0, 0, 0, 0, 0, 0, 0, 0)
sum_ie = (86583718, 372528824, 373294451, 645684063, 112220581, 692852209, 155456985, 797128860, 90816748, 860285882, 927414960, 354738543, 109331171, 293255632, 535113200, 308540755, 121186627, 608385704, 438932459, 359477183, 824071951, 0, 0, 0, 0, 0, 0, 0, 0, 0)
class Convolution:
    # ROOT = 3
    def __init__(self):
        pass
        # global MOD
        # MOD = mod
        # # self.MOD = mod

    def __butterfly(self, arr):
        n = len(arr)
        h = (n - 1).bit_length()
        # MOD = self.MOD

        for ph in range(1, h + 1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            now = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = arr[i + offset]
                    r = arr[i + offset + p] * now
                    arr[i + offset] = (l + r) % MOD
                    arr[i + offset + p] = (l - r) % MOD
                now *= sum_e[(~s & -~s).bit_length() - 1]
                now %= MOD

    def __butterfly_inv(self, arr):
        n = len(arr)
        h = (n - 1).bit_length()
        # MOD = self.MOD

        for ph in range(1, h + 1)[::-1]:
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            inow = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = arr[i + offset]
                    r = arr[i + offset + p]
                    arr[i + offset] = (l + r) % MOD
                    arr[i + offset + p] = (MOD + l - r) * inow % MOD
                inow *= sum_ie[(~s & -~s).bit_length() - 1]
                inow %= MOD

    def __call__(self, a, b):
        n = len(a)
        m = len(b)
        # MOD = self.MOD
        if not n or not m: return []
        if min(n, m) <= 50:
            if n < m:
                n, m = m, n
                a, b = b, a
            res = [0] * (n + m - 1)
            for i in range(n):
                for j in range(m):
                    res[i + j] += a[i] * b[j]
                    res[i + j] %= MOD
            return res
        z = 1 << (n + m - 2).bit_length()
        a += [0] * (z - n)
        b += [0] * (z - m)
        self.__butterfly(a)
        self.__butterfly(b)
        for i in range(z):
            a[i] *= b[i]
            a[i] %= MOD
        self.__butterfly_inv(a)
        a = a[:n + m - 1]
        iz = pow(z, MOD - 2, MOD)
        for i in range(n + m - 1):
            a[i] *= iz
            a[i] %= MOD
        return a




def atcoder_practice2_f():
    import sys
    input = sys.stdin.buffer.readline

    N, M = map(int, input().split())
    A = [int(x) for x in input().split()]
    B = [int(x) for x in input().split()]

    convolution = Convolution()
    print(*convolution(A, B))


if __name__ == "__main__":
    atcoder_practice2_f()

Submission status: Accepted

Solution code:
MOD = 998244353

def butterfly(arr):
    n = len(arr)
    h = (n - 1).bit_length()
    sum_e = (911660635, 509520358, 369330050, 332049552, 983190778, 123842337, 238493703, 975955924, 603855026, 856644456, 131300601,
                842657263, 730768835, 942482514, 806263778, 151565301, 510815449, 503497456, 743006876, 741047443, 56250497, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    for ph in range(1, h + 1):
        w = 1 << (ph - 1)
        p = 1 << (h - ph)
        now = 1
        for s in range(w):
            offset = s << (h - ph + 1)
            for i in range(p):
                l = arr[i + offset]
                r = arr[i + offset + p] * now
                arr[i + offset] = (l + r) % MOD
                arr[i + offset + p] = (l - r) % MOD
            now *= sum_e[(~s & -~s).bit_length() - 1]
            now %= MOD

def butterfly_inv(arr):
    n = len(arr)
    h = (n - 1).bit_length()
    sum_ie = (86583718, 372528824, 373294451, 645684063, 112220581, 692852209, 155456985, 797128860, 90816748, 860285882, 927414960, 354738543,
                109331171, 293255632, 535113200, 308540755, 121186627, 608385704, 438932459, 359477183, 824071951, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    for ph in range(1, h + 1)[::-1]:
        w = 1 << (ph - 1)
        p = 1 << (h - ph)
        inow = 1
        for s in range(w):
            offset = s << (h - ph + 1)
            for i in range(p):
                l = arr[i + offset]
                r = arr[i + offset + p]
                arr[i + offset] = (l + r) % MOD
                arr[i + offset + p] = (MOD + l - r) * inow % MOD
            inow *= sum_ie[(~s & -~s).bit_length() - 1]
            inow %= MOD



def convolution(a, b):
    n = len(a)
    m = len(b)
    if not n or not m: return []
    if min(n, m) <= 50:
        if n < m:
            n, m = m, n
            a, b = b, a
        res = [0] * (n + m - 1)
        for i in range(n):
            for j in range(m):
                res[i + j] += a[i] * b[j]
                res[i + j] %= MOD
        return res
    z = 1 << (n + m - 2).bit_length()
    a += [0] * (z - n)
    b += [0] * (z - m)
    butterfly(a)
    butterfly(b)
    for i in range(z):
        a[i] *= b[i]
        a[i] %= MOD
    butterfly_inv(a)
    a = a[:n + m - 1]
    iz = pow(z, MOD - 2, MOD)
    for i in range(n + m - 1):
        a[i] *= iz
        a[i] %= MOD
    return a


def atcoder_practice2_f():
    import sys
    input = sys.stdin.buffer.readline

    N, M = map(int, input().split())
    A = [int(x) for x in input().split()]
    B = [int(x) for x in input().split()]

    print(*convolution(A, B))


if __name__ == "__main__":
    atcoder_practice2_f()