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):

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 = math.ceil(n**0.5)
        if b**2 == n:
            print(0,0,b,0,0,n//b)
            return
        x = b
        m = b*x-n
        if m < 0:
            b += 1
            x += 1
            m = b*x-n
        n = m
        for a in range(int(n**0.5)+1)[::-1]:
            if n%a == 0:
                print(0,0,x,n//a,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()

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

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

FACT=[1]
for i in range(1,4*10**6+1):
    FACT.append(FACT[-1]*i%mod)

FACT_INV=[pow(FACT[-1],mod-2,mod)]
for i in range(4*10**6,0,-1):
    FACT_INV.append(FACT_INV[-1]*i%mod)

FACT_INV.reverse()

def Combi(a,b):
    if 0<=b<=a:
        return FACT[a]*FACT_INV[b]*FACT_INV[a-b]%mod
    else:
        return 0


ANS=Combi(M*3+N-1,N-1)

for i in range(M*2+1,M*3+1):
    k=M*3-i
    ANS=(ANS-N*Combi(k+(N-1)-1,N-2))%mod


for i in range(M+1,N+1):
    if (M*3-i)%2==1:
        continue
    k=(M*3-i)//2
    ANS=(ANS-Combi(N,i)*Combi(k+N-1,N-1))%mod

print(ANS)

Solution (Python):
MOD = 998244353
list_size = 2500100

f_list = [1] * list_size
f_r_list = [1] * list_size

for i in range(list_size - 1):
	f_list[i + 1] = (f_list[i] * (i + 2)) % MOD

f_r_list[-1] = pow(f_list[-1], MOD - 2, MOD)

for i in range(2, list_size + 1):
	f_r_list[-i] = (f_r_list[-i + 1] * (list_size + 2 - i)) % MOD

def comb(n, r):
	if n < r:
		return 0
	elif n == 0 or r == 0 or n == r:
		return 1
	else:
		return (((f_list[n - 1] * f_r_list[n - r - 1]) % MOD) * f_r_list[r - 1]) % MOD 

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

Solution (Python):
class comb():
  F = [1, 1]
  Fi = [1, 1]
  I = [0, 1]
  def __init__(self, num, mod):
    self.MOD = mod
    for i in range(2, num + 1):
      self.F.append((self.F[-1] * i) % mod)
      self.I.append(mod - self.I[mod % i] * (mod // i) % mod)
      self.Fi.append(self.Fi[-1] * self.I[i] % mod)
  def com(self, n, k):
    if n < k: return 0
    if n < 0 or k < 0: return 0
    return self.F[n] * (self.Fi[k] * self.Fi[n - k] % self.MOD) % self.MOD


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

com = comb(M * 3 + N, MOD)

t = com.com(M * 3 + N - 1, N - 1)
t -= com.com(M - 1 + N - 1, N -1) * N

while t < 0:
  t += MOD
if N > M:
  for k in range(M + 2, N + 1, 2):
    i = (M * 3 - k) // 2
    t -= com.com(N, k) * com.com(N - 1 + i, i) % MOD
  while t < 0:
    t += MOD

print(t)

Solution (Python):
from sys import exit, setrecursionlimit, stdin, stderr
from functools import reduce
from itertools import *
from collections import *
from bisect import *
from heapq import *
import operator
import math
from fractions import gcd
from random import randint

setrecursionlimit(10**8)

def input():
  return stdin.readline().strip()

def read():
  return int(input())
 
def reads():
  return [int(x) for x in input().split()]

MAX = 3 * 10**6
MOD = 998244353

fact = [1] * MAX
for i in range(1, MAX):
  fact[i] = (fact[i-1] * i) % MOD
invfact = [0] * MAX
invfact[-1] = pow(fact[-1], MOD-2, MOD)
for i in range(MAX-2, -1, -1):
  invfact[i] = invfact[i+1] * (i+1) % MOD

def comb(n, k):
  return fact[n] * invfact[n-k] * invfact[k] % MOD if 0 <= k <= n else 0

def H(n, k):
  return comb(n+k-1, k-1)

N, M = reads()

ans =0
for k in range(M%2, M+1, 2):
  d = comb(N, k) * (
    H((3 * M - k) // 2, N)
          - k * H((M-k)//2, N)
          - (N-k) * H((M-2-k)//2, N)
    )
  
  ans = (ans + d) % MOD

print(ans)

Solution (Python):
mod = 998244353
frac = [1]*3000000
for i in range(2,3000000):
    frac[i] = i * frac[i-1]%mod
fraci = [None]*3000000
fraci[-1] = pow(frac[-1], mod -2, mod)
for i in range(-2, -3000000-1, -1):
    fraci[i] = fraci[i+1] * (3000000+ i + 1) % mod
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return frac[a]*fraci[b]*fraci[a-b]%mod


N, M = map(int, input().split())
A = (comb(3*M+N-1, 3*M) - N*(comb(N+M-2, M-1)))%mod
B = 0
if N >= M+1:
    for j in range(M+1, 3*M+1):
        if (3*M-j)%2 == 1:
            continue
        k = (3*M-j)//2
        B = (B+comb(N, j)*comb(k+N-1,k))%mod
        
print((A-B)%mod)

Solution (Python):
P = 998244353
fa = [1]
for i in range(1, 3000100):
    fa.append(fa[-1]*i%P)
fainv = [pow(fa[-1], P-2, P)%P]
for i in range(1, 3000100)[::-1]:
    fainv.append(fainv[-1]*i%P)
fainv = fainv[::-1]

def C(a, b):
    return fa[a]*fainv[a-b]*fainv[b]%P

N, M = map(int, input().split())
if N == 2:
    print(M+1)
else:
    ans = C(3*M+N-1, N-1)
    for K in range(2*M+1, 3*M+1):
        ans = (ans - N*C(3*M-K+N-2, N-2)) % P
    for K in range(M):
        if 3*M-K*2 <= N:
            ans = (ans - C(K+N-1, N-1)*C(N, 3*M-K*2)) % P
    print(ans)