Problem:
 Score : 900 points    Problem Statement We have a sequence of N integers: x=(x_0,x_1,\cdots,x_{N-1}) .
Initially, x_i=0 for each i ( 0 \leq i \leq N-1 ).  Snuke will perform the following operation exactly  M times:   Choose two distinct indices i, j ( 0 \leq i,j \leq N-1,\ i \neq j ).
Then, replace x_i with x_i+2 and x_j with x_j+1 .   Find the number of different sequences that can result after M operations.
Since it can be enormous, compute the count modulo 998244353 .      Constraints  2 \leq N \leq 10^6  1 \leq M \leq 5 \times 10^5  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M       Output Print the number of different sequences that can result after M operations, modulo 998244353 .        Sample Input 1 2 2      Sample Output 1 3  After two operations, there are three possible outcomes:   x=(2,4)  x=(3,3)  x=(4,2)   For example, x=(3,3) can result after the following sequence of operations:   First, choose i=0,j=1 , changing x from (0,0) to (2,1) .  Second, choose i=1,j=0 , changing x from (2,1) to (3,3) .        Sample Input 2 3 2      Sample Output 2 19       Sample Input 3 10 10      Sample Output 3 211428932       Sample Input 4 100000 50000      Sample Output 4 3463133

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

U = 4*10**6
MOD = 998244353
 
fact = [1]*(U+1)
fact_inv = [1]*(U+1)
 
for i in range(1,U+1):
    fact[i] = (fact[i-1]*i)%MOD
fact_inv[U] = pow(fact[U], MOD-2, MOD)
 
for i in range(U,0,-1):
    fact_inv[i-1] = (fact_inv[i]*i)%MOD
 
def perm(n, k):
    if k < 0 or k > n:
        return 0
    z = fact[n]
    z *= fact_inv[n-k]
    z %= MOD
    return z

def comb(n, k):
    if k < 0 or k > n:
        return 0
    z = fact[n]
    z *= fact_inv[k]
    z %= MOD
    z *= fact_inv[n-k]
    z %= MOD
    return z

n, m = map(int, input().split())
ans = 0
for k in range(m%2, min(n, m)+1, 2):
  ans += comb((3*m-k)//2+n-1, n-1) * comb(n, k)
  ans %= MOD
cnt = 0
for k in range(m):
  cnt += comb(k+n-2, n-2)
  cnt %= MOD
cnt *= n
cnt %= MOD
ans -= cnt
ans %= MOD
print(ans)

Solution (Python):
class Combination:
    """階乗とその逆元のテーブルをO(N)で事前作成し、組み合わせの計算をO(1)で行う"""
    def __init__(self, n, MOD):
        self.fact = [1]
        for i in range(1, n + 1):
            self.fact.append(self.fact[-1] * i % MOD)
        self.inv_fact = [0] * (n + 1)
        self.inv_fact[n] = pow(self.fact[n], MOD - 2, MOD)
        for i in reversed(range(n)):
            self.inv_fact[i] = self.inv_fact[i + 1] * (i + 1) % MOD
        self.MOD = MOD

    def inverse(self, k):
        """kの逆元を求める O(1)"""
        return (self.inv_fact[k] * self.fact[k - 1]) % self.MOD

    def factorial(self, k):
        """k!を求める O(1)"""
        return self.fact[k]

    def inverse_factorial(self, k):
        """k!の逆元を求める O(1)"""
        return self.inv_fact[k]

    def permutation(self, k, r):
        """kPrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r]) % self.MOD

    def combination(self, k, r):
        """kCrを求める O(1)"""
        if k < r:
            return 0
        return (self.fact[k] * self.inv_fact[k - r] * self.inv_fact[r]) % self.MOD


n, m = map(int, input().split())
MOD = 998244353
comb = Combination(10 ** 7, MOD)


ans = 0
for odd_cnt in range(m + 1):
    if odd_cnt > n:
        continue
    if (3 * m - odd_cnt) % 2 == 1:
        continue
    ptn = 1
    box = n
    odd_ball = odd_cnt
    ptn *= comb.combination(box, odd_ball)

    even_pair_ball = (3 * m - odd_cnt) // 2
    ptn *= comb.combination(even_pair_ball + box - 1, even_pair_ball)

    ans += ptn
    ptn %= MOD

for i in range(1, m + 1):
    ball = m - i
    box = n - 1
    ans -= comb.combination(ball + box - 1, ball) * n
    ans %= MOD

print(ans)

Solution (Python):
n,m=map(int,input().split())
o,f,i=998244353,[1],1
while i<3*m+n:f+=[f[-1]*i%o];i+=1
c=lambda x,y=n-1:f[x]*pow(f[x-y]*f[y],o-2,o)
a=c(-1)-n*c(n+m-2)
while~-n>m<i-n:m+=2;a-=c(n,m)*c(n-1+i-m>>1)
print(a%o)

Solution (Python):
n,m = map(int,input().split())
mod = 998244353





max = n+2*m + 100
fac, finv, inv = [0]*max, [0]*max, [0]*max

def comInit(max):
    fac[0] = fac[1] = 1
    finv[0] = finv[1] = 1
    inv[1] = 1

    for i in range(2,max):
      fac[i] = fac[i-1]* i% mod
      inv[i] = mod - inv[mod%i] * (mod // i) % mod
      finv[i] = finv[i-1] * inv[i] % mod

comInit(max)


def com(n,k):
    if(n < k):
        return 0
    if( (n<0) | (k < 0)):
        return 0
    return fac[n] * (finv[k] * finv[n-k] % mod) % mod


a = 0
for x in range(min(m,n)+1):
    if(3*m-x)%2==1:
        continue
    y = (3*m-x)//2
    a += com(n,x)*fac[y+n-1]*finv[y]*finv[n-1]
    a %= mod

b = fac[n-1+m-1] * finv[n-1] * finv[m-1] * n
b %= mod

ans = a-b
ans %= mod

print(ans)

Solution (Python):
N, M = map(int, input().split())

mod = 998244353
def cmb(n, r, mod=mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

NN = 4*10**6 
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 

for i in range( 2, NN + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )


A = []
tmp = 0
for i in range(2*M+N+2):
    tmp = (tmp + cmb(N-2+i, N-2)) % mod
    A.append(tmp)

ans = 0
for k in range(min(N, M)+1):
    if (M-k) % 2 != 0: continue
    l = (3*M-k)//2
    a = cmb(N-1+l, N-1) * cmb(N, k)
    b = N * (cmb(N-1, k-1)*A[l-M])
    if l >= M+1:
        b += N * cmb(N-1, k)*A[l-M-1]
    
    ans = (ans + a - b) % mod 
print(ans)

Solution (Python):
mod = 998244353
eps = 10**-9


def main():
    import sys
    input = sys.stdin.readline

    
    
    nmax = 25 * 10 ** 5 + 10  
    fac = [0] * nmax
    finv = [0] * nmax
    inv = [0] * nmax
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    finv[1] = 1
    inv[1] = 1
    for i in range(2, nmax):
        fac[i] = fac[i - 1] * i % mod
        inv[i] = mod - inv[mod % i] * (mod // i) % mod
        finv[i] = finv[i - 1] * inv[i] % mod

    def comb(n, r):
        if n < r:
            return 0
        else:
            return (fac[n] * ((finv[r] * finv[n - r]) % mod)) % mod

    N, M = map(int, input().split())
    ans = comb(3 * M + N - 1, N-1)
    for i in range(2 * M + 1, 3 * M+1):
        ans = (ans - (N * comb(3*M - i + N-2, N-2))%mod)%mod
    for i in range(M+2, 3 * M + 1, 2):
        if i > N:
            break
        j = (3 * M - i) // 2
        ans = (ans - (comb(N, i) * comb(j + N-1, N-1))%mod)%mod
    print(ans)


if __name__ == '__main__':
    main()

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
from pprint import pprint
from copy import deepcopy
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor
from operator import mul
from functools import reduce
from pprint import pprint


sys.setrecursionlimit(2147483647)
INF = 10 ** 15
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 998244353

n, m = LI()
a = n + 3 * m
fac = [1] * (a + 1)
inv = [1] * (a + 1)
for j in range(1, a + 1):
    fac[j] = fac[j-1] * j % mod

inv[a] = pow(fac[a], mod-2, mod)
for j in range(a-1, -1, -1):
    inv[j] = inv[j+1] * (j+1) % mod

def comb(x, y):
    if y > x or x < 0 or y < 0:
        return 0
    return fac[x] * inv[x - y] * inv[y] % mod

ans = 0
for i in range(min(m + 1, n + 1)):
    if (m - i) % 2:
        continue
    ans += comb(n, i) * comb((m * 3 - i) // 2 + n - 1, n - 1)
    ans %= mod

print((ans - n * comb(m - 1 + n - 1, n - 1)) % mod)

Solution (Python):
n,m=map(int,input().split())
o,f,i=998244353,[1],1
while i<3*m+n:f+=[f[-1]*i%o];i+=1
c=lambda x,y=n-1:f[x]*pow(f[x-y]*f[y]%o,o-2,o)
a=c(-1)-n*c(n+m-2)
while~-n>m<i-n:m+=2;a-=c(n,m)*c(n-1+i-m>>1)
print(a%o)

Solution (Python):
import sys
input=sys.stdin.readline
N,M=map(int,input().split())
mod=998244353

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


g1=[1,1]
g2=[1,1]
for i in range(2,3*M//2+N+2):
    num_1=g1[-1]*i%mod
    g1.append(num_1)
    g2.append(g2[-1]*inv[i]%mod)
    
def cmb(n,r):
    return g1[n]*g2[r]*g2[n-r]%mod




ans=0
for i in range(min(N+1,M+1)):
    if M%2!=i%2:
        continue
    m=(3*M-i)//2
    num=(cmb(m+N-1,N-1)-i*cmb(m-M+N-1,N-1))%mod
    if m>=M+1:
        num=(num-(N-i)*cmb(m-M+N-2,N-1))%mod
    ans=(ans+cmb(N,i)*num)%mod

print(ans)

Solution (Python):
n,m=map(int,input().split())
w,u=3*m,n-1
o,f,i=998244353,[1],1
while i<w+n:f+=[f[-1]*i%o];i+=1
c=lambda n,k=u:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
a=c(w+u)-n*c(n+m-2)
while~-n>m<w:m+=2;a-=c(n,m)*c(2*u+w-m>>1)
print(a%o)

Solution (Python):
n,m=map(int,input().split())
w,u=3*m,n-1
o,f,i=998244353,[1],1
while i<w+n:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
a=c(w+u,u)-n*c(n+m-2,u)
while~-n>m<w:m+=2;a-=c(n,m)*c(2*u+w-m>>1,u)
print(a%o)

Solution (Python):
n,m=map(int,input().split())
w,u=3*m,n-1
o,f,i=998244353,[1],1
while i<w+n:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
a=c(w+u,u)-n*c(n+m-2,u)
while n>-~m<w:m+=2;a-=c(n,m)*c(2*u+w-m>>1,u)
print(a%o)

Solution (Python):
n,m=map(int,input().split())
w,u=3*m,n-1
o,f,i=998244353,[1],1
while i<w+n:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
a=c(w+u,u)-n*c(n+m-2,u)
while-~m<min(w,n):m+=2;a-=c(n,m)*c(2*u+w-m>>1,u)
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<o>>8:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split());w,u=3*m,n-1;a=c(w+u,u)-n*c(n+m-2,u)
while-~m<min(w,n):m+=2;a-=c(n,m)*c(2*u+w-m>>1,u)
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<o>>8:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
while-~m<min(w,n):m+=2;a-=c(n,m)*c(2*u+w-m>>1,u)
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<o/399:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
while-~m<min(w,n):m+=2;a-=c(n,m)*c((w-m)//2+u,u)
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<o/399:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
while-~m<min(w,n):m+=2;a-=c(n,m)*c((w-m)//2+u,u)
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<o/399:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
m+=1
if w-m&1:m+=1
while m<=min(w,n):a-=c(n,m)*c((w-m)//2+u,u);m+=2
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<o/399:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u,j=3*m,n-1,m+1
a=c(w+u,u)-n*c(n+m-2,u)
if w-j&1:j+=1
while j<=min(w,n):a-=c(n,j)*c((w-j)//2+u,u);j+=2
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<o/399:f+=[f[-1]*i%o];i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
if u>m:
 for i in range(m+1,min(w,n)+1):
  if w-i+1&1:a-=c(n,i)*c((w-i)//2+u,u) 
print(a%o)

Solution (Python):
o,f,i=998244353,[1],1
while i<2500100:f.append(f[-1]*i%o);i+=1
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
if u>m:
 for i in range(m+1,min(w,n)+1):
  if w-i+1&1:a-=c(n,i)*c((w-i)//2+u,u) 
print(a%o)

Solution (Python):
o,f=998244353,[1]
for i in range(1,2500100):f.append(f[-1]*i%o)
c=lambda n,k:f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
if u>m:
 for i in range(m+1,min(w,n)+1):
  if w-i+1&1:a-=c(n,i)*c((w-i)//2+u,u) 
print(a%o)

Solution (Python):
o,f=998244353,[1]
for i in range(1,2500100):f.append(f[-1]*i%o)
def c(n,k):return f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w,u=3*m,n-1
a=c(w+u,u)-n*c(n+m-2,u)
if u>m:
 for i in range(m+1,min(w,n)+1):
  if w-i+1&1:a-=c(n,i)*c((w-i)//2+u,u)
print(a%o)

Solution (Python):
o=998244353
f=[1]
for i in range(1,2500100):f.append(f[-1]*i%o)
def c(n,k):return f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w=3*m
u=n-1
a=c(w+u,u)-n*c(n+m-2,u)
if u>m:
 for i in range(m+1,min(w,n)+1):
  if w-i+1&1:a-=c(n,i)*c((w-i)//2+u,u)
print(a%o)

Solution (Python):
o=998244353
f=[1]
for i in range(1,2500100):f.append(f[-1]*i%o)
def c(n,k):return f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)
n,m=map(int,input().split())
w=3*m
u=n-1
a=c(w+u,u)-n*c(n+m-2,u)
if u>m:
 for i in range(m+1,min(w,n)+1):
  if w-i+1&1:a=(a-c(n,i)*c((w-i)//2+u,u))%o
print(a%o)

Solution (Python):
o=998244353
f=[1]
for i in range(1,2500100):f.append(f[-1]*i%o)
def c(n,k):return f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)%o
n,m=map(int,input().split())
w=3*m
u=n-1
a=c(w+u,u)-n*c(n+m-2,u)
if u>m:
 for i in range(m+1,min(w,n)+1):
  if w-i+1&1:a=(a-c(n,i)*c((w-i)//2+u,u))%o
print(a%o)

Solution (Python):
o=998244353
r=2500100
f=[1]
for i in range(1,r):f.append(f[-1]*i%o)
def c(n,k):return f[n]*pow(f[n-k],o-2,o)*pow(f[k],o-2,o)%o
n,m=map(int,input().split())
w=3*m
a=(c(n+w-1,n-1)-n*c(n+m-2,n-1))%o
if n>m+1:
 for i in range(m+1,min(w,n)+1):
  if w-i&1:continue
  a=(a-c(n,i)*c((w-i)//2+n-1,n-1))%o
print(a)

Solution (Python):
o=998244353
r=2500100
f=[1]
for i in range(1,r):f.append(f[-1]*i%o)
def v(x):return pow(x,o-2,o)
def c(n,k):return f[n]*v(f[n-k])*v(f[k])%o
n,m=map(int,input().split())
w=3*m
a=(c(n+w-1,n-1)-n*c(n+m-2,n-1))%o
if n>m+1:
 for i in range(m+1,min(w,n)+1):
  if (w-i)%2:continue
  x=(w-i)//2
  a=(a-c(n,i)*c(x+n-1,n-1))%o
print(a)

Solution (Python):
o=998244353
r=2500100
f=[1]
for i in range(1,r):
 f.append(f[-1]*i%o)
def v(x):
 return pow(x,o-2,o)
def c(n,k):
 return f[n]*v(f[n-k])*v(f[k])%o
n,m=map(int,input().split())
w=3*m
a=(c(n+w-1,n-1)-n*c(n+m-2,n-1))%o
if n>m+1:
 for i in range(m+1,min(w,n)+1):
  if (w-i)%2:continue
  x=(w-i)//2
  a=(a-c(n,i)*c(x+n-1,n-1))%o
print(a)

Solution (Python):
o = 998244353
r = 2500100
f = [1]
for i in range(1,r):
 f.append(f[-1]*i%o)
def v(x):
 return pow(x,o-2,o)
def c(n,k):
 return f[n]*v(f[n-k])*v(f[k])%o
n,m = map(int,input().split())
a = (c(n+3*m-1,n-1)-n*c(n+m-2,n-1))%o
if n>m+1:
 for i in range(m+1,min(3*m,n)+1):
  if (3*m-i)%2:continue
  x = (3*m-i)//2
  a = (a-c(n,i)*c(x+n-1,n-1))%o
print(a)

Solution (Python):
mod = 998244353
rng = 2500100
fctr = [1]
finv = [1]
for i in range(1,rng):
  fctr.append(fctr[-1]*i%mod)
def inv(x):
  return pow(x,mod-2,mod)
def cmb(n,k):
  if n<0 or k<0:
    return 0
  else:
    return fctr[n]*inv(fctr[n-k])*inv(fctr[k])%mod

n,m = map(int,input().split())
ans = cmb(n+3*m-1,n-1)
ans = (ans-n*cmb(n+m-2,n-1))%mod
if n>m+1:
  for i in range(m+1,min(3*m,n)+1):
    if (3*m-i)%2 == 1:
      continue
    x = (3*m-i)//2
    ans = (ans-cmb(n,i)*cmb(x+n-1,n-1))%mod
print(ans)

Solution (Python):

"""
Writer: SPD_9X2
https://atcoder.jp/contests/agc036/tasks/agc036_c

可能性:
3MをN個に分割する分け方から、存在しない分け方を引く

存在しない分け方とは？
→奇数個だけ置かれている場所の数がM個より多い場合

k(M+1～N)までに関して、1をk箇所に分配し、その後2個ずつ好きな場所に分配する方法を数えればいい
→解けたかな？

→まだ駄目な条件が抜けてる
最大値が2Mより大きいならダメ
↑の部分で引いてないのは、奇数個置かれている場所はM個以下だが、最大値が2Mより大きい分け方…

(k個の1のおき方) * {(1の中から2Mを超えるように置く場所の選び方)*残りの2の分配法}
+ (k個の1のおき方) * {(0の中から2Mを超えるように置く場所の選び方)*残りの2の分配法}

をk(0以上M以下に関して求める)

→解けた！
"""

def modfac(n, MOD):
 
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv
    return factorials, invs


def modnCr(n,r,mod,fac,inv): 
    if n < r:
        return 0
    return fac[n] * inv[n-r] * inv[r] % mod

mod = 998244353
N,M = map(int,input().split())

fac,inv = modfac(N+3*M+10,mod)

ans = modnCr(3*M+N-1,N-1,mod,fac,inv)


for k in range(M+1,N+1):

    

    if (3*M-k) % 2 == 0:
        ans -= modnCr(N,k,mod,fac,inv) * modnCr((3*M-k)//2 + N-1 , N-1 , mod,fac,inv)
        ans %= mod

for k in range(M+1):

    if (3*M-k) % 2 == 0:
        
        amari = (3*M - (k+2*M))//2

        now = modnCr(N,k,mod,fac,inv) * k * modnCr(amari+N-1,N-1, mod,fac,inv)

        amari -= 1
        now +=modnCr(N,k,mod,fac,inv)*(N-k)*modnCr(amari+N-1,N-1, mod,fac,inv)

        ans -= now
        ans %= mod

    

print (ans)

Solution (Python):
from collections import defaultdict
class Combination():
    def __init__(self, n, mod):
        self.n = n
        self.mod = mod
        self.fact = self.make_fact(n)
        self.fact_inv = self.make_fact_inv(n)
    def make_fact(self, n):
        res = [1]*(n+1)
        for i in range(1, n+1):
            res[i] = res[i-1]*i%self.mod
        return res
    def make_fact_inv(self, n):
        fact_inv = [1]*(n+1)
        fact_inv[n] = pow(self.fact[n], self.mod-2, self.mod)
        for i in range(n, 0, -1):
            fact_inv[i-1] = fact_inv[i]*i%self.mod
        return fact_inv
    def comb(self, m, k):
        if m<0 or k<0 or m<k:
            return 0
        return self.fact[m]*self.fact_inv[k]*self.fact_inv[m-k]%self.mod
mod = 998244353
n, m = map(int, input().split())
C = Combination(n+m*3, mod)
ans = 0
for i in range(min(n, m)+1):
    if m%2 != i%2:
        continue
    res = 0
    r = (m*3-i)//2
    res+=C.comb(r+n-1, r)
    res-=C.comb(r-m-1+n-1, n-1)*(n-i)
    res-=C.comb(r-m+n-1, n-1)*i
    res*=C.comb(n, i)
    res%-mod
    ans+=res
    ans%=mod
print(ans)

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
from bisect import bisect_left, bisect_right
import random
from itertools import permutations, accumulate, combinations
import sys
import string



INF = float('inf')
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS(): return sys.stdin.readline().split()
def S(): return sys.stdin.readline().strip()
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 998244353


n, m = LI()
total = m * 3

fac = [1] * (total+n+1)
inv = [1] * (total+n+1)
for i in range(total+n):
    fac[i+1] = fac[i]*(i+1)%mod

inv[total+n]=pow(fac[-1], mod-2, mod)
for j in range(total+n-1, -1, -1):
    inv[j]=inv[j+1]*(j+1)%mod



def comb(n, r):
    if r > n:
        return 0
    return fac[n] * inv[n - r] * inv[r] % mod


ans = comb(total+n-1, n-1)
for i in range(m + 2, min(n + 1, total + 1)):
    if (total - i) % 2 == 0:
        ans -= comb(n, i) * comb(n + (total - i) // 2 - 1, n - 1) % mod
        ans %= mod


ret = 0
for i in range(m):
    ret = (ret + comb(i + n - 2, n - 2)) % mod


ans -= (ret * n) % mod
print(ans % mod)

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

import numpy as np

N,M = map(int,input().split())
MOD = 998244353

def cumprod_mod(arr):
    L = len(arr)
    Lsq = int(L**.5+1)
    arr = np.resize(arr,Lsq**2).reshape(Lsq,Lsq)
    for n in range(1,Lsq):
        arr[:,n] *= arr[:,n-1]; arr[:,n] %= MOD
    for n in range(1,Lsq):
        arr[n] *= arr[n-1,-1]; arr[n] %= MOD
    arr = arr.ravel()[:L]
    return arr

U = 5*10**6
x = np.arange(U,dtype=np.int64)
x[0] = 1
fact = cumprod_mod(x)

x = np.arange(U,0,-1,dtype=np.int64)
x[0] = pow(int(fact[-1]),MOD-2,MOD)
fact_inv = cumprod_mod(x)[::-1]

def comb(n,k):
    return fact[n] * fact_inv[k] % MOD * fact_inv[n-k] % MOD

x = comb(3*M+N-1,N-1) 
y = N * comb(M+N-2,N-1) 


odd_cnt = np.arange(M+2,3*M+1,2)
S = (3*M - odd_cnt)//2
bl = (odd_cnt <= N)&(S>=0)
odd_cnt = odd_cnt[bl]
S = S[bl]
arr = fact[N] * fact_inv[odd_cnt] % MOD * fact_inv[N-odd_cnt] % MOD

arr *= fact[S+N-1] * fact_inv[N-1] % MOD * fact_inv[S] % MOD
arr %= MOD
z = arr.sum() % MOD

answer = x-y-z
answer %= MOD
print(answer)

Solution (Python):
N, M = map(int, input().split())
MOD = 998244353

if N < M :
    one = N - (N - M) % 2
    two = (3 * M - one) // 2
else :
    one, two = M, M
    
L = M + M // 2 + N - 1

fac = [0] * (L + 1)
inv = [0] * (L + 1)
fac[0], inv[0] = 1, 1

for i in range(1, L + 1) :
    fac[i] = fac[i-1] * i % MOD
    inv[i] = pow(fac[i], MOD-2, MOD)

def comb(n, r) :
    if n <= 0 or r < 0 :
        return 0
    return fac[n]*inv[n-r]*inv[r]%MOD

ret = 0
while one >= 0 :
    ret += comb(N, one) * comb(two+N-1, N-1) % MOD
    ret %= MOD
    one -= 2
    two += 1
 
ret -= comb(3*M+N-1-(2*M+1), N-1) * N
print(ret % MOD)

Solution (Python):
def prepare(n, MOD):
    f = 1
    factorials = [1] * (n + 1)
    for m in range(1, n + 1):
        f = f * m % MOD
        factorials[m] = f
    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv = inv * m % MOD
        invs[m - 1] = inv

    return factorials, invs


def get_nhr(n, facts, invs, MOD):
    def nhr(r):
        return facts[n + r - 1] * invs[r] * invs[n - 1] % MOD

    return nhr


def solve(n, m):
    MOD = 998244353
    facts, invs = prepare(n + 3 * m, MOD)
    fn = facts[n]
    nhr = get_nhr(n, facts, invs, MOD)
    nm2 = min(n, m - 1)
    ans = 0
    for odd in range(m % 2, min(m, n) + 1, 2):
        tmp = nhr((3 * m - odd) // 2)
        if odd > 0:
            tmp -= odd * nhr((m - odd) // 2)
        if odd < nm2:
            tmp -= (n - odd) * nhr((m - odd - 2) // 2)
        ans = (ans + fn * invs[odd] * invs[n - odd] % MOD * tmp) % MOD
    return ans


n, m = list(map(int, input().split()))
print(solve(n, m))

Solution (Python):
from functools import lru_cache


def prepare(n, MOD):
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv

    return factorials, invs


def get_nhr(n, facts, invs):
    @lru_cache(maxsize=None)
    def nhr(r):
        return facts[n + r - 1] * invs[r] * invs[n - 1] % MOD

    return nhr


MOD = 998244353
n, m = list(map(int, input().split()))
facts, invs = prepare(n + 3 * m, MOD)
nhr = get_nhr(n, facts, invs)
ans = 0
for odd in range(m % 2, min(m, n) + 1, 2):
    tmp = nhr((3 * m - odd) // 2)
    if odd > 0:
        tmp -= odd * nhr((m - odd) // 2)
    if odd < n and odd <= m - 2:
        tmp -= (n - odd) * nhr((m - odd - 2) // 2)
    ans = (ans + facts[n] * invs[odd] * invs[n - odd] % MOD * tmp) % MOD
print(ans)

Solution (Python):
N, M = map(int, input().split())

def cmb(n, r, mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

mod = 998244353
n = 10**6*3
g1 = [1, 1]
g2 = [1, 1] 
inverse = [0, 1] 

for i in range( 2, n + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )

ans = cmb(M*3+N-1, N-1, mod)
for i in range(M*2+1, M*3+1):
    ans -= cmb(M*3-i+(N-1)-1, N-2, mod) * N - mod*mod
    ans = (ans+mod)%mod

for i in range(M+1, min(N, 3*M)+1):
    if (3*M - i) % 2 == 0:
        two = (3*M - i) // 2
        ans -= cmb(two+N-1, N-1, mod) * cmb(N, i, mod) - mod*mod
        ans = (ans+mod)%mod
print(ans)

Solution (Python):
import os
import sys

import numpy as np

if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")

sys.setrecursionlimit(2147483647)
INF = float("inf")
IINF = 10 ** 18
MOD = 998244353

N, M = list(map(int, sys.stdin.readline().split()))


def get_factorials(max, mod=None):
    """
    階乗 0!, 1!, 2!, ..., max!
    :param int max:
    :param int mod:
    :return:
    """
    ret = np.ones(max + 1, dtype=int)
    n = 1
    for i in range(1, max + 1):
        n *= i
        n %= mod
        ret[i] = n
    return ret


def mod_invs(max, mod):
    """
    逆元のリスト 0 から max まで
    :param max:
    :param mod:
    :return:
    """
    invs = [1] * (max + 1)
    for x in range(2, max + 1):
        invs[x] = (-(mod // x) * invs[mod % x]) % mod
    return invs


factorials = get_factorials(M * 3 // 2 + N, MOD)
finvs = np.ones(M * 3 // 2 + N + 1, dtype=int)
p = 1
for i, inv in enumerate(mod_invs(M * 3 // 2 + N, MOD)):
    p = p * inv % MOD
    finvs[i] = p


def ncr(n, r, mod=None):
    """
    scipy.misc.comb または scipy.special.comb と同じ
    組み合わせの数 nCr
    """
    ret = factorials[n] * finvs[r] % MOD * finvs[n - r] % mod
    ret[n < r] = 0
    return ret



odds = np.arange(M % 2, M + 1, 2)
halves1 = (M * 3 - odds) // 2

c1 = ncr(halves1 + N - 1, N - 1, MOD) * ncr(N, odds, MOD) % MOD




halves2 = (M - odds) // 2

c2a = ncr(halves2 + N - 1, N - 1, MOD) * ncr(N, odds, MOD) % MOD

c2b = ncr(halves2 + N - 2, N - 2, MOD) * ncr(N - 1, odds, MOD) % MOD

c2 = (c2a - c2b) * N % MOD

print((c1 - c2).sum() % MOD)

Solution (Python):
import os
import sys

import numpy as np

if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")


def main():
    sys.setrecursionlimit(2147483647)
    INF = float("inf")
    IINF = 10 ** 18
    MOD = 998244353

    N, M = list(map(int, sys.stdin.readline().split()))

    def get_factorials(max, mod=None):
        """
        階乗 0!, 1!, 2!, ..., max!
        :param int max:
        :param int mod:
        :return:
        """
        ret = [1]
        n = 1
        if mod:
            for i in range(1, max + 1):
                n *= i
                n %= mod
                ret.append(n)
        else:
            for i in range(1, max + 1):
                n *= i
                ret.append(n)
        return ret

    def mod_invs(max, mod):
        """
        逆元のリスト 0 から max まで
        :param max:
        :param mod:
        :return:
        """
        invs = [1] * (max + 1)
        for x in range(2, max + 1):
            invs[x] = (-(mod // x) * invs[mod % x]) % mod
        return invs

    factorials = np.array(get_factorials(M * 3 // 2 + N, MOD), dtype=int)
    finvs = np.ones(M * 3 // 2 + N + 1, dtype=int)
    p = 1
    for i, inv in enumerate(mod_invs(M * 3 // 2 + N, MOD)):
        p = p * inv % MOD
        finvs[i] = p

    def ncr(n, r, mod=None):
        """
        scipy.misc.comb または scipy.special.comb と同じ
        組み合わせの数 nCr
        :param int n:
        :param int r:
        :param int mod: 3 以上の素数であること
        :rtype: int
        """
        ret = factorials[n] * finvs[r] % MOD * finvs[n - r] % mod
        ret[n < r] = 0
        return ret

    
    odds = np.arange(M % 2, M + 1, 2)
    halves = (M * 3 - odds) // 2

    oc1 = ncr(N, odds, MOD)
    oc2 = ncr(N - 1, odds, MOD)

    
    c1 = ncr(halves + N - 1, N - 1, MOD) * oc1 % MOD

    
    
    
    halves = (M - odds) // 2
    
    c2a = ncr(halves + N - 1, N - 1, MOD) * oc1 % MOD
    
    
    c2b = ncr(halves + N - 2, N - 2, MOD) * oc2 % MOD

    c2 = (c2a - c2b) * N % MOD
    print((c1 - c2).sum() % MOD)


if __name__ == '__main__':
    main()

Solution (Python):
import os
import sys

if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")

sys.setrecursionlimit(2147483647)
INF = float("inf")
IINF = 10 ** 18
MOD = 998244353

N, M = list(map(int, sys.stdin.readline().split()))


def get_factorials(max, mod=None):
    """
    階乗 0!, 1!, 2!, ..., max!
    :param int max:
    :param int mod:
    :return:
    """
    ret = [1]
    n = 1
    if mod:
        for i in range(1, max + 1):
            n *= i
            n %= mod
            ret.append(n)
    else:
        for i in range(1, max + 1):
            n *= i
            ret.append(n)
    return ret


def mod_invs(max, mod):
    """
    逆元のリスト 0 から max まで
    :param max:
    :param mod:
    :return:
    """
    invs = [1] * (max + 1)
    for x in range(2, max + 1):
        invs[x] = (-(mod // x) * invs[mod % x]) % mod
    return invs


factorials = get_factorials(M * 3 // 2 + N, MOD)
finvs = []
inv = 1
for i in mod_invs(M * 3 // 2 + N, MOD):
    inv = inv * i % MOD
    finvs.append(inv)


def ncr(n, r, mod=None):
    """
    scipy.misc.comb または scipy.special.comb と同じ
    組み合わせの数 nCr
    :param int n:
    :param int r:
    :param int mod: 3 以上の素数であること
    :rtype: int
    """
    if n < r:
        return 0

    
    return factorials[n] * finvs[r] * finvs[n - r] % mod





c1 = 0
for odds in range(M * 3 % 2, M + 1, 2):
    half = (M * 3 - odds) // 2
    
    c1 += ncr(half + N - 1, N - 1, MOD) * ncr(N, odds, MOD)
    c1 %= MOD

















c2a = 0

c2b = 0
for odds in range(M % 2, M + 1, 2):
    half = (M - odds) // 2
    c2a += ncr(half + N - 1, N - 1, MOD) * ncr(N, odds, MOD)
    c2a %= MOD
    c2b += ncr(half + N - 2, N - 2, MOD) * ncr(N - 1, odds, MOD)
    c2b %= MOD

c2 = (c2a - c2b) * N
print((c1 - c2) % MOD)

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def I(): return int(sys.stdin.readline())
def LS():return [list(x) for x in sys.stdin.readline().split()]
def S(): return list(sys.stdin.readline())[:-1]
def IR(n):
    return [I() for i in range(n)]
def LIR(n):
    return [LI() for i in range(n)]
def SR(n):
    return [S() for i in range(n)]
def LSR(n):
    return [LS() for i in range(n)]

sys.setrecursionlimit(1000000)


def A():
    n = I()
    if n <= 1000000000:
        print(0,0,1,0,0,n)
    else:
        b = 1000000000
        a = 1
        y = -n%1000000000
        x = (n+a*y)//b
        print(0,0,x,y,a,b)
        return

    return


def B():
    n,k = LI()
    a = LI()
    s = []
    f = [1]*1000000
    for i in range(k):
        for j in range(n):
            aj = a[j]
            if f[aj]:
                s.append(aj)
                f[aj] = 0
            else:
                while 1:
                    x = s.pop(-1)
                    f[x] = 1
                    if x == aj:
                        break
        print(i+1,s)
    return


mod = 998244353
def C():
    def comb(a,b):
        if a < b:
            return 1
        return fact[a]*inv[b]*inv[a-b]%mod
    n,m = LI()
    MA = 3*m+max(n,2*m)-1
    fact = [1]*(MA+1)
    for i in range(MA):
        fact[i+1] = fact[i]*(i+1)%mod
    inv = [1]*(MA+1)
    inv[MA] = pow(fact[MA],mod-2,mod)
    for i in range(MA)[::-1]:
        inv[i] = inv[i+1]*(i+1)%mod
    ans = 0
    for k in range(min(m,n)+1):
        if (3*m-k)%2 == 0:
            ans += comb(n,k)*comb((3*m-k)//2+n-1,n-1) 
            ans %= mod
    ans -= n*comb(3*m-(2*m+1)+n-1,n-1)
    ans %= mod
    print(ans)

    return


def D():
    return


def E():
    n = I()

    return


def F():
    n = I()

    return


if __name__ == "__main__":
    C()