#include <bits/stdc++.h>

using i64 = long long;
using u64 = unsigned long long;
using u32 = unsigned;
// TODO: Dynamic ModInt

template<typename T>
constexpr T power(T a, u64 b) {
    T res {1};
    for (; b != 0; b /= 2, a *= a) {
        if (b % 2 == 1) {
            res *= a;
        }
    }
    return res;
}

template<u32 P>
constexpr u32 mulMod(u32 a, u32 b) {
    return 1ULL * a * b % P;
}

template<u64 P>
constexpr u64 mulMod(u64 a, u64 b) {
    u64 res = a * b - u64(1.L * a * b / P - 0.5L) * P;
    res %= P;
    return res;
}

template<typename U, U P>
requires std::unsigned_integral<U>
struct ModIntBase {
public:
    constexpr ModIntBase() : x {0} {}
    
    template<typename T>
    requires std::integral<T>
    constexpr ModIntBase(T x_) : x {norm(x_ % T {P})} {}
    
    constexpr static U norm(U x) {
        if ((x >> (8 * sizeof(U) - 1) & 1) == 1) {
            x += P;
        }
        if (x >= P) {
            x -= P;
        }
        return x;
    }
    
    constexpr U val() const {
        return x;
    }
    
    constexpr ModIntBase operator-() const {
        ModIntBase res;
        res.x = norm(P - x);
        return res;
    }
    
    constexpr ModIntBase inv() const {
        return power(*this, P - 2);
    }
    
    constexpr ModIntBase &operator*=(const ModIntBase &rhs) & {
        x = mulMod<P>(x, rhs.val());
        return *this;
    }
    
    constexpr ModIntBase &operator+=(const ModIntBase &rhs) & {
        x = norm(x + rhs.x);
        return *this;
    }
    
    constexpr ModIntBase &operator-=(const ModIntBase &rhs) & {
        x = norm(x - rhs.x);
        return *this;
    }
    
    constexpr ModIntBase &operator/=(const ModIntBase &rhs) & {
        return *this *= rhs.inv();
    }
    
    friend constexpr ModIntBase operator*(ModIntBase lhs, const ModIntBase &rhs) {
        lhs *= rhs;
        return lhs;
    }
    
    friend constexpr ModIntBase operator+(ModIntBase lhs, const ModIntBase &rhs) {
        lhs += rhs;
        return lhs;
    }
    
    friend constexpr ModIntBase operator-(ModIntBase lhs, const ModIntBase &rhs) {
        lhs -= rhs;
        return lhs;
    }
    
    friend constexpr ModIntBase operator/(ModIntBase lhs, const ModIntBase &rhs) {
        lhs /= rhs;
        return lhs;
    }
    
    friend constexpr std::ostream &operator<<(std::ostream &os, const ModIntBase &a) {
        return os << a.val();
    }
    
    friend constexpr bool operator==(ModIntBase lhs, ModIntBase rhs) {
        return lhs.val() == rhs.val();
    }
    
    friend constexpr bool operator!=(ModIntBase lhs, ModIntBase rhs) {
        return lhs.val() != rhs.val();
    }
    
    friend constexpr bool operator<(ModIntBase lhs, ModIntBase rhs) {
        return lhs.val() < rhs.val();
    }
    
private:
    U x;
};

template<u32 P>
using ModInt = ModIntBase<u32, P>;

template<u64 P>
using ModInt64 = ModIntBase<u64, P>;

constexpr u32 P = 1000000007;
using Z = ModInt<P>;

struct Comb {
    int n;
    std::vector<Z> _fac;
    std::vector<Z> _invfac;
    std::vector<Z> _inv;
    
    Comb() : n{0}, _fac{1}, _invfac{1}, _inv{0} {}
    Comb(int n) : Comb() {
        init(n);
    }
    
    void init(int m) {
        if (m <= n) return;
        _fac.resize(m + 1);
        _invfac.resize(m + 1);
        _inv.resize(m + 1);
        
        for (int i = n + 1; i <= m; i++) {
            _fac[i] = _fac[i - 1] * i;
        }
        _invfac[m] = _fac[m].inv();
        for (int i = m; i > n; i--) {
            _invfac[i - 1] = _invfac[i] * i;
            _inv[i] = _invfac[i] * _fac[i - 1];
        }
        n = m;
    }
    
    Z fac(int m) {
        if (m > n) init(2 * m);
        return _fac[m];
    }
    Z invfac(int m) {
        if (m > n) init(2 * m);
        return _invfac[m];
    }
    Z inv(int m) {
        if (m > n) init(2 * m);
        return _inv[m];
    }
    Z binom(int n, int m) {
        if (n < m || m < 0) return 0;
        return fac(n) * invfac(m) * invfac(n - m);
    }
} comb;

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    
    int A, B, C, x;
    std::cin >> A >> B >> C >> x;
    
    int S = A + B + C;
    
    if (S % 2 != 0 || std::max({A, B, C}) > S / 2) {
        std::cout << 0 << "\n";
        return 0;
    }
    
    int M = std::max({A, B, C});
    
    std::vector<Z> iw(M + 1);
    for (int i = 0; i <= M; i++) {
        iw[i] = comb.binom(2 * i + 1, i) * comb.inv(2 * i + 1);
    }
    std::vector<Z> w(M + 1);
    w[0] = 1;
    for (int i = 1; i <= M; i++) {
        for (int j = 0; j < i; j++) {
            w[i] -= w[j] * iw[i - j];
        }
    }
    w[0] = 0;
    for (int i = 1; i <= M; i++) {
        w[i] *= -1;
    }
    
    std::vector<Z> y(S / 2 + 1);
    
    std::vector dp(M + 1, std::vector<Z>(M + 1));
    dp[0][0] = 1;
    for (int x = M; x >= 1; x--) {
        for (int i = M; i >= 0; i--) {
            for (int j = i / x; j >= 0; j--) {
                Z pw = 1;
                for (int a = 1; i + a * x <= M; a++) {
                    pw *= w[x];
                    dp[i + a * x][j + a] += dp[i][j] * pw * comb.invfac(a);
                }
            }
        }
    }
    
    std::vector dpy(S / 2 + 1, std::vector<Z>(M + 1));
    for (int x = 0; x <= S / 2; x++) {
        for (int i = 0; i <= M; i++) {
            for (int j = i; j >= 0; j--) {
                dpy[x][i] = dpy[x][i] * x + dp[i][j];
            }
        }
    }
    for (int aa = 0; aa <= A; aa++) {
        for (int bb = 0; bb <= B; bb++) {
            int cc = S / 2 - aa - bb;
            if (cc < 0 || cc > C) {
                continue;
            }
            
            std::vector<Z> f(S / 2 + 1);
            for (int x = 0; x <= S / 2; x++) {
                f[x] = dpy[x][aa] * dpy[x][bb] * dpy[x][cc];
            }
            
            Z sum = 0;
            for (int ab = 0; aa + ab <= A && ab <= bb; ab++) {
                int ac = A - aa - ab;
                int cb = bb - ab;
                int ca = C - cc - cb;
                int ba = aa - ca;
                int bc = B - bb - ba;
                
                if (ac < 0 || cb < 0 || ca < 0 || ba < 0 || bc < 0) {
                    continue;
                }
                
                sum += comb.binom(aa, ba) * comb.binom(bb, cb) * comb.binom(cc, ac);
            }
            
            for (int x = 0; x <= S / 2; x++) {
                y[x] += f[x] * sum;
            }
        }
    }
    
    std::vector<Z> prod(S / 2 + 1), g(S / 2 + 1);
    g[0] = y[0];
    prod[0] = 1;
    for (int i = 1; i <= S / 2; i++) {
        Z cur = 0;
        for (int j = i - 1; j >= 0; j--) {
            prod[j + 1] += prod[j];
            prod[j] *= -(i - 1);
            cur = cur * i + g[j];
        }
        Z w = (y[i] - cur) * comb.invfac(i);
        for (int j = 0; j <= i; j++) {
            g[j] += prod[j] * w;
        }
    }
    
    Z ans = 0;
    for (int i = 0; i <= S / 2; i++) {
        g[i] *= comb.fac(i);
        ans += power(Z(i), x) * g[i];
    }
    std::cout << ans << "\n";
    
    return 0;
}