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

struct uns_t {} uns;
template <typename Element, typename Head, typename ...Args>
auto vec(Element init, Head arg, Args ...args) {
    if constexpr (sizeof...(Args) == 0) return std::vector(arg, init);
    else return std::vector(arg, vec(init, args...));
}
template <typename Element, typename Head, typename ...Args>
auto vec(uns_t, Head arg, Args ...args) {
    return vec(Element(), arg, args...);
}

tuple<int64_t, int64_t, int64_t> exgcd(int64_t a, int64_t b) {
    if (b != 0) {
        auto [xp, yp, gcd] = exgcd(b, a % b);
        return { yp, xp - a / b * yp, gcd };
    }
    else {
        return { 1, 0, a };
    }
}

int main() {
    constexpr int64_t mod = 1e9 + 7;

    int t; cin >> t;
    while (t--) {
        int64_t n, k, h, y; cin >> n >> k >> h >> y;
        tie(n, k, h) = make_tuple(
            min({ n, k, h }),
            n + k + h - min({ n, k, h }) - max({ n, k, h }),
            max({ n, k, h })
        );

        int64_t ans = 0;

        auto [n1, n2, gcd] = exgcd(n, k);
        auto np = n / gcd, kp = k / gcd;
        for (int i = 0; i * h <= y; ++i) {
            auto c = y - i * h;

            if (c == 0) {
                ++ans;
                ans %= mod;
            }
            else {
                ans += c % n == 0;
                ans %= mod;
                ans += c % k == 0;
                ans %= mod;
            }

            if (c % gcd) continue;
            auto a = n1 * (c / gcd), b = n2 * (c / gcd);

            auto X = [&](auto n) { return a - kp * n; };
            auto Y = [&](auto n) { return b + np * n; };

            int64_t xmax = 0, ymin = 0;
            if (a == 0) {
                xmax = -1;
            }
            if (0 < a) {
                xmax = a / kp - (a % kp == 0);
            }
            if (a < 0) {
                xmax = a / kp - 1;
            }

            if (b == 0) {
                ymin = 1;
            }
            if (0 < b) {
                ymin = -b / np + (b % np == 0);
            }
            if (b < 0) {
                ymin = -b / np + 1;
            }

            if (ymin <= xmax) {
                ans += xmax - ymin + 1;
                ans %= mod;
            }
        }

        cout << ans << endl;
    }
}