#include <bits/stdc++.h>

#define _overload(_1,_2,_3,name,...) name
#define _rep(i,n) _range(i,0,n)
#define _range(i,a,b) for(int i=int(a);i<int(b);++i)
#define rep(...) _overload(__VA_ARGS__,_range,_rep,)(__VA_ARGS__)

#define _rrep(i,n) _rrange(i,n,0)
#define _rrange(i,a,b) for(int i=int(a)-1;i>=int(b);--i)
#define rrep(...) _overload(__VA_ARGS__,_rrange,_rrep,)(__VA_ARGS__)

#define _all(arg) begin(arg),end(arg)
#define uniq(arg) sort(_all(arg)),(arg).erase(unique(_all(arg)),end(arg))
#define getidx(ary,key) lower_bound(_all(ary),key)-begin(ary)
#define clr(a,b) memset((a),(b),sizeof(a))
#define bit(n) (1LL<<(n))
#define popcount(n) (__builtin_popcountll(n))

using namespace std;

template<class T>bool chmax(T &a, const T &b) { return (a < b) ? (a = b, 1) : 0;}
template<class T>bool chmin(T &a, const T &b) { return (b < a) ? (a = b, 1) : 0;}

using ll = long long;
using R = long double;
const R EPS = 1e-9L; // [-1000,1000]->EPS=1e-8 [-10000,10000]->EPS=1e-7
inline int sgn(const R& r) {return (r > EPS) - (r < -EPS);}
inline R sq(R x) {return sqrt(max(x, 0.0L));}

const int dx[8] = {1, 0, -1, 0, 1, -1, -1, 1};
const int dy[8] = {0, 1, 0, -1, 1, 1, -1, -1};

// Problem Specific Parameter:


template<int MOD> struct ModInt {
	static const int Mod = MOD;
	unsigned x;
	ModInt(): x(0) {}
	ModInt(signed sig) {int sigt = sig % MOD; if (sigt < 0) sigt += MOD; x = sigt; }
	ModInt(signed long long sig) { int sigt = sig % MOD; if (sigt < 0) sigt += MOD; x = sigt; }
	int getval() const { return int(x); }
	inline int extgcd(int a, int b, int& x, int& y) {x = 1, y = 0; int g = a; if (b != 0) g = extgcd(b, a % b, y, x), y -= a / b * x; return g;}


	ModInt &operator+=(ModInt that) { if ((x += that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator-=(ModInt that) { if ((x += MOD - that.x) >= MOD) x -= MOD; return *this; }
	ModInt &operator*=(ModInt that) { x = 1LL * x * that.x % MOD; return *this; }
	ModInt &operator/=(ModInt that) { int y, d; extgcd(that.x, MOD, y, d); return x = 1LL * x * y % MOD; return *this;}

	ModInt operator+(ModInt that) const { return ModInt(*this) += that; }
	ModInt operator-(ModInt that) const { return ModInt(*this) -= that; }
	ModInt operator*(ModInt that) const { return ModInt(*this) *= that; }
	ModInt operator/(ModInt that) const { return ModInt(*this) /= that; }
};

using mint = ModInt<1000000007>;

const int limit = 100000;
mint I[limit + 10];
void inv_table(int n) {I[1] = 1LL; rep(i, 2, n + 1) I[i] = I[mint::Mod % i] * (mint::Mod - mint::Mod / i);}

mint F[limit + 10], IF[limit + 10];
void fact_table(int n) {F[0] = IF[0] = 1LL; rep(i, 1, n + 1) F[i] = F[i - 1] * i, IF[i] = IF[i - 1] * I[i];}

inline mint comb(int n, int k) {return F[n] * IF[k] * IF[n - k];}


void init() {
	inv_table(limit);
	fact_table(limit);
}

mint dp[5010][5010][2];

int main(void) {
	init();

	int n, a, b, c;
	cin >> n >> a >> b >> c;

	dp[0][0][0] = mint(1);

	rep(block, n + 1)rep(three, n + 1)rep(state, 2) {
		if (dp[block][three][state].getval() == 0) continue;

		if (state == 0) {
			// 2222
			// 232 ... 3333 ... 232 first
			dp[block][three][1] += dp[block][three][0];
		} else {
			// 232 ... 3333 ... 232 middle
			dp[block][three + 1][1] += dp[block][three][1];
			// 232 ... 3333 ... 232 last
			dp[block + 1][three][0] += dp[block][three][1];
		}
	}

	// ca 2222  cb 232232 cc 3333 cd 333333 ce 3113 cf 11

	// a = ce + cf
	// b = 2 * ca + 2 * cb
	// c = cb + 2 * cc + 3 * cd + ce

	// block = ca + cb
	// three = cb + 2 * cc

	mint ans;
	rep(three, n + 1)rep(ce, n + 1) {
		if (b % 2 != 0) continue;
		const int block = b / 2;
		if (dp[block][three][0].getval() == 0) continue;
		if (a - ce < 0) continue;
		const int cf = a - ce;
		if (c - three - ce < 0 or (c - three - ce) % 3 != 0) continue;
		const int cd = (c - three - ce) / 3;

		mint cur = dp[block][three][0];
		mint coef = comb(block + cd + ce + cf, block) * comb(cd + ce + cf, cd) * comb(ce + cf, ce);
		ans += coef * cur;
	}

	cout << ans.getval() << endl;
	return 0;
}


Measured memory_bytes: 197632