Example

Input

3 1 0 1


Output

18
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> P;
typedef pair<int,P> P1;
typedef pair<P,P> P2;
#define pu push
#define pb push_back
#define mp make_pair
#define eps 1e-7
#define INF 1000000000
#define mod 1000000007
#define fi first
#define sc second
#define rep(i,x) for(int i=0;i<x;i++)
#define repn(i,x) for(int i=1;i<=x;i++)
#define SORT(x) sort(x.begin(),x.end())
#define ERASE(x) x.erase(unique(x.begin(),x.end()),x.end())
#define POSL(x,v) (lower_bound(x.begin(),x.end(),v)-x.begin())
#define POSU(x,v) (upper_bound(x.begin(),x.end(),v)-x.begin())
int n;
ll x,y,z,A,B;
ll modpow(ll x,ll n){
	ll res=1;
	while(n>0){
		if(n&1) res=res*x%mod;
		x=x*x%mod;
		n>>=1;
	}
	return res;
}
ll F[2005],R[2005];
void make(){
	F[0] = 1;
	for(int i=1;i<2005;i++) F[i] = F[i-1]*i%mod;
	for(int i=0;i<2005;i++) R[i] = modpow(F[i],mod-2);
}
ll C(int a,int b){
	return F[a]*R[b]%mod*R[a-b]%mod;
}
ll gcdd(ll a,ll b){
	if(a < b) swap(a,b);
	if(b == 0) return a;
	return gcdd(b,a%b);
}
ll e[1005][1005];
void make2(){
    for(int j=1;j<=1000;j++){
        e[0][j] = 1;
        for(int i=2;i<=1000;i+=2){
            e[i][j] = e[i-2][j];
            e[i][j] = e[i][j]*C(i,2)%mod*1LL*j%mod;
            e[i][j] = e[i][j]*modpow(i/2,mod-2)%mod;
        }
    }
}
int main(){
	make();
	make2();
	cin>>n>>x>>y>>z;
	A = x+z, B = y;
	if(A < B) swap(A,B);
	if(x+y+z == 0){
	    ll cnt = 0;
		ll cur = 1;
		for(int i=0;i<=n/2;i++){
			cnt = (cnt+cur)%mod;
			cur = cur*C(n-i*2,2)%mod;
			cur = cur*modpow(i+1,mod-2)%mod;
		}
		cout << (cnt*modpow(n,n)%mod+mod)%mod << endl;
	}
	else if(A == B){
		ll cnt = 0;
		ll cur = 1;
		for(int i=0;i<=n/2;i++){
			cnt = (cnt+cur)%mod;
			cur = cur*C(n-i*2,2)%mod;
			cur = cur*modpow(i+1,mod-2)%mod;
		}
		cout << (cnt*F[n]%mod+mod)%mod << endl;
	}
	else{
		ll S = A-B;
		//T^S = P^2
		ll dp[1005][2] = {};
		int cur = 0, nxt = 1;
		dp[0][0] = 1;
		for(int i=1;i<=n;i++){
		    for(int j=0;j<1005;j++) dp[j][nxt] = 0;
			vector<int>ok;
			for(int a=1;i*a<=n;a++){
				int cyc = i*a;
				ll gcd = gcdd(1LL*cyc,S);
				if(gcd == a) ok.pb(i*a);
			}
			if(ok.empty()) continue;
			ll t[1005][2] = {};
			int cc = 0, nn = 1;
			t[0][0] = 1;
			for(int j=0;j<ok.size();j++){
				rep(k,1005) t[k][nn] = 0;
				for(int k=0;k<=n;k++){
					if(t[k][cc] == 0) continue;
					ll mul = 1;
					for(int l=0;k+ok[j]*l<=n;l++){
						//add l ko
						t[k+ok[j]*l][nn] += t[k][cc]*mul%mod;
						int nw = k+ok[j]*l;
						for(int w=0;w<ok[j]/i;w++){
							mul = mul*C(nw+i,i)%mod*F[i-1]%mod;
							nw += i;
						}
						mul = mul*R[ok[j]/i]%mod;
						mul = mul*modpow(l+1,mod-2)%mod;
						mul = mul*F[ok[j]/i-1]%mod*modpow(i,ok[j]/i-1)%mod;
					}
				}
				swap(cc,nn);
				for(int k=0;k<=n;k++) t[k][cc] = (t[k][cc]%mod+mod)%mod;
			}
			ll add[1005]={};
			for(int j=0;j<=n;j+=i){
				if(i%2 == 0){
					if((j/i)%2 == 1) add[j] = 0;
					else{
					    add[j] = t[j][cc]*e[j/i][i]%mod;
					}
				}
				else{
				    int cur = j/i;
				    ll sum = 0,cnt = 1,f = 1;
				    for(int j=0;j<=cur/2;j++){
				        sum += cnt*f%mod;
				        sum %= mod;
				        cnt = cnt*C(cur-2*j,2)%mod;
				        cnt = cnt*modpow(j+1,mod-2)%mod;
				        f = f*i%mod;
				    }
					add[j] = t[j][cc]*sum%mod;
				}
			}
			for(int j=0;j<=n;j++){
				if(dp[j][cur] == 0) continue;
				for(int k=0;j+k<=n;k+=i){
					dp[j+k][nxt] += dp[j][cur] * add[k] % mod * C(j+k,k) % mod;
				}
			}
			swap(cur,nxt);
			for(int j=0;j<=n;j++) dp[j][cur] = (dp[j][cur]%mod+mod)%mod;
		}
		cout << dp[n][cur] << endl;
	}
}
