Example

Input

3 1 0 1


Output

18
#include <bits/stdc++.h>
using namespace std;

using ll=long long;
#define int ll

#define FOR(i,a,b) for(int i=int(a);i<int(b);i++)
#define REP(i,b) FOR(i,0,b)
#define MP make_pair
#define PB push_back
#define ALL(x) x.begin(),x.end()
#define REACH cerr<<"reached line "<<__LINE__<<endl
#define DBG(x) cerr<<"line "<<__LINE__<<" "<<#x<<":"<<x<<endl

using pii=pair<int,int>;
using vi=vector<int>;
using ld=long double;

template<class T,class U>
ostream& operator<<(ostream& os,const pair<T,U>& p){
	os<<"("<<p.first<<","<<p.second<<")";
	return os;
}

template<class T>
ostream& operator <<(ostream& os,const vector<T>& v){
	os<<"[";
	REP(i,(int)v.size()){
		if(i)os<<",";
		os<<v[i];
	}
	os<<"]";
	return os;
}

int read(){
	int i;
	scanf("%lld",&i);
	return i;
}

void print(int x){
	printf("%lld",x);
}

void printSpace(){
	printf(" ");
}

void printEoln(){
	printf("\n");
}

string readString(){
	static char buf[3341919];
	scanf("%s",buf);
	return string(buf);
}

char* readCharArray(){
	static char buf[3341919];
	static int bufUsed=0;
	char* ret=buf+bufUsed;
	scanf("%s",ret);
	bufUsed+=strlen(ret)+1;
	return ret;
}

template<class T,class U>
void chmax(T& a,U b){
	if(a<b)
		a=b;
}

template<class T,class U>
void chmin(T& a,U b){
	if(a>b)
		a=b;
}

template<class T>
T Sq(const T& t){
	return t*t;
}

const int mod=1000000007;
template<class T,class U>
void add(T& a,U b){
	a=((int)a+b)%mod;
}

template<class T,class U>
void sub(T& a,U b){
	a=((ll)a-b%mod+mod)%mod;
}

template<class T,class U>
void mult(T& a,U b){
	a=((ll)a*b)%mod;
}

int modPow(int a,int p){
	int s=1;
	REP(i,30){
		if((p>>i)&1)
			mult(s,a);
		mult(a,a);
	}
	return s;
}

int modInv(int a){
	return modPow(a,mod-2);
}

int gcd(int a,int b){
	return b?gcd(b,(a%b)):a;
}

signed main(){
	int n=read();
	int a,b;
	{
		int x=read(),y=read(),z=read();
		a=x+z;
		b=y;
		if(a<b)swap(a,b);
	}
	if(a==0&&b==0){
		vi dp(n+1,0);
		dp[0]=1;
		FOR(i,1,n+1){
			add(dp[i],dp[i-1]);
			if(i>=2)
				add(dp[i],dp[i-2]*(i-1));
		}
		REP(i,n)
			mult(dp[n],n);
		cout<<dp[n]<<endl;
		return 0;
	}
	vector<vi> waf(n+1);
	FOR(i,1,n+1)
		waf[i/gcd(i,a-b)].PB(i);
	vi ans(n+1,0);
	ans[0]=1;
	FOR(t,1,n+1)if(!waf[t].empty()){
		vi c[2]={vi(n+1,0),vi(n+1,0)};
		c[0][0]=1;
		c[1][0]=1;
		for(auto v:waf[t]){
			int can=0;
			REP(leap,v)
				if((a+b+2*leap)%v==0&&(a+leap)%gcd(v,a-b)==0&&(b+leap)%gcd(v,a-b)==0)
					can++;
			assert(can<2);
			vi nx(n+1,0);
			int fact=1;
			for(int i=0;i*v<=n;i++){
				int x=modPow(modInv(v),i)*modInv(fact)%mod;
				for(int j=0;j<=n;j+=t)if(j+i*v<=n)
					add(nx[j+i*v],c[can][j]*x);
				mult(fact,i+1);
			}
			swap(nx,c[can]);
		}
		vi kyuri(n+1,0),nx(n+1,0);
		vector<vi> dp(n/t+1,vi(n/t+1,0));
		dp[0][0]=1;
		REP(p,n/t+1)REP(q,n/t+1)if(p+q<=n/t){
			if(p>0){
				add(dp[p][q],dp[p-1][q]);
				if(p>=2)
					add(dp[p][q],dp[p-2][q]*(p-1)*t);
				if(q>0)
					add(dp[p][q],dp[p-1][q-1]*q*t);
			}else{
				if(q>=2)
					add(dp[p][q],dp[p][q-2]*((q-1)*t%mod));
			}
			add(kyuri[(p+q)*t],c[1][p*t]*c[0][q*t]%mod*dp[p][q]);
		}
		for(int i=0;i<=n;i+=t)
			REP(j,n+1)if(j+i<=n)
				add(nx[j+i],ans[j]*kyuri[i]);
		swap(nx,ans);
	}
	FOR(i,1,n+1)
		mult(ans[n],i);
	cout<<ans[n]<<endl;
}