#include #include #include #include #include #include #include #include #include #include #include #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; templatebool chmax(T &a, const T &b) { if(a < b){ a = b; return 1; } return 0; } templatebool chmin(T &a, const T &b) { if(a > b){ a = b; return 1; } return 0; } template inline int sz(T &a) { return a.size(); } using ll = long long; using ld = long double; using pi = pair; using pl = pair; using vi = vector; using vvi = vector; using vl = vector; using vvl = vector; const int inf = numeric_limits::max(); const ll infll = numeric_limits::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 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}; } // Modint // modint で宣言 template struct modint{ long long x; long long mod = MOD; modint(long long x=0):x(x%MOD){} modint& operator+=(const modint a){ if((x+=a.x)>=MOD) x-=MOD; return *this; } modint& operator-=(const modint a){ if((x += MOD-a.x)>=MOD) x-=MOD; return *this; } modint& operator*=(const modint a){ (x*=a.x)%=MOD; return *this; } modint operator+(const modint a) const{ modint res(*this); return res+=a; } modint operator-(const modint a) const{ modint res(*this); return res-=a; } modint operator*(const modint a) const{ modint res(*this); return res*=a; } modint pow(long long t) const{ if(!t) return 1; modint a = pow(t>>1); a*=a; if(t&1) a*=*this; return a; } //for prime mod modint inv() const{ return pow(MOD-2); } modint& operator/=(const modint a){ return (*this) *= a.inv(); } modint operator/(const modint a) const{ modint res(*this); return res/=a; } }; using mint = modint<1000000007>; int main() { 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 g = lcm(x,y); mint res; vl dp(y+1, -1); while(n > 0) { if(dp[n%y] >= 0) { res += mint(dp[n%y]); n -= z; continue; } pl p = CRT(n%y, y, 0, x); if (p.second != -1) { ll r = p.first; if(n-r >= 0) { ll a = (r/x)/(g/x), b = ((n-r)/y)/(g/y); res += mint(1 + a + b); dp[n%y] = 1+a+b; } else dp[n%y] = 0; } n -= z; } cout << res.x << "\n"; } return 0; }