#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <utility>
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <cstdio>
#include <limits>
#define rep(i,n) for(int i = 0; i < n; ++i)
#define rep1(i,n) for(int i = 1; i <= n; ++i)
using namespace std;
template<class T>bool chmax(T &a, const T &b) { if(a < b){ a = b; return 1; } return 0; }
template<class T>bool chmin(T &a, const T &b) { if(a > b){ a = b; return 1; } return 0; }
template<class T> inline int  sz(T &a) { return a.size(); }
using ll = long long; using ld = long double;
using pi = pair<int,int>; using pl = pair<ll,ll>;
using vi = vector<int>; using vvi = vector<vi>;
using vl = vector<ll>; using vvl = vector<vl>;
const int inf = numeric_limits<int>::max();
const ll infll = numeric_limits<ll>::max();
// ユークリッドの互除法で最大公約数を求める
ll gcd(ll a,ll b){
  if(b==0) return a;
  return gcd(b,a%b);
}
// 最小公倍数
ll lcm(ll a,ll b){
  return a/gcd(a,b)*b;
}
// 拡張ユークリッドの互除法
// ap + bq = gcd(a, b) となる (p, q) を求める.
// 返り値はgcd(a, b)
ll ExtGcd(ll a, ll b, ll &p, ll &q) {  
  if (b == 0) {
    p = 1; q = 0;
    return a;
  }  
  ll d = ExtGcd(b, a%b, q, p);  
  q -= a/b * p;  
  return d;  
}

// 中国剰余定理
// x = b1 (mod m1), x = b2 (mod m2), m1とm2は互いに素, を満たすようなxを求める. 
// 返り値を (r, m) とすると解は x = r (mod m)
// 解なしの場合は (0, -1) を返すよ
pair<ll, ll> CRT(ll b1, ll m1, ll b2, ll m2) {
  ll p, q;
  ll d = ExtGcd(m1, m2, p, q); // p is inv of m1/d (mod. m2/d)
  if ((b2 - b1) % d != 0) return {0, -1};
  ll m = m1 * (m2/d); // lcm of (m1, m2)
  ll tmp = (b2 - b1) / d * p % (m2/d);
  ll r = ((b1 + m1 * tmp) % m + m) % m;
  return {r, m};
}


int main()
{
  ll MOD = 1e+9+7;
  int t; cin >> t;
  while(t-- > 0) {
    vl a(3); rep(i,3) cin >> a[i];
    ll n; cin >> n;
    sort(a.begin(),a.end());
    ll x = a[0], y = a[1], z = a[2];
    ll res = 0;
    while(n > 0) {
      pl p = CRT(n%y, y, 0, x);
      if(p.second != -1) {
        ll r = p.first;
        ll g = lcm(x,y);      
        if(n-r >= 0) {
          ll a = r/x, b = (n-r)/y;
          ll aa = g/x, bb = g/y;
          res = (res + 1 + a/aa + b/bb) % MOD;
        }
      }
      n -= z;
    }
    cout << res << "\n";
  }
  return 0;
}