#include <bits/stdc++.h>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
#define chmin(x,y) x = min((x),(y));
#define chmax(x,y) x = max((x),(y));
using namespace std;
using ll = long long ;
using P = pair<int,int> ;
using pll = pair<long long,long long>;
const int INF = 1e9;
const long long LINF = 1e17;
const int MOD = 1000000007;
//const int MOD = 998244353;
const double PI = 3.14159265358979323846;

template<int mod> struct ModInt{
    long long x=0; 
    constexpr ModInt(long long x=0):x((x%mod+mod)%mod){}
    constexpr ModInt operator+(const ModInt& r)const{return ModInt(*this)+=r;}
    constexpr ModInt operator-(const ModInt& r)const{return ModInt(*this)-=r;}
    constexpr ModInt operator*(const ModInt& r)const{return ModInt(*this)*=r;}
    constexpr ModInt operator/(const ModInt& r)const{return ModInt(*this)/=r;}
    constexpr ModInt& operator+=(const ModInt& r){ if((x+=r.x)>=mod) x-=mod; return *this;}
    constexpr ModInt& operator-=(const ModInt& r){ if((x-=r.x)<0) x+=mod; return *this;}
    constexpr ModInt& operator*=(const ModInt& r){ if((x*=r.x)>=mod) x%=mod; return *this;}
    constexpr ModInt& operator/=(const ModInt& r){ return *this*=r.inv();}
    ModInt inv() const {
        long long s=x,sx=1,sy=0,t=mod,tx=0,ty=1;
        while(s%t!=0){
            long long temp=s/t,u=s-t*temp,ux=sx-temp*tx,uy=sy-temp*ty;
            s=t;sx=tx;sy=ty;
            t=u;tx=ux;ty=uy;
        }
        return ModInt(tx);
    }
    ModInt pow(long long n) const {
        ModInt a=1;
        ModInt b=*this;
        while(n>0){
            if(n&1) a*=b;
            b*=b;
            n>>=1;
        }
        return a;
    }
    friend constexpr ostream& operator<<(ostream& os,const ModInt<mod>& a) {return os << a.x;}
    friend constexpr istream& operator>>(istream& is,ModInt<mod>& a) {return is >> a.x;}
};
using mint = ModInt<MOD>;

ll gcd(ll i,ll j){
    if(j == 0) return i;
    return gcd(j,i%j);
}

//a*x + b*y =gcd(a,b)を満たすx、yをもとめる
pair<long long,long long> extgcd(long long a,long long b){
    long long s=a,sx=1,sy=0,t=b,tx=0,ty=1;
    while(s%t!=0){
        long long temp=s/t,u=s-t*temp,ux=sx-temp*tx,uy=sy-temp*ty;
        s=t;sx=tx;sy=ty;
        t=u;tx=ux;ty=uy;
    }
    return {tx,ty};
}

void solve(){
    ll a,b,c,y;
    cin >> a >> b >> c >> y;
    if(a < b) swap(a,b);
    if(a < c) swap(a,c);
    ll g = gcd(b,c);
    ll b_ = b / g;
    ll c_ = c / g;
    pll p = extgcd(b_,c_);
    mint ans = 0;
    //cout << a << endl;
    //cout << b_ << " " << c_ << endl;
    //cout << p.first << " " << p.second << endl;
    for(ll i=0;i*a<=y;i++){
        //cout << "i"  << i << endl;
        ll x = y - i*a;
        if(x % g != 0) continue;
        ll x_ = x / g;
        pll now = p;
        now.first *= x_;
        now.second *= x_;
        //cout << now.first << " " << now.second << endl;
        if(now.first < 0){
            ll k = (-now.first +c_ - 1)/c_;
            now.first += c_*k;
            now.second -= b_*k;
        }
        //cout << now.first << " " << now.second << endl;
        if(now.second < 0){
            ll k = (-now.second +b_ - 1)/b_;
            now.second += b_*k;
            now.first -= c_*k;
        }
        //cout << now.first << " " << now.second << endl;
        if(now.first < 0 || now.second < 0) continue;
        ans += 1;
        ans += now.first/c_ + now.second/b_;
    }
    cout << ans << endl;
    return;
}

int main(){
    int t;
    cin >> t;
    while(t--){
        solve();
    }
    return 0;
}