#764981 (C++23) No.1936 Rational Approximation

提出ソース
結果

問題	No.1936 Rational Approximation
コンテスト
ユーザー	hitonanode
提出日時	2022-05-22 09:56:24
言語	C++23 (gcc 13.3.0 + boost 1.87.0)
結果	AC
実行時間	2 ms / 2,000 ms
コード長	8,054 bytes
コンパイル時間	1,410 ms
コンパイル使用メモリ	170,284 KB
実行使用メモリ	6,944 KB
最終ジャッジ日時	2024-09-20 12:25:42
合計ジャッジ時間	1,918 ms
ジャッジサーバーID （参考情報）	judge4 / judge1
このコードへのチャレンジ
（要ログイン）
ファイルパターン	結果
sample	AC * 1
other	AC * 14
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <deque>
#include <forward_list>
#include <fstream>
#include <functional>
#include <iomanip>
#include <ios>
#include <iostream>
#include <limits>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <string>
#include <tuple>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
using namespace std;
using lint = int64_t;
using pint = pair<int, int>;
using plint = pair<lint, lint>;
struct fast_ios { fast_ios(){ cin.tie(nullptr), ios::sync_with_stdio(false), cout << fixed << setprecision(20); }; } fast_ios_;
#define ALL(x) (x).begin(), (x).end()
#define FOR(i, begin, end) for(int i=(begin),i##_end_=(end);i<i##_end_;i++)
#define IFOR(i, begin, end) for(int i=(end)-1,i##_begin_=(begin);i>=i##_begin_;i--)
#define REP(i, n) FOR(i,0,n)
#define IREP(i, n) IFOR(i,0,n)
template <typename T, typename V>
void ndarray(vector<T>& vec, const V& val, int len) { vec.assign(len, val); }
template <typename T, typename V, typename... Args> void ndarray(vector<T>& vec, const V& val, int len, Args... args) { vec.resize(len), for_each(begin(vec), end(vec), [&](T& v) { ndarray(v, val, args...); }); }
template <typename T> bool chmax(T &m, const T q) { return m < q ? (m = q, true) : false; }
template <typename T> bool chmin(T &m, const T q) { return m > q ? (m = q, true) : false; }
// int floor_lg(__int128_t x) { return x <= 0 ? -1 : 63 - __builtin_clzll(x); }
template <typename T1, typename T2> pair<T1, T2> operator+(const pair<T1, T2> &l, const pair<T1, T2> &r) { return make_pair(l.first + r.first, l.second + r.second); }
template <typename T1, typename T2> pair<T1, T2> operator-(const pair<T1, T2> &l, const pair<T1, T2> &r) { return make_pair(l.first - r.first, l.second - r.second); }
template <typename T> vector<T> sort_unique(vector<T> vec) { sort(vec.begin(), vec.end()), vec.erase(unique(vec.begin(), vec.end()), vec.end()); return vec; }
template <typename T> int arglb(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::lower_bound(v.begin(), v.end(), x)); }
template <typename T> int argub(const std::vector<T> &v, const T &x) { return std::distance(v.begin(), std::upper_bound(v.begin(), v.end(), x)); }
template <typename T> istream &operator>>(istream &is, vector<T> &vec) { for (auto &v : vec) is >> v; return is; }
template <typename T> ostream &operator<<(ostream &os, const vector<T> &vec) { os << '['; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <typename T, size_t sz> ostream &operator<<(ostream &os, const array<T, sz> &arr) { os << '['; for (auto v : arr) os << v << ','; os << ']'; return os; }
#if __cplusplus >= 201703L
template <typename... T> istream &operator>>(istream &is, tuple<T...> &tpl) { std::apply([&is](auto &&... args) { ((is >> args), ...);}, tpl); return is; }
template <typename... T> ostream &operator<<(ostream &os, const tuple<T...> &tpl) { os << '('; std::apply([&os](auto &&... args) { ((os << args << ','), ...);}, tpl); return os << ')'; }
#endif
template <typename T> ostream &operator<<(ostream &os, const deque<T> &vec) { os << "deq["; for (auto v : vec) os << v << ','; os << ']'; return os; }
template <typename T> ostream &operator<<(ostream &os, const set<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T, typename TH> ostream &operator<<(ostream &os, const unordered_set<T, TH> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T> ostream &operator<<(ostream &os, const multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T> ostream &operator<<(ostream &os, const unordered_multiset<T> &vec) { os << '{'; for (auto v : vec) os << v << ','; os << '}'; return os; }
template <typename T1, typename T2> ostream &operator<<(ostream &os, const pair<T1, T2> &pa) { os << '(' << pa.first << ',' << pa.second << ')'; return os; }
template <typename TK, typename TV> ostream &operator<<(ostream &os, const map<TK, TV> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
template <typename TK, typename TV, typename TH> ostream &operator<<(ostream &os, const unordered_map<TK, TV, TH> &mp) { os << '{'; for (auto v : mp) os << v.first << "=>" << v.second << ','; os << '}'; return os; }
#ifdef HITONANODE_LOCAL
const string COLOR_RESET = "\033[0m", BRIGHT_GREEN = "\033[1;32m", BRIGHT_RED = "\033[1;31m", BRIGHT_CYAN = "\033[1;36m", NORMAL_CROSSED = "\033[0;9;37m", RED_BACKGROUND = "\033[1;41m", NORMAL_FAINT = "\033[0;2m";
#define dbg(x) cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << endl
#define dbgif(cond, x) ((cond) ? cerr << BRIGHT_CYAN << #x << COLOR_RESET << " = " << (x) << NORMAL_FAINT << " (L" << __LINE__ << ") " << __FILE__ << COLOR_RESET << endl : cerr)
#else
#define dbg(x) 0
#define dbgif(cond, x) 0
#endif

// \sum_{i=0}^{n-1} floor((ai + b) / m)
// 0 <= n < 2e32
// 1 <= m < 2e32
// 0 <= a, b < m
// Complexity: O(lg(m))
__int128_t floor_sum(__int128_t n, __int128_t m, __int128_t a, __int128_t b) {
    auto safe_mod = [](__int128_t x, __int128_t m) -> __int128_t {
        x %= m;
        if (x < 0) x += m;
        return x;
    };
    auto floor_sum_unsigned = [](__uint128_t n, __uint128_t m, __uint128_t a,
                                 __uint128_t b) -> __uint128_t {
        __uint128_t ans = 0;
        while (true) {
            if (a >= m) {
                ans += n * (n - 1) / 2 * (a / m);
                a %= m;
            }
            if (b >= m) {
                ans += n * (b / m);
                b %= m;
            }

            __uint128_t y_max = a * n + b;
            if (y_max < m) break;
            // y_max < m * (n + 1)
            // floor(y_max / m) <= n
            n = (__uint128_t)(y_max / m);
            b = (__uint128_t)(y_max % m);
            std::swap(m, a);
        }
        return ans;
    };

    __uint128_t ans = 0;
    if (a < 0) {
        __uint128_t a2 = safe_mod(a, m);
        ans -= 1ULL * n * (n - 1) / 2 * ((a2 - a) / m);
        a = a2;
    }
    if (b < 0) {
        __uint128_t b2 = safe_mod(b, m);
        ans -= 1ULL * n * ((b2 - b) / m);
        b = b2;
    }
    return ans + floor_sum_unsigned(n, m, a, b);
}


int main() {
    lint P, Q;
    cin >> P >> Q;

    const __int128_t D = 1000000007LL * 1000000009;

    const lint cnt = floor_sum(Q, Q, P, 0);
    dbg(cnt);

    lint lo = 0, hi = D;
    // lo / D 以下にあるやつは足りない hi / D 以下にあるやつは足りる

    while (hi - lo > 1) {
        lint c = (lo + hi) / 2;
        (floor_sum(Q, D, c, 0) < cnt ? lo : hi) = c;
    }
    dbg(lo);

    // P / Q と lo / D で差分が出るところは？
    lint xsame = 0, xdiff = Q - 1;
    while (xdiff - xsame > 1) {
        const lint x = (xdiff + xsame) / 2;
        if (floor_sum(x + 1, Q, P, 0) - floor_sum(x + 1, D, lo, 0)) {
            xdiff = x;
        } else {
            xsame = x;
        }
    }

    lint ret_l_den = xdiff, ret_l_num = (P * xdiff - 1) / Q;

    lo = 0, hi = D;
    // lo / D 以下にあるやつは以下 hi / D 以下にあるやつは過剰

    while (hi - lo > 1) {
        lint c = (lo + hi) / 2;
        (floor_sum(Q, D, c, 0) <= cnt ? lo : hi) = c;
    }
    dbg(lo);
    xsame = 0, xdiff = Q - 1;

    while (xdiff - xsame > 1) {
        const lint x = (xdiff + xsame) / 2;
        if (floor_sum(x + 1, Q, P, 0) - floor_sum(x + 1, D, hi, 0)) {
            xdiff = x;
        } else {
            xsame = x;
        }
    }

    lint ret_r_den = xdiff, ret_r_num = P * xdiff / Q + 1;

    dbg(make_tuple(ret_l_num, ret_l_den, ret_r_num, ret_r_den));

    cout << ret_l_den + ret_l_num + ret_r_den + ret_r_num << endl;
}
yukicoder

結果

ソースコード
