Problem:
 Score : 1500 points    Problem Statement  Consider all integers between 1 and 2N , inclusive.
Snuke wants to divide these integers into N pairs such that:   Each integer between 1 and 2N is contained in exactly one of the pairs.  In exactly A pairs, the difference between the two integers is 1 .  In exactly B pairs, the difference between the two integers is 2 .  In exactly C pairs, the difference between the two integers is 3 .   Note that the constraints guarantee that N = A + B + C , thus no pair can have the difference of 4 or more.  Compute the number of ways to do this, modulo 10^9+7 .      Constraints  1 ≤ N ≤ 5000  0 ≤ A, B, C  A + B + C = N         Input The input is given from Standard Input in the following format:  N  A  B  C       Output Print the answer.        Sample Input 1 3 1 2 0      Sample Output 1 2  There are two possibilities: 1-2, 3-5, 4-6 or 1-3, 2-4, 5-6 .       Sample Input 2 600 100 200 300      Sample Output 2 522158867

Solution (Python):
def inverse(k):
    return fast_pow(k, MOD - 2)


def comb(n, k):
    return FACT[n] * INVERSE_FACT[n - k] * INVERSE_FACT[k] % MOD


def fast_pow(x, y):
    if y == 0:
        return 1
    p = fast_pow(x, y // 2) % MOD
    p = p * p % MOD
    if y % 2:
        p = p * x % MOD
    return p


MOD = 1000000007

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

FACT = [0] * (n + 1)
INVERSE_FACT = [0] * (n + 1)

FACT[0] = 1
INVERSE_FACT[0] = 1
for i in range(1, n + 1):
    FACT[i] = FACT[i - 1] * i % MOD
    INVERSE_FACT[i] = inverse(FACT[i])

ans = 0
if b % 2 != 0:
    print(0)
    exit()

for y in range(0, min(a + 1, c + 1)):
    for z in range(0, (c - y) // 3 + 1):
        x = a - y
        rest3 = c - y - 3 * z

        if b == 0 and rest3 != 0:
            continue

        cur_ans = FACT[x + y + z] * INVERSE_FACT[x] % MOD
        cur_ans = cur_ans * INVERSE_FACT[y] % MOD
        cur_ans = cur_ans * INVERSE_FACT[z] % MOD
        cur_ans = cur_ans * comb(x + y + z + 1 + b // 2 - 1, b // 2) % MOD
        cur_ans = cur_ans * comb(b // 2 + rest3 - 1, rest3) % MOD

        ans = (ans + cur_ans) % MOD

print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 7)

"""
・2,3^n,2
・3(1)
・1
・333
"""

import numpy as np

MOD = 10**9 + 7

N,A,B,C = map(int,readline().split())

if B&1:
    print(0)
    exit()

def cumprod(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
    return arr.ravel()[:L]

U = 10**5
x = np.arange(U,dtype=np.int64); x[0] = 1
fact = cumprod(x)
x = np.arange(U,0,-1,dtype=np.int64); x[0] = pow(int(fact[-1]),MOD-2,MOD)
fact_inv = cumprod(x)[::-1]

"""
・2,3^n,2：B2個
・3113：A-n個
・11：n個
・333333：m個
"""

B2 = B//2

answer = 0
for m in range(C//3+1):
    if B2 == 0:
        a = C-3*m; b = m; c = A-a
        if a < 0 or c < 0:
            continue
        x = fact[a+b+c] * fact_inv[a] % MOD
        x *= fact_inv[b]; x %= MOD
        x *= fact_inv[c]; x %= MOD
        answer += x
        continue
    n_min = max(0,A-C+3*m)
    
    x = fact[B2+A+m] * fact_inv[m] % MOD * fact_inv[B2] % MOD
    x *= fact_inv[A-n_min::-1]; x %= MOD
    x *= fact_inv[n_min:A+1]; x %= MOD
    
    
    y = fact[B2-1+C-3*m-A+n_min:B2+C-3*m] * fact_inv[B2-1] % MOD
    y *= fact_inv[C-3*m-A+n_min:C-3*m+1]; y %= MOD
    answer += (x*y%MOD).sum()

answer %= MOD

print(answer)