#include <bits/stdc++.h>

using namespace std;

const int N = 1005, P = 1e9 + 7;

int ksm(int x, int y = P - 2) {
    int res = 1;
    while (y) {
        if (y & 1) res = 1ll * res * x % P;
        x = 1ll * x * x % P;
        y >>= 1;
    }
    return res;
}

int fac[N], iFac[N], f[N][N], tf[N][N], g[N];

int C(int x, int y) {
    return 1ll * fac[x] * iFac[y] % P * iFac[x - y] % P;
}

int Catalan(int x) {
    --x;
    return 1ll * C(x + x, x) * ksm(x + 1) % P;
}

typedef vector<int> poly;

poly operator + (poly x, poly y) {
    if (x.size() < y.size()) swap(x, y);
    poly res = x;
    for (int i = 0; i < (int)y.size(); ++i) {
        res[i] = (res[i] + y[i]) % P;
    }
    return res;
}

poly operator * (poly x, poly y) {
    poly res(x.size() + y.size() - 1);
    for (int i = 0; i < (int)x.size(); ++i) {
        for (int j = 0; j < (int)y.size(); ++j) {
            res[i + j] = (res[i + j] + 1ll * x[i] * y[j]) % P;
        }
    }
    return res;
}

signed main() {
    int x, y, z, v; scanf("%d%d%d%d", &x, &y, &z, &v);
    int n = x + y + z, m = n / 2, L = max({x, y, z});
    if ((n & 1) || L > m) return 0 & puts("0");
    for (int i = fac[0] = iFac[0] = 1; i <= n; ++i) {
        fac[i] = 1ll * fac[i - 1] * i % P;
        iFac[i] = ksm(fac[i]);
    }
    f[0][0] = 1;
    for (int i = 1; i <= L; ++i) {
        for (int j = 1; j <= i; ++j) {
            for (int k = 1; k <= i; ++k) {
                f[i][j] = (f[i][j] + 1ll * f[i - k][j - 1] * Catalan(k)) % P;
            }
        }
    }
    for (int i = 0; i <= L; ++i) {
        for (int j = 0; j <= m; ++j) {
            for (int k = 0, now = 1; k <= i; ++k) {
                tf[i][j] = (tf[i][j] + 1ll * f[i][k] * now % P * iFac[k]) % P;
                now = 1ll * now * j % P;
            }
        }
    }
    for (int a = 0; a <= x; ++a) {
        for (int b = 0; b <= y; ++b) {
            int c = m - a - b;
            if (c > z || c < 0) continue;
            int coef = 0;
            for (int cab = 0; cab <= x - a && cab <= b; ++cab) {
                int cac = x - a - cab, ccb = b - cab, cca = z - c - ccb, cba = a - cca, cbc = c - cac;
                if (cab >= 0 && cac >= 0 && ccb >= 0 && cca >= 0 && cba >= 0 && cbc >= 0) {
                    coef = (coef + 1ll * C(cba + cca, cba) * C(cab + ccb, cab) % P * C(cac + cbc, cac)) % P;
                }
            }
            for (int i = 0; i <= m; ++i) {
                g[i] = (g[i] + 1ll * tf[a][i] * tf[b][i] % P * tf[c][i] % P * coef) % P;
            }
        }
    }
    poly sum;
    for (int i = 0; i <= m; ++i) {
        poly res = {g[i]};
        for (int j = 0; j <= m; ++j) {
            if (i != j) res = res * (poly){(int)(-1ll * j * ksm(i - j) % P), ksm(i - j)};
        }
        sum = sum + res;
    }
    int res = 0;
    for (int i = 1; i <= m; ++i) res = (res + 1ll * fac[i] * sum[i] % P * ksm(i, v)) % P;
    printf("%d\n", (res % P + P) % P);
    return 0;
}