MOD = 10**9 + 7

def main():
    import sys
    N, R, G, B = map(int, sys.stdin.readline().split())
    sum_global = R + G + B
    if sum_global > N:
        print(0)
        return
    
    max_fact = 6000  # Sufficient for comb(rem +k, k) where rem +k can be up to N + sum_global
    # Precompute factorial and inverse factorial
    fact = [1] * (max_fact + 1)
    for i in range(1, max_fact + 1):
        fact[i] = fact[i-1] * i % MOD
    inv_fact = [1] * (max_fact + 1)
    inv_fact[max_fact] = pow(fact[max_fact], MOD-2, MOD)
    for i in range(max_fact -1, -1, -1):
        inv_fact[i] = inv_fact[i+1] * (i+1) % MOD
    
    pow2 = [1] * (3000 * 3 + 1)
    for i in range(1, 3000 *3 +1):
        pow2[i] = pow2[i-1] * 2 % MOD
    
    def comb(n, k):
        if n <0 or k <0 or n <k:
            return 0
        return fact[n] * inv_fact[k] % MOD * inv_fact[n -k] % MOD
    
    answer = 0
    
    max_p = sum_global // 2
    for p in range(0, max_p + 1):
        s = sum_global - 2 * p
        if s < 0:
            continue
        k = s + p
        # compute rem_blanks = N - sum_global - k + 1
        rem = N - sum_global - k + 1
        if rem <0:
            continue
        
        c_comb = comb(rem +k, k)
        if c_comb ==0:
            continue
        
        total =0
        a_min = max(0, p - B)
        a_max = min(p, R, G)
        if a_min > a_max:
            continue
        
        for a in range(a_min, a_max +1):
            remaining = p - a
            # Compute valid b range
            b_min = max(0, p - G)
            b_max = min(R -a, remaining)
            if b_min > b_max:
                continue
            if b_max <0:
                continue
            
            # iterate b and c = remaining -b
            for b in range(b_min, b_max +1):
                c = remaining - b
                # Check constraints on c
                if c <0:
                    continue
                # Check a + c <= G --> a + c = a + (remaining -b) = a + p -a -b = p -b <= G --> p -b <=G --> b >= p - G (already considered in b_min)
                # Check b + c <= B --> c <= B -b
                if b + c > B:
                    continue
                # Compute single runs
                s_r = R - a - b
                s_g = G - a - c
                s_b = B - b - c
                if s_r <0 or s_g <0 or s_b <0:
                    continue
                
                # Calculate denominator: s_r! s_g! s_b! a! b! c!
                denom = fact[s_r] * fact[s_g] % MOD
                denom = denom * fact[s_b] % MOD
                denom = denom * fact[a] % MOD
                denom = denom * fact[b] % MOD
                denom = denom * fact[c] % MOD
                inv_denom = pow(denom, MOD-2, MOD)
                
                term = fact[k] * inv_denom % MOD
                term = term * pow2[a + b + c] % MOD
                total = (total + term) % MOD
        
        answer = (answer + total * c_comb) % MOD
    
    print(answer % MOD)

if __name__ == "__main__":
    main()