結果

問題 No.3398 Accuracy of Integer Division Approximate Function 2
コンテスト
ユーザー 👑 Mizar
提出日時 2025-11-02 15:14:29
言語 C++17
(gcc 13.3.0 + boost 1.89.0)
結果
AC  
実行時間 98 ms / 2,000 ms
コード長 6,621 bytes
記録
コンパイル時間 4,875 ms
コンパイル使用メモリ 385,644 KB
実行使用メモリ 7,720 KB
最終ジャッジ日時 2025-12-04 23:30:26
合計ジャッジ時間 6,380 ms
ジャッジサーバーID
（参考情報） 		judge2 / judge1
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 20
権限があれば一括ダウンロードができます

ソースコード

diff #
raw source code 

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

// 商を無限大方向に切り下げる除算(剰余が0でない時の符号は除数の符号に一致)
template <typename T>
constexpr std::pair<T, T> floor_divmod(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 <typename T>
constexpr std::pair<T, T> euclid_divmod(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};
}

/**
    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)
*/
bigint mwf(bigint n, bigint m, bigint a, bigint b, bigint c, bigint d) {
    assert(n > 0 && m > 0);
    bigint sum_acc = 0, max_acc = b * euclid_divmod(d, m).first;
    while(true) {
        auto [q1, r1] = euclid_divmod(c, m);
        c = r1;
        a = a + b * q1;
        auto [q2, r2] = euclid_divmod(d, m);
        d = r2;
        sum_acc = sum_acc + b * q2;
        assert(0 <= c && c < m && 0 <= d && d < m);
        max_acc = mp::max(max_acc, sum_acc);
        auto y_max = (c * (n - 1) + d) / m;
        if(y_max == 0) {
            return mp::max(max_acc, sum_acc + a * (n - 1));
        }
        if(a >= 0) {
            max_acc = mp::max(max_acc, sum_acc + a * (n - 1) + b * y_max);
        } else {
            sum_acc = sum_acc + a + b;
        }
        n = y_max;
        d = m - d - 1;
        std::swap(a, b);
        std::swap(c, m);
    }
}

/**
 * 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 bigint
 * @param n bigint
 * @param m bigint
 * @param a bigint
 * @param b bigint
 * @param c bigint
 * @param d bigint
 * @return bool
 */
bool mwf_leq(bigint z, bigint n, bigint m, bigint a, bigint b, bigint c, bigint d) {
    assert(n > 0 && m > 0);
    bigint sum_acc = -z;
    while (true) {
        auto [q1, r1] = euclid_divmod(c, m);
        c = r1;
        a = a + b * q1;
        auto [q2, r2] = euclid_divmod(d, m);
        d = r2;
        sum_acc = sum_acc + b * q2;
        assert(0 <= c && c < m && 0 <= d && d < m);
        if (sum_acc > 0) {
            return false;
        }
        auto y_max = (c * (n - 1) + d) / m;
        if (y_max == 0) {
            return sum_acc + a * (n - 1) <= 0;
        }
        if (sum_acc + mp::max(a * (n - 1), bigint(0)) + mp::max(b * y_max, bigint(0)) <= 0) {
            return true;
        }
        if (a >= 0) {
            if (sum_acc + a * (n - 1) + b * y_max > 0) {
                return false;
            }
        } else {
            sum_acc = sum_acc + a + b;
        }
        n = y_max;
        d = m - d - 1;
        std::swap(a, b);
        std::swap(c, m);
    }
}

bigint mwf_lr(bigint l, bigint r, bigint m, bigint a, bigint b, bigint c, bigint d) {
    assert(l < r && m > 0);
    bigint n = r - l;
    auto [quot, rem] = euclid_divmod(bigint(c * l + d), m);
    return a * l + b * quot + mwf(n, m, a, b, c, rem);
}

bool mwf_lr_leq(bigint z, bigint l, bigint r, bigint m, bigint a, bigint b, bigint c, bigint d) {
    assert(l < r && m > 0);
    bigint n = r - l;
    auto [quot, rem] = euclid_divmod(bigint(c * l + d), m);
    return mwf_leq(z - a * l - b * quot, n, m, a, b, c, rem);
}

/**
 * x_min(D, A, B, K) を半開区間二分探索 [0, A'BK+2) で求めます（解なしは -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 bigint
 * @param A bigint
 * @param B bigint
 * @param K bigint
 * @return bigint 
 */
bigint compute_xmin_leq(bigint D, bigint A, bigint B, bigint K) {
    assert(D > 0 && A > 0 && B > 0);
    if (K < 0) {
        return 0;
    }
    bigint gcd_DA = mp::gcd(D, A);
    bigint Dred = D / gcd_DA, Ared = A / gcd_DA;
    bigint Mred = (Ared * B) / Dred, Rred = (Ared * B) % Dred;
    bigint Tdelta = B * K;
    // 解なしをパラメータを用いて判定
    if (Rred == 0 && Dred * K + 1 >= Ared) {
        return -1;
    }
    // [0, hi) の半開区間、緩い上界 A'BK+1 を包括する hi = A'BK+2 を設定
    bigint lo = 0, hi = Ared * B * K + 2;
    // F(hi) > T の不変条件を確認
    assert(!mwf_lr_leq(Tdelta, lo, hi, Ared, B, -Mred, Dred, 0));
    // F(lo) <= T, F(hi) > T の不変条件で u_min を二分探索
    while (lo + 1 < hi) {
        bigint mid = (lo + hi) / 2;
        if (mwf_lr_leq(Tdelta, lo, mid, Ared, B, -Mred, Dred, 0)) {
            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';
    }
}
0