Example

Input

3 1 0 1


Output

18
#include <iostream>
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <string>
#include <vector>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <bitset>
#include <numeric>
#include <utility>
#include <iomanip>
#include <algorithm>
#include <functional>
using namespace std;

const int MOD = 1000000007;
inline long long mod(long long a, long long m) { return (a % m + m) % m; }

struct Fp {
    long long val;
    
    Fp() : val(0) {} 
    Fp(long long val_) { this->val = mod(val_, MOD); }
    Fp operator = (long long val_) { this->val = mod(val_, MOD); return *this; }
    inline Fp operator - () { return mod(-val, MOD); }
    inline const Fp& operator += (const Fp &x);
    inline const Fp& operator -= (const Fp &x);
    inline const Fp& operator *= (const Fp &x);
    inline const Fp& operator /= (const Fp &x);
};

ostream &operator << (ostream &os, Fp x) { return os << x.val;}
istream &operator >> (istream &is, Fp &x) { is >> x; return is; }
bool operator == (Fp x, Fp y) { return mod(x.val, MOD) == mod(y.val, MOD); }

inline Fp operator + (Fp x, Fp y) { return mod(x.val + y.val, MOD); } 
inline Fp operator - (Fp x, Fp y) { return mod(x.val - y.val, MOD); } 
inline Fp operator * (Fp x, Fp y) { return mod(x.val * y.val, MOD); } 
inline Fp operator / (Fp x, Fp y) {
    long long a = y.val, b = MOD, u = 1, v = 0;
    while (b) {
        long long t = a/b;
        a -= t*b; swap(a, b);
        u -= t*v; swap(u, v);
    }
    return x * u;
}
inline Fp pow(Fp a, long long n) {
    if (n == 0) return Fp(1);
    Fp t = pow(a, n/2);
    t = t * t; if (n & 1) t = t * a;
    return t;
}
inline const Fp& Fp::operator += (const Fp &x) {*this = *this + x; return *this;}
inline const Fp& Fp::operator -= (const Fp &x) {*this = *this - x; return *this;}
inline const Fp& Fp::operator *= (const Fp &x) {*this = *this * x; return *this;}
inline const Fp& Fp::operator /= (const Fp &x) {*this = *this / x; return *this;}

int N;
long long X, Y, Z, p;
const int MAX = 1050;

vector<long long> divisor(long long n) {
    vector<long long> res;
    for (long long i = 1LL; i*i <= n && i <= N; ++i) {
        if (n%i == 0LL) {
            res.push_back(i);
            long long temp = n/i;
            if (i != temp && temp <= N) res.push_back(temp);
        }
    }
    sort(res.begin(), res.end());
    return res;
}

long long gcd(long long a, long long b) {
    if (b == 0) return a;
    else return gcd(b, a % b);
}

vector<long long> alt[2];
Fp Com[MAX][MAX];
Fp Fact[MAX];

void pre_calc() {
    memset(Com, 0, sizeof(Com));
    Com[0][0] = 1;
    for (int i = 1; i <= N; ++i) {
        Com[i][0] = 1;
        for (int j = 1; j <= N; ++j) {
            Com[i][j] = Com[i-1][j-1] + Com[i-1][j];
        }
    }
    Fact[0] = Fact[1] = 1;
    for (int i = 2; i <= N; ++i) {
        Fact[i] = Fact[i-1] * i;
    }
}

Fp sdp[2][MAX][MAX];
Fp memo[2][MAX][MAX];
Fp dp[MAX][MAX];

int main() {
    cin >> N >> X >> Y >> Z;
    long long p[2];
    p[0] = 2;
    p[1]= abs(X - Y + Z);
    for (int i = 0; i < 2; ++i) if (p[i] != 0) alt[i] = divisor(p[i]);

    pre_calc();
        
    if (p[1] == 0) {
        Fp gfact = pow(Fp(N), N);

        memset(sdp, 0, sizeof(sdp)); 
        sdp[0][0][0] = 1;
        for (long long r = 0; r <= N; ++r) {
            for (long long ai = 0; ai < alt[0].size(); ++ai) {
                long long a = alt[0][ai];
                    
                for (long long k = 0; k <= N; ++k) {
                    if (r + a*k > N) break;
                        
                    sdp[0][r+a*k][ai+1] += sdp[0][r][ai] * pow(Fact[a-1], k) 
                    * Com[r+a*k][a*k] * Fact[a*k] / (pow(Fact[a], k) * Fact[k]);
                }
            }
        }
        cout << gfact * sdp[0][N][alt[0].size()] << endl;
    }
    else {
        memset(memo, 0, sizeof(memo));
        for (int it = 0; it < 2; ++it) {
            for (long long l = 1; l <= N; ++l) {
                memset(sdp, 0, sizeof(sdp)); 
                sdp[it][0][0] = 1;
                for (long long r = 0; l * r <= N; ++r) {
                    for (long long ai = 0; ai < alt[it].size(); ++ai) {
                        long long a = alt[it][ai];
                        
                        sdp[it][r][ai+1] += sdp[it][r][ai];
                        if (gcd(p[it]/a, l) > 1) continue;
                        for (long long k = 1; k <= N; ++k) {
                            if (l * (r + a*k) > N) break;
                            
                            sdp[it][r+a*k][ai+1] += sdp[it][r][ai] * pow(Fact[a-1] * pow(Fp(l), a-1), k) 
                            * Com[r+a*k][a*k] * Fact[a*k] / (pow(Fact[a], k) * Fact[k]);
                        }
                    }
                }
                for (long long r = 0; l * r <= N; ++r) {
                    memo[it][l][r] = sdp[it][r][alt[it].size()];
                }
            }
        }
        
        memset(dp, 0, sizeof(dp));
        dp[0][0] = 1;
        for (long long l = 1; l <= N; ++l) {
            for (long long num = 0; num <= N; ++num) {
                for (long long r = 0; r <= N; ++r) {
                    if (num+l*r > N) break;
                    dp[l][num+l*r] += dp[l-1][num] * Com[num+l*r][l*r] * Fact[l*r] / Fact[r] / pow(Fp(l), r)
                    * memo[0][l][r] * memo[1][l][r];
                }
            }
        }
        cout << dp[N][N] << endl;
    }  
    return 0;
}
