bigint.hpp

Go to the documentation of this file.

#pragma once
 
#include "core/numeric.hpp"

 
#include <cassert>
#include <cmath>
#include <compare>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
#include "numeric.hpp"
 
// Marks the rarely-taken, out-of-line limb-path helpers cold so the int128 fast
// paths stay small enough to inline into hot callers. GCC and clang honour the
// hint; other compilers (e.g. MSVC on the Boost int128 fallback) ignore it.
#if defined(__GNUC__) || defined(__clang__)
#define PGL_BIGINT_COLD [[gnu::noinline, gnu::cold]]
#else
#define PGL_BIGINT_COLD
#endif
 
namespace pgl {
 
namespace detail {

class LimbStore {
    pgl::int128* data_ = nullptr;   
 
    std::size_t count() const { return static_cast<std::size_t>(data_[-1]); }
 
    void allocate(std::size_t n) {
        pgl::int128* base = new pgl::int128[n + 1];
        base[0] = static_cast<pgl::int128>(n);
        data_ = base + 1;
    }
 
    void copyFrom(const pgl::int128* src, std::size_t n) {
        allocate(n);
        for (std::size_t i = 0; i < n; ++i) {
            data_[i] = src[i];
        }
    }
 
    void release() {
        if (data_) {
            delete[] (data_ - 1);
            data_ = nullptr;
        }
    }
 
public:
    LimbStore() = default;
    LimbStore(const LimbStore& o) {
        if (o.data_) {
            copyFrom(o.data_, o.count());
        }
    }
    LimbStore(LimbStore&& o) noexcept : data_(o.data_) { o.data_ = nullptr; }
    LimbStore& operator=(const LimbStore& o) {
        if (this != &o) {
            release();
            if (o.data_) {
                copyFrom(o.data_, o.count());
            }
        }
        return *this;
    }
    LimbStore& operator=(LimbStore&& o) noexcept {
        if (this != &o) {
            release();
            data_ = o.data_;
            o.data_ = nullptr;
        }
        return *this;
    }
    ~LimbStore() { release(); }

    LimbStore& operator=(const std::vector<pgl::int128>& v) {
        release();
        if (!v.empty()) {
            copyFrom(v.data(), v.size());
        }
        return *this;
    }
 
    bool empty() const { return data_ == nullptr; }
    std::size_t size() const { return data_ ? count() : 0; }
    pgl::int128& operator[](std::size_t i) { return data_[i]; }
    const pgl::int128& operator[](std::size_t i) const { return data_[i]; }
    const pgl::int128* begin() const { return data_; }
    const pgl::int128* end() const { return data_ ? data_ + count() : nullptr; }
    void clear() { release(); }

    void assign(std::size_t n, pgl::int128 val) {
        release();
        if (n) {
            allocate(n);
            for (std::size_t i = 0; i < n; ++i) {
                data_[i] = val;
            }
        }
    }
};
 
}  // namespace detail

class BigInt {
private:
    // alignas(16): on the x86_64-windows-msvc target clang gives __int128 a
    // *preferred* alignment of 16 (used by alignof and by aligned-move codegen,
    // e.g. movaps) but an *ABI* alignment of only 8 (used to lay out aggregates
    // and stack slots). The mismatch means a Point/Rational holding a BigInt can
    // land on an 8-byte-aligned address while the generated code stores small_
    // with movaps, faulting. Pinning the storage to 16 forces the ABI alignment
    // to match the codegen on every toolchain (a no-op where __int128 is already
    // 16-aligned, e.g. the System V ABI, and on the Boost fallback).
    alignas(16) pgl::int128 small_ = 0; 
    detail::LimbStore limbs_;          
    bool negative_ = false;            
 
    using Limbs = std::vector<pgl::int128>;
 
    friend struct std::hash<BigInt>;

    static constexpr int kLimbBits = 62;

    static pgl::int128 base() { return pgl::int128(1) << kLimbBits; }

    static pgl::int128 limbMask() { return base() - 1; }

    static pgl::int128 int128Max() { return pgl::detail::numeric_limits<pgl::int128>::max(); }
 
    // --- magnitude helpers (operate on normalized little-endian limb vectors) ---

    static void trim(Limbs& v) {
        while (!v.empty() && v.back() == 0) {
            v.pop_back();
        }
    }

    static int cmpMag(const Limbs& a, const Limbs& b) {
        if (a.size() != b.size()) {
            return a.size() < b.size() ? -1 : 1;
        }
        for (std::size_t i = a.size(); i-- > 0;) {
            if (a[i] != b[i]) {
                return a[i] < b[i] ? -1 : 1;
            }
        }
        return 0;
    }

    static Limbs addMag(const Limbs& a, const Limbs& b) {
        Limbs r;
        const std::size_t n = a.size() < b.size() ? b.size() : a.size();
        r.reserve(n + 1);
        pgl::int128 carry = 0;
        for (std::size_t i = 0; i < n; ++i) {
            pgl::int128 sum = carry;
            if (i < a.size()) sum += a[i];
            if (i < b.size()) sum += b[i];
            r.push_back(sum & limbMask());
            carry = sum >> kLimbBits;
        }
        if (carry != 0) {
            r.push_back(carry);
        }
        return r;
    }

    static void subMagInPlace(Limbs& a, const Limbs& b) {
        pgl::int128 borrow = 0;
        for (std::size_t i = 0; i < a.size(); ++i) {
            pgl::int128 d = a[i] - borrow - (i < b.size() ? b[i] : pgl::int128(0));
            if (d < 0) {
                d += base();
                borrow = 1;
            } else {
                borrow = 0;
            }
            a[i] = d;
        }
        trim(a);
    }

    static Limbs subMag(const Limbs& a, const Limbs& b) {
        Limbs r = a;
        subMagInPlace(r, b);
        return r;
    }

    static void shiftLeftOneMag(Limbs& v) {
        pgl::int128 carry = 0;
        for (pgl::int128& limb : v) {
            const pgl::int128 next = limb >> (kLimbBits - 1);
            limb = ((limb << 1) & limbMask()) | carry;
            carry = next;
        }
        if (carry != 0) {
            v.push_back(carry);
        }
    }

    static Limbs mulMag(const Limbs& a, const Limbs& b) {
        if (a.empty() || b.empty()) {
            return {};
        }
        Limbs r(a.size() + b.size(), pgl::int128(0));
        for (std::size_t i = 0; i < a.size(); ++i) {
            pgl::int128 carry = 0;
            for (std::size_t j = 0; j < b.size(); ++j) {
                pgl::int128 cur = r[i + j] + a[i] * b[j] + carry;
                r[i + j] = cur & limbMask();
                carry = cur >> kLimbBits;
            }
            r[i + b.size()] += carry;
        }
        trim(r);
        return r;
    }

    static int topBit(pgl::int128 x) {
        int n = 0;
        while (x > 1) {
            x >>= 1;
            ++n;
        }
        return n;
    }

    static std::size_t bitLengthMag(const Limbs& v) {
        if (v.empty()) {
            return 0;
        }
        return (v.size() - 1) * kLimbBits + static_cast<std::size_t>(topBit(v.back())) + 1;
    }

    static bool testBitMag(const Limbs& v, std::size_t i) {
        const std::size_t limb = i / kLimbBits;
        const std::size_t off = i % kLimbBits;
        if (limb >= v.size()) {
            return false;
        }
        return ((v[limb] >> static_cast<int>(off)) & 1) != 0;
    }

    static std::pair<Limbs, Limbs> divmodMag(const Limbs& n, const Limbs& d) {
        if (cmpMag(n, d) < 0) {
            return {Limbs(), n};
        }
        Limbs q(n.size(), pgl::int128(0));
        if (d.size() == 1) {
            const pgl::int128 divisor = d[0];
            pgl::int128 remainder = 0;
            for (std::size_t i = n.size(); i-- > 0;) {
                // remainder < divisor < 2^62, so this is below 2^124.
                const pgl::int128 current = (remainder << kLimbBits) | n[i];
                q[i] = current / divisor;
                remainder = current % divisor;
            }
            trim(q);
            Limbs r;
            if (remainder != 0) {
                r.push_back(remainder);
            }
            return {std::move(q), std::move(r)};
        }
        Limbs r;
        r.reserve(d.size() + 1);
        for (std::size_t bit = bitLengthMag(n); bit-- > 0;) {
            shiftLeftOneMag(r);
            if (testBitMag(n, bit)) {
                if (r.empty()) {
                    r.push_back(pgl::int128(0));
                }
                r[0] |= pgl::int128(1);
            }
            if (cmpMag(r, d) >= 0) {
                subMagInPlace(r, d);
                q[bit / kLimbBits] |= pgl::int128(1) << static_cast<int>(bit % kLimbBits);
            }
        }
        trim(q);
        trim(r);
        return {std::move(q), std::move(r)};
    }
 
    // --- conversions between the small and limb representations ---

    Limbs magToLimbs() const {
        if (!limbs_.empty()) {
            return Limbs(limbs_.begin(), limbs_.end());
        }
        Limbs v;
        pgl::int128 x = small_;
        while (x > 0) {
            v.push_back(x & limbMask());
            x >>= kLimbBits;
        }
        return v;
    }

    static bool limbsFitInt128(const Limbs& v, pgl::int128& out) {
        pgl::int128 acc = 0;
        const pgl::int128 maxv = int128Max();
        for (std::size_t i = v.size(); i-- > 0;) {
            if (acc > maxv / base()) {
                return false;
            }
            acc *= base();
            if (acc > maxv - v[i]) {
                return false;
            }
            acc += v[i];
        }
        out = acc;
        return true;
    }

    void setFromLimbs(Limbs v, bool neg) {
        trim(v);
        if (v.empty()) {
            small_ = 0;
            limbs_.clear();
            negative_ = false;
            return;
        }
        pgl::int128 fitted;
        if (limbsFitInt128(v, fitted)) {
            small_ = fitted;
            limbs_.clear();
        } else {
            small_ = 0;
            limbs_ = std::move(v);
        }
        negative_ = neg;
    }

    static BigInt fromSmall(pgl::int128 magnitude, bool neg) {
        BigInt b;
        b.small_ = magnitude;
        b.negative_ = (magnitude != 0) && neg;
        return b;
    }

    int compareMag(const BigInt& o) const {
        if (limbs_.empty() && o.limbs_.empty()) {
            if (small_ == o.small_) return 0;
            return small_ < o.small_ ? -1 : 1;
        }
        return compareMagGeneral(o);
    }

    PGL_BIGINT_COLD
    int compareMagGeneral(const BigInt& o) const {
        // A limb store only ever holds a magnitude past what an int128 can, so
        // having one at all settles the order against an inline value, and two
        // stores compare by width and then from the top limb down: nothing
        // needs copying out.
        if (limbs_.empty()) {
            return -1;
        }
        if (o.limbs_.empty()) {
            return 1;
        }
        const std::size_t n = limbs_.size(), m = o.limbs_.size();
        if (n != m) {
            return n < m ? -1 : 1;
        }
        for (std::size_t i = n; i-- > 0;) {
            if (limbs_[i] != o.limbs_[i]) {
                return limbs_[i] < o.limbs_[i] ? -1 : 1;
            }
        }
        return 0;
    }

    std::string magToDecimalString() const {
        Limbs v = magToLimbs();
        if (v.empty()) {
            return "0";
        }
        const Limbs divisor = {pgl::int128(int64_t(1000000000000000000))};   // 10^18
        std::vector<int64_t> chunks;
        while (!v.empty()) {
            auto [quotient, remainder] = divmodMag(v, divisor);
            chunks.push_back(remainder.empty() ? int64_t(0) : static_cast<int64_t>(remainder[0]));
            v = std::move(quotient);
        }
        std::string s = std::to_string(chunks.back());
        for (std::size_t i = chunks.size() - 1; i-- > 0;) {
            std::string part = std::to_string(chunks[i]);
            s += std::string(18 - part.size(), '0');
            s += part;
        }
        return s;
    }

    PGL_BIGINT_COLD
    static BigInt addGeneral(const BigInt& a, const BigInt& b) {
        Limbs av = a.magToLimbs();
        Limbs bv = b.magToLimbs();
        BigInt r;
        if (a.negative_ == b.negative_) {
            r.setFromLimbs(addMag(av, bv), a.negative_);
        } else {
            const int c = cmpMag(av, bv);
            if (c > 0) {
                r.setFromLimbs(subMag(av, bv), a.negative_);
            } else if (c < 0) {
                r.setFromLimbs(subMag(bv, av), b.negative_);
            }
            // c == 0 leaves r as the default zero.
        }
        return r;
    }

    static pgl::int128 quotSmall(pgl::int128 a, pgl::int128 b) {
        const std::uint64_t u64max = pgl::detail::numeric_limits<std::uint64_t>::max();
        if (a <= u64max && b <= u64max) {
            return pgl::int128(static_cast<std::uint64_t>(a) / static_cast<std::uint64_t>(b));
        }
        return a / b;
    }

    static pgl::int128 remSmall(pgl::int128 a, pgl::int128 b) {
        const std::uint64_t u64max = pgl::detail::numeric_limits<std::uint64_t>::max();
        if (a <= u64max && b <= u64max) {
            return pgl::int128(static_cast<std::uint64_t>(a) % static_cast<std::uint64_t>(b));
        }
        return a % b;
    }

    PGL_BIGINT_COLD
    static void divmod(const BigInt& a, const BigInt& b, BigInt& q, BigInt& r) {
        if (b.isZero()) {
            throw std::domain_error("pgl::BigInt: division by zero");
        }
        const bool quotientNegative = a.negative_ != b.negative_;
        const bool remainderNegative = a.negative_;   // remainder follows the dividend
        if (a.limbs_.empty() && b.limbs_.empty()) {
            q = fromSmall(quotSmall(a.small_, b.small_), quotientNegative);
            r = fromSmall(remSmall(a.small_, b.small_), remainderNegative);
            return;
        }
        auto [qm, rm] = divmodMag(a.magToLimbs(), b.magToLimbs());
        q.setFromLimbs(std::move(qm), quotientNegative);
        r.setFromLimbs(std::move(rm), remainderNegative);
    }

    PGL_BIGINT_COLD
    pgl::int128 lowBitsInt128() const {
#if defined(__SIZEOF_INT128__)
        __uint128_t m = 0;
        for (std::size_t i = 0; i < limbs_.size() && i < 3; ++i) {
            m += static_cast<__uint128_t>(limbs_[i]) << (kLimbBits * i);
        }
        return static_cast<pgl::int128>(negative_ ? -m : m);
#else
        pgl::int128 m = 0;
        pgl::int128 weight = 1;
        for (std::size_t i = 0; i < limbs_.size() && i < 3; ++i) {
            m += limbs_[i] * weight;
            weight *= base();
        }
        return negative_ ? -m : m;
#endif
    }
 
public:
    BigInt() = default;

    BigInt(pgl::detail::extended_integral auto value) {
        pgl::int128 x = static_cast<pgl::int128>(value);
        if (x == 0) {
            return;
        }
        if (x == pgl::detail::numeric_limits<pgl::int128>::min()) {
            // |min| == 2^127 does not fit in a signed int128, so spell it out.
            negative_ = true;
            limbs_.assign(3, pgl::int128(0));
            limbs_[2] = pgl::int128(1) << 3;   // 8 * 2^124 == 2^127
        } else {
            negative_ = x < 0;
            small_ = x < 0 ? -x : x;
        }
    }

    template <std::floating_point Float>
    explicit BigInt(Float value) {
        assert(std::isfinite(value)
               && "pgl::BigInt: cannot construct from non-finite floating point");
        const bool neg = value < 0;
        const Float mag = std::trunc(neg ? -value : value);  // |value|, fraction dropped
        if (mag == 0) {
            return;   // zero (also covers -0.0 and tiny |value| < 1)
        }
        // Decompose the magnitude exactly as sig * 2^shift, where sig is the
        // integer significand (at most `digits` bits, hence always within an
        // int128) and `shift` places it at the right binary exponent.
        int exponent;
        const Float frac = std::frexp(mag, &exponent);   // mag == frac * 2^exponent
        const int digits = pgl::detail::numeric_limits<Float>::digits;
        const pgl::int128 sig = static_cast<pgl::int128>(std::ldexp(frac, digits));
        const int shift = exponent - digits;             // mag == sig * 2^shift
        if (shift <= 0) {
            // mag < 2^digits <= 2^64, so it fits in a single int128 limb.
            small_ = sig >> (-shift);   // exact: the dropped low bits are zero
            negative_ = neg;
            return;
        }
        // shift > 0: lay sig down at bit offset `shift` in the base-2^62 store.
        // sig < 2^digits and the in-limb offset is < 62, so the shifted chunk
        // still fits in an int128 before it is split across limbs.
        Limbs v(static_cast<std::size_t>(shift / kLimbBits), pgl::int128(0));
        pgl::int128 x = sig << (shift % kLimbBits);
        while (x > 0) {
            v.push_back(x & limbMask());
            x >>= kLimbBits;
        }
        setFromLimbs(std::move(v), neg);
    }
 
    // --- inspectors ---

    bool isZero() const { return limbs_.empty() && small_ == 0; }

    bool isOne() const { return limbs_.empty() && !negative_ && small_ == 1; }

    bool isNegative() const { return negative_; }

    bool fitsInt128() const { return limbs_.empty(); }

    bool fitsInt64() const {
        return limbs_.empty() && small_ < (pgl::int128(1) << 63);
    }

    bool fitsLimbs(std::size_t limbs) const { return limbs_.size() <= limbs; }

    int sign() const {
        if (isZero()) return 0;
        return negative_ ? -1 : 1;
    }

    BigInt abs() const {
        BigInt r = *this;
        r.negative_ = false;
        return r;
    }
 
    // --- conversions ---

    explicit operator pgl::int128() const {
        if (limbs_.empty()) {
            return negative_ ? -small_ : small_;
        }
        return lowBitsInt128();
    }

    template <std::signed_integral T>
    explicit operator T() const {
        return static_cast<T>(static_cast<pgl::int128>(*this));
    }
 
    explicit operator bool() const { return !isZero(); }

    explicit operator long double() const {
        long double d = 0;
        if (limbs_.empty()) {
            d = static_cast<long double>(small_);
        } else {
            const long double b = static_cast<long double>(base());
            for (std::size_t i = limbs_.size(); i-- > 0;) {
                d = d * b + static_cast<long double>(limbs_[i]);
            }
        }
        return negative_ ? -d : d;
    }
    explicit operator double() const { return static_cast<double>(static_cast<long double>(*this)); }
    explicit operator float() const { return static_cast<float>(static_cast<long double>(*this)); }
 
    // --- arithmetic ---
 
    BigInt operator-() const {
        BigInt r = *this;
        if (!r.isZero()) {
            r.negative_ = !r.negative_;
        }
        return r;
    }
 
    BigInt operator+() const { return *this; }
 
    friend BigInt operator+(const BigInt& a, const BigInt& b) {
        if (a.limbs_.empty() && b.limbs_.empty()) {
            const pgl::int128 ma = a.small_;
            const pgl::int128 mb = b.small_;
            if (a.negative_ == b.negative_) {
                if (ma <= int128Max() - mb) {
                    return fromSmall(ma + mb, a.negative_);
                }
                // magnitude overflow: fall through to the limb path
            } else if (ma >= mb) {
                return fromSmall(ma - mb, a.negative_);
            } else {
                return fromSmall(mb - ma, b.negative_);
            }
        }
        return addGeneral(a, b);
    }
 
    friend BigInt operator-(const BigInt& a, const BigInt& b) {
        if (a.limbs_.empty() && b.limbs_.empty()) {
            const pgl::int128 ma = a.small_;
            const pgl::int128 mb = b.small_;
            // a - b == a + (-b); replicate operator+'s small path with b's sign
            // flipped, so we avoid materialising the negated temporary -b.
            if (a.negative_ != b.negative_) {
                if (ma <= int128Max() - mb) {
                    return fromSmall(ma + mb, a.negative_);
                }
            } else if (ma >= mb) {
                return fromSmall(ma - mb, a.negative_);
            } else {
                return fromSmall(mb - ma, !a.negative_);
            }
        }
        return subGeneral(a, b);
    }

    PGL_BIGINT_COLD
    static BigInt subGeneral(const BigInt& a, const BigInt& b) { return a + (-b); }
 
    friend BigInt operator*(const BigInt& a, const BigInt& b) {
        const bool neg = a.negative_ != b.negative_;
        if (a.limbs_.empty() && b.limbs_.empty()) {
            const pgl::int128 ma = a.small_;
            const pgl::int128 mb = b.small_;
#if defined(__SIZEOF_INT128__)
            // Native int128: both magnitudes are non-negative, so a signed-overflow
            // check on the product is exactly the "does it still fit" test we want,
            // and __builtin_mul_overflow performs it without a 128-bit division.
            pgl::int128 prod;
            if (!__builtin_mul_overflow(ma, mb, &prod)) {
                return fromSmall(prod, neg);
            }
#else
            // Boost fallback (no native int128, no overflow builtin): guard the
            // product with a division, matching operator+'s overflow style.
            if (ma == 0 || mb <= int128Max() / ma) {
                return fromSmall(ma * mb, neg);
            }
#endif
            // magnitude overflow: fall through to the limb path
        }
        return mulGeneral(a, b, neg);
    }

    PGL_BIGINT_COLD
    static BigInt mulGeneral(const BigInt& a, const BigInt& b, bool neg) {
        BigInt out;
        out.setFromLimbs(mulMag(a.magToLimbs(), b.magToLimbs()), neg);
        return out;
    }
 
    friend BigInt operator/(const BigInt& a, const BigInt& b) {
        if (a.limbs_.empty() && b.limbs_.empty()) {
            if (b.small_ == 0) {
                throw std::domain_error("pgl::BigInt: division by zero");
            }
            // The quotient is all the caller wants, so skip computing a remainder.
            return fromSmall(quotSmall(a.small_, b.small_), a.negative_ != b.negative_);
        }
        BigInt q;
        BigInt r;
        divmod(a, b, q, r);
        return q;
    }
 
    friend BigInt operator%(const BigInt& a, const BigInt& b) {
        if (a.limbs_.empty() && b.limbs_.empty()) {
            if (b.small_ == 0) {
                throw std::domain_error("pgl::BigInt: division by zero");
            }
            // The remainder is all the caller wants, so skip computing a quotient.
            return fromSmall(remSmall(a.small_, b.small_), a.negative_);
        }
        BigInt q;
        BigInt r;
        divmod(a, b, q, r);
        return r;
    }
 
    BigInt& operator+=(const BigInt& o) { return *this = *this + o; }
    BigInt& operator-=(const BigInt& o) { return *this = *this - o; }
    BigInt& operator*=(const BigInt& o) { return *this = *this * o; }
    BigInt& operator/=(const BigInt& o) { return *this = *this / o; }
    BigInt& operator%=(const BigInt& o) { return *this = *this % o; }
 
    BigInt& operator++() { return *this += BigInt(1); }
    BigInt operator++(int) { BigInt tmp = *this; ++(*this); return tmp; }
    BigInt& operator--() { return *this -= BigInt(1); }
    BigInt operator--(int) { BigInt tmp = *this; --(*this); return tmp; }
 
    // --- comparison ---
 
    std::strong_ordering operator<=>(const BigInt& o) const {
        if (negative_ != o.negative_) {
            return negative_ ? std::strong_ordering::less : std::strong_ordering::greater;
        }
        int c = compareMag(o);
        if (negative_) {
            c = -c;   // among negatives, the larger magnitude is the smaller value
        }
        if (c < 0) return std::strong_ordering::less;
        if (c > 0) return std::strong_ordering::greater;
        return std::strong_ordering::equal;
    }
 
    bool operator==(const BigInt& o) const {
        return negative_ == o.negative_ && compareMag(o) == 0;
    }
 
    // --- comparison against floating point ---
    //
    // These are exact for every finite magnitude: rather than casting *this to a
    // (lossy) double, the float is split at its truncation point with the same
    // round-toward-zero rule the floating-point constructor uses, and the dropped
    // fractional part breaks any tie on the integer part. Being members, they
    // also cover the reversed forms (`f == b`, `f < b`, ...) via C++20 rewriting.

    template <std::floating_point Float>
    bool operator==(Float f) const {
        if (!std::isfinite(f) || f != std::trunc(f)) {
            return false;   // NaN/inf or a fractional value can never equal an integer
        }
        return *this == BigInt(f);
    }

    template <std::floating_point Float>
    std::partial_ordering operator<=>(Float f) const {
        if (std::isnan(f)) {
            return std::partial_ordering::unordered;
        }
        if (std::isinf(f)) {
            return f > 0 ? std::partial_ordering::less : std::partial_ordering::greater;
        }
        // Compare against the truncated integer part first; if they differ, that
        // settles it. BigInt(f) discards f's fraction exactly as std::trunc does.
        const Float t = std::trunc(f);
        if (std::strong_ordering c = *this <=> BigInt(f); c != 0) {
            return c;
        }
        // Integer parts are equal, so the leftover fraction decides: a positive
        // fraction makes f larger than this integer, a negative one smaller.
        if (f == t) {
            return std::partial_ordering::equivalent;
        }
        return f > t ? std::partial_ordering::less : std::partial_ordering::greater;
    }
 
    // --- stream I/O ---
 
    friend std::ostream& operator<<(std::ostream& os, const BigInt& b) {
        if (b.negative_) {
            os << '-';
        }
        if (b.limbs_.empty()) {
            os << b.small_;   // uses the pgl::int128 stream operator
        } else {
            os << b.magToDecimalString();
        }
        return os;
    }
 
    friend std::istream& operator>>(std::istream& is, BigInt& b) {
        std::istream::sentry sentry(is);
        if (!sentry) {
            return is;
        }
        bool neg = false;
        int c = is.peek();
        if (c == '+' || c == '-') {
            neg = c == '-';
            is.get();
            c = is.peek();
        }
        if (c < '0' || c > '9') {
            is.setstate(std::ios::failbit);
            return is;
        }
        BigInt result;
        const BigInt ten(10);
        for (; c >= '0' && c <= '9'; c = is.peek()) {
            is.get();
            result = result * ten + BigInt(c - '0');
        }
        b = neg ? -result : result;
        return is;
    }
};

inline BigInt abs(const BigInt& v) { return v.abs(); }
 
}  // namespace pgl
 
#undef PGL_BIGINT_COLD