#include // assert #include // cin, cout, ios #include // swap, pair #include namespace mp = boost::multiprecision; using bigint = boost::multiprecision::cpp_int; // 通常の Truncating 除算における商と剰余の組を返す(商はゼロ方向に切り捨て) template constexpr std::pair divmod(T a, T b) { T q = a / b; T r = a % b; return {q, r}; } // 商を無限大方向に切り下げる除算(剰余が0でない時の符号は除数の符号に一致) template constexpr std::pair divmod_floor(T a, T b) { // 標準の truncating 除算を使う T q = a / b; T r = a % b; // もし符号が食い違っていたら 1 調整 if ((r != 0) && ((r > 0) != (b > 0))) { q -= 1; r += b; } return {q, r}; } // 剰余が非負になる除算(ユークリッド除算) template constexpr std::pair divmod_euclid(T a, T b) { // 標準の truncating 除算を使う T q = a / b; T r = a % b; // 剰余が負なら 1 調整 if (r < 0) { if (b > 0) { q -= 1; r += b; } else { q += 1; r -= b; } } return {q, r}; } template T gcd(T a, T b){while(b != 0){a %= b; std::swap(a, b);}return a;} template bool chmin(T& a, T b){if(a > b){a = b; return true;} return false;} template bool chmax(T& a, T b){if(a < b){a = b; return true;} return false;} /** Max Weighted Floor (mwf) の非再帰実装。 mwf(n,m,a,b,c,d) = max_{0 <= x < n} a*x + b*floor((c*x + d)/m) 前提: - n > 0, m > 0 計算量/メモリ: - 時間: O(log m)(ユークリッド互除法的再帰による構造縮約) - 追加メモリ: O(1) */ template T mwf(T n, T m, T a, T b, T c, T d) { assert(n > 0 && m > 0); // 正規化: c,d を m で割った余りに変換し、a,b を更新、 sum_acc,max_acc を初期化 // c,d が負の場合にも剰余を非負にするため、ユークリッド除算を使う const auto [qc, rc] = divmod_euclid(c, m); c = rc; a = a + b * qc; const auto [qd, rd] = divmod_euclid(d, m); d = rd; T sum_acc = b * qd; T max_acc = sum_acc; while(true) { assert(0 <= c && c < m && 0 <= d && d < m); // 0 ≤ x < n における y = floor((c*x+d)/m) の最大値 y_max を計算 const T y_max = (c * (n - 1) + d) / m; // y_max >= 0 // 現在の小問題における x = 0, n-1 のときの値を max_acc に反映 const T rval = sum_acc + a * (n - 1) + b * y_max; if(max_acc < sum_acc) { max_acc = sum_acc; } if(max_acc < rval) { max_acc = rval; } // x = 0, n-1 のいずれかで最大値を取るのが確定したら終了 if(y_max == 0 || a == 0 || b == 0 || (a > 0 && b > 0) || (a < 0 && b < 0)) { return max_acc; } // 小問題へのパラメータ変換 if(a < 0) { sum_acc += a + b; } n = y_max; d = m - d - 1; std::swap(a, b); std::swap(c, m); assert(0 < n && 0 < m && 0 <= c && 0 <= d); // 正規化: c,d を m で割った余りに変換し、a,b,sum_acc を更新 // c,d,m は非負なので通常の剰余で良い const std::pair divmod_c = divmod(c, m); c = divmod_c.second; a = a + b * divmod_c.first; const std::pair divmod_d = divmod(d, m); d = divmod_d.second; sum_acc = sum_acc + b * divmod_d.first; } } /** * Max Weighted Floor (mwf) の非再帰実装。 * mwf(n,m,a,b,c,d) = max_{0 <= x < n} a*x + b*floor((c*x + d)/m) * * 前提: * * n > 0, m > 0 * * 返り値: mwf(n,m,a,b,c,d) <= z なら True、そうでなければ False を返す。 * * 計算量/メモリ: * * 時間: O(log m)(ユークリッド互除法的再帰による構造縮約) * * 追加メモリ: O(1) * * @param z T * @param n T * @param m T * @param a T * @param b T * @param c T * @param d T * @return bool */ template bool mwf_leq(T z, T n, T m, T a, T b, T c, T d) { assert(0 < n && 0 < m); const auto [qc, rc] = divmod_euclid(c, m); c = rc; a += b * qc; const auto [qd, rd] = divmod_euclid(d, m); d = rd; T sum_acc = -z + b * qd; auto sat = [](T v) { T zero = 0; return v > 0 ? v : zero; }; while (true) { assert(0 < n && 0 < m && 0 <= c && c < m && 0 <= d && d < m); if (sum_acc > 0) { return false; } if (a <= 0 && b <= 0) { return true; } T n1 = n - 1; auto y_max = (c * n1 + d) / m; if (y_max == 0) { return sum_acc + a * n1 <= 0; } if (a >= 0) { if (sum_acc + a * n1 + b * y_max > 0) { return false; } if (b >= 0) { return true; } } else { sum_acc = sum_acc + a + b; } n = y_max; d = m - d - 1; std::swap(a, b); std::swap(c, m); assert(0 < n && 0 < m && 0 <= c && 0 <= d); // 正規化: c,d を m で割った余りに変換し、a,b,sum_acc を更新 // c,d,m は非負なので通常の剰余で良い const auto [qc, rc] = divmod(c, m); c = rc; a += b * qc; const auto [qd, rd] = divmod(d, m); d = rd; sum_acc += b * qd; } } template T mwf_lr(T l, T r, T m, T a, T b, T c, T d) { assert(l < r && m > 0); T n = r - l; const auto [quot, rem] = divmod_euclid(T(c * l + d), m); return a * l + b * quot + mwf(n, m, a, b, c, rem); } template bool mwf_lr_leq(T z, T l, T r, T m, T a, T b, T c, T d) { assert(l < r && m > 0); T n = r - l; const auto [quot, rem] = divmod_euclid(T(c * l + d), m); T z_adjusted = z - a * l - b * quot; return mwf_leq(z_adjusted, n, m, a, b, c, rem); } /** * Δ(D,A,B,x) = floor(x/D) - floor( (floor(x/A) * floor(A*B/D)) / B ) において、 * u_min(D,A,B,K) = min { u >= 0 | Δ(D,A,B,u*D) > K } を半開区間二分探索 [0, A'BK+2) で求め、 * x_min(D,A,B,K) = min { x >= 0 | Δ(D,A,B,x) > K } = u_min(D,A,B,K)*D を返します(解なしは -1)。 * * 前提: * * D > 0, A > 0, B > 0, K >= 0(整数) * 手順概要: * 1) 既約化: g = gcd(D, A), D' = D/g, A' = A/g * 2) (M', R') = divmod(A' * B, D')(A'B = D'*M' + R') * 3) 閾値 T_Δ = B*K を設定 * 4) E(u) = B*u - M'*floor(D'u / A') * 5) F(N) = max_{0 <= u < N} E(u) を mwf で評価(N > 0, m = A' > 0) * 6) 区間 [0, A'BK+2) で述語 [F(u) <= T_Δ] を二分探索し、 * F(u) <= T_Δ となる最大の u 、つまり T_Δ < E(u) となる最小の u を特定。x = D*u を返す。 * * 備考: * * R' = 0 かつ D'K + 1 >= A' のときは解が存在しないため -1 を返します。 * * 解が存在する場合、 u_min は必ず [0, A'BK+2) の範囲に存在します。 * @param D T * @param A T * @param B T * @param K T * @return T */ template T compute_xmin_leq(T D, T A, T B, T K) { assert(D > 0 && A > 0 && B > 0); if (K < 0) { return 0; } T gcd_DA = gcd(D, A); T Dred = D / gcd_DA, Ared = A / gcd_DA; T Mred = (Ared * B) / Dred, Rred = (Ared * B) % Dred; T Tdelta = B * K; T Mred_neg = -Mred; T zero = 0; // 解なしをパラメータを用いて判定 if (Rred == 0 && Dred * K + 1 >= Ared) { return -1; } // [0, hi) の半開区間、緩い上界 A'BK+1 を包括する hi = A'BK+2 を設定 T lo = 0, hi = Ared * B * K + 2; // F(hi) > T の不変条件を確認 assert(!mwf_lr_leq(Tdelta, lo, hi, Ared, B, Mred_neg, Dred, zero)); // F(lo) <= T, F(hi) > T の不変条件で u_min を二分探索 while (lo + 1 < hi) { T mid = (lo + hi) / 2; if (mwf_lr_leq(Tdelta, lo, mid, Ared, B, Mred_neg, Dred, zero)) { lo = mid; } else { hi = mid; } } // lo = u_min, hi = lo + 1 return D * lo; } int main() { std::cin.tie(nullptr); std::ios::sync_with_stdio(false); int t; bigint d, a, b, k, ans; std::cin >> t; for(int i = 0; i < t; i++) { std::cin >> d >> a >> b >> k; assert(1 <= d && 1 <= a && 1 <= b && 0 <= k); ans = compute_xmin_leq(d, a, b, k); std::cout << ans << '\n'; } }