#1160836 (C++17) No.3398 Accuracy of Integer Division Approximate Function 2

• 問題: No.3398 Accuracy of Integer Division Approximate Function 2
• ユーザー: 👑 Mizar
• 提出日時: 2026-04-20 20:16:40
• 言語: C++17 (gcc 15.2.0 + boost 1.89.0)
• 結果: AC
• 実行時間: 10 ms / 2,000 ms
• コード長: 4,857 bytes
• コンパイル時間: 4,570 ms
• コンパイル使用メモリ: 376,204 KB
• 実行使用メモリ: 6,400 KB
• 最終ジャッジ日時: 2026-04-20 20:16:47
• 合計ジャッジ時間: 5,317 ms
• ジャッジサーバーID （参考情報）: judge2_1 / judge1_0

ファイルパターン	結果
sample	AC * 3
other	AC * 20

ソースコード

#include <cassert> // assert
#include <iostream> // cin, cout, ios
#include <utility> // swap, pair
#include <boost/multiprecision/cpp_int.hpp>
namespace mp = boost::multiprecision;
using bigint = boost::multiprecision::int256_t;

// 剰余が非負になる除算（ユークリッド除算）
template<typename T>
constexpr std::pair<T, T> divmod_euclid(T a, T b) {
    // 標準の truncating 除算を使う
    T q, r;
    mp::divide_qr(a, b, q, r);
    // 剰余が負なら 1 調整
    if (r < 0) {
        if (b > 0) {
            q -= 1;
            r += b;
        } else {
            q += 1;
            r -= b;
        }
    }
    return {q, r};
}

template<typename T> bool chmin(T& a, T b){if(a > b){a = b; return true;} return false;}
template<typename T> bool chmax(T& a, T b){if(a < b){a = b; return true;} return false;}

template<typename T>
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;
    d += c * l;
    T s = a * l;
    T t = s + b * divmod_euclid(d, m).first;
    while (true) {
        const auto [qc, rc] = divmod_euclid(c, m);
        a += b * qc;
        c = rc;
        const auto [qd, rd] = divmod_euclid(d, m);
        s += b * qd;
        d = rd;
        n -= 1;
        T y_max = (c * n + d) / m;
        chmax(t, s);
        chmax(t, s + a * n + b * y_max);
        if (y_max == 0 || (a >= 0 && b >= 0) || (a <= 0 && b <= 0)) {
            return t;
        }
        if (a <= 0) {
            s += a + b;
        }
        n = y_max;
        d = m - d - 1;
        std::swap(a, b);
        std::swap(c, m);
    }
}

template<typename T>
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;
    d += c * l;
    T s = a * l - z;
    if (s + b * divmod_euclid(d, m).first > 0) {
        return false;
    }
    while (true) {
        const auto [qd, rd] = divmod_euclid(d, m);
        s += b * qd;
        d = rd;
        if (s > 0) {
            return false;
        }
        const auto [qc, rc] = divmod_euclid(c, m);
        a += b * qc;
        c = rc;
        n -= 1;
        T y_max = (c * n + d) / m;
        if (s + a * n + b * y_max > 0) {
            return false;
        }
        if (y_max == 0 || (a >= 0 && b >= 0) || (a <= 0 && b <= 0)) {
            return true;
        }
        if (a <= 0) {
            s += a + b;
        }
        n = y_max;
        d = m - d - 1;
        std::swap(a, b);
        std::swap(c, m);
    }
}

/**
 * Δ(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<typename T>
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 g = mp::gcd(D, A);
    T m, r;
    mp::divide_qr(A * B, D, m, r);
    // 解なしをパラメータを用いて判定
    if (r == 0 && D * K + g >= A) {
        return -1;
    }
    T m_neg = -m;
    T threthold = B * K;
    T zero = 0;
    T lo = K + 1, hi;
    if (r == 0) {
        hi = A / g;
    } else {
        chmax(lo, (A * B * K - (A - g) * m) / r + 1);
        hi = (A * B * K) / r + 2;
    }
    // F(lo) <= T, F(hi) > T の不変条件で u_min を二分探索
    while (lo + 1 < hi) {
        T mid = (lo + hi) / 2;
        if (mwf_lr_leq<T>(threthold, lo, mid, A, B, m_neg, D, 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<bigint>(d, a, b, k);
        std::cout << ans << '\n';
    }
}