numeric.hpp

Go to the documentation of this file.

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

 
#include <array>
#include <cmath>
#include <compare>
#include <concepts>
#include <cstdint>
#include <istream>
#include <limits>
#include <numeric>
#include <ostream>
#include <type_traits>
 
#if !defined(__SIZEOF_INT128__)
// Toolchains without the 128-bit extension (e.g. MSVC) fall back to Boost.
// Boost is only pulled in here, never when __int128_t is available.
#include <boost/multiprecision/cpp_int.hpp>
#endif
 
namespace pgl {

template <class T>
class Rational;

class BigInt;

#if defined(__SIZEOF_INT128__)
using int128 = __int128_t;
#else
// Boost's own `int128_t` is sign-magnitude with a *128-bit magnitude*, so its
// max() is 2^128 - 1 — wider than native __int128 (two's complement, max
// 2^127 - 1). Using it would silently give pgl::int128 a different range on the
// fallback than on the native path, breaking the shared overflow guards and
// invariants. A 127-bit magnitude reproduces native's max() exactly; the only
// difference is min() is -(2^127 - 1) rather than -2^127, an extreme value the
// exact-arithmetic intermediates never depend on.
using int128 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<
    127, 127, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
#endif
 
#if defined(__SIZEOF_INT128__)
inline std::ostream& operator<<(std::ostream& stream, const int128& value) {
    // The magnitude is taken unsigned: the most negative value has no negation
    // in the signed type, and forming one there is undefined behavior that
    // optimizers exploit -- a `value == -value` guard for it folds away.
    __uint128_t magnitude = value < 0 ? -static_cast<__uint128_t>(value)
                                      : static_cast<__uint128_t>(value);
    char digits[41];   // 39 decimal digits, a sign and a terminator
    char* first = digits + sizeof(digits) - 1;
    *first = '\0';
    do {
        *--first = static_cast<char>('0' + static_cast<int>(magnitude % 10));
        magnitude /= 10;
    } while (magnitude != 0);
    if (value < 0) {
        *--first = '-';
    }
    return stream << first;
}

inline std::istream& operator>>(std::istream& stream, int128& value) {
    std::istream::sentry sentry(stream);
    if (!sentry) {
        return stream;
    }
    bool neg = false;
    int c = stream.peek();
    if (c == '+' || c == '-') {
        neg = c == '-';
        stream.get();
        c = stream.peek();
    }
    if (c < '0' || c > '9') {
        stream.setstate(std::ios::failbit);
        return stream;
    }
    // Accumulated unsigned, where the most negative value's magnitude exists
    // and wrapping is defined; a magnitude past the type wraps, as before.
    __uint128_t result = 0;
    for (; c >= '0' && c <= '9'; c = stream.peek()) {
        stream.get();
        result = result * 10 + static_cast<unsigned>(c - '0');
    }
    value = static_cast<int128>(neg ? -result : result);
    return stream;
}
#endif
 
}  // namespace pgl
 
#if !defined(__SIZEOF_INT128__)
// The Boost.Multiprecision fallback for pgl::int128 (used when the native
// __int128 extension is unavailable, e.g. MSVC) deliberately omits mixed
// operations with floating-point types, whereas native __int128 supports them
// through the usual arithmetic conversions. These free operators restore that
// double interop so pgl::int128 stays a drop-in replacement: each one converts
// the 128-bit value to double and operates in double, exactly as the built-in
// path would.
//
// They are defined in namespace boost::multiprecision (the namespace of
// pgl::int128's underlying type) so that argument-dependent lookup finds them
// everywhere -- crucially inside the `requires(Int x, Float g){ x * g; }` probe
// in Rational's float constructor, where MSVC only consults ADL, not ordinary
// lookup. A single operator<=> covers all four relational operators in both
// argument orders, and a single operator== covers ==/!=, avoiding the C++20
// reversed-candidate ambiguity that defining both directions would create.
namespace boost::multiprecision {
inline std::partial_ordering operator<=>(const pgl::int128& a, double b) {
    return static_cast<double>(a) <=> b;
}
inline bool operator==(const pgl::int128& a, double b) {
    return static_cast<double>(a) == b;
}
inline double operator+(const pgl::int128& a, double b) { return static_cast<double>(a) + b; }
inline double operator+(double a, const pgl::int128& b) { return a + static_cast<double>(b); }
inline double operator-(const pgl::int128& a, double b) { return static_cast<double>(a) - b; }
inline double operator-(double a, const pgl::int128& b) { return a - static_cast<double>(b); }
inline double operator*(const pgl::int128& a, double b) { return static_cast<double>(a) * b; }
inline double operator*(double a, const pgl::int128& b) { return a * static_cast<double>(b); }
inline double operator/(const pgl::int128& a, double b) { return static_cast<double>(a) / b; }
inline double operator/(double a, const pgl::int128& b) { return a / static_cast<double>(b); }
}  // namespace boost::multiprecision
#endif
 
namespace pgl::detail {

template <class T>
struct numeric_limits : std::numeric_limits<T> {};
 
#if defined(__SIZEOF_INT128__)
template <>
struct numeric_limits<__int128_t> {
    static constexpr bool is_specialized = true;
    static constexpr bool is_signed = true;
    static constexpr bool is_integer = true;
    static constexpr bool is_exact = true;
    static constexpr bool is_bounded = true;
    static constexpr int digits = 127;     // value bits of a signed 128-bit int
    static constexpr int digits10 = 38;
    static constexpr __int128_t min() noexcept {
        return static_cast<__int128_t>(static_cast<__uint128_t>(1) << 127);
    }
    static constexpr __int128_t lowest() noexcept { return min(); }
    static constexpr __int128_t max() noexcept {
        return static_cast<__int128_t>(~(static_cast<__uint128_t>(1) << 127));
    }
};
#endif
 
template <typename T>
concept extended_integral =
    std::integral<T>
    || std::same_as<std::remove_cv_t<T>, pgl::int128>;

template <class A, class B>
constexpr std::strong_ordering strongOrder(const A& a, const B& b) {
    if constexpr (requires { std::strong_order(a, b); }) {
        if constexpr (std::is_floating_point_v<A> || std::is_floating_point_v<B>) {
            if (a == b) {
                return std::strong_ordering::equal;
            }
        }
        return std::strong_order(a, b);
    } else {
        if (a < b) return std::strong_ordering::less;
        if (b < a) return std::strong_ordering::greater;
        return std::strong_ordering::equal;
    }
}

template <class A, class B>
constexpr auto threeWay(const A& a, const B& b) {
    if constexpr (requires { a <=> b; }) {
        return a <=> b;
    } else {
        return a < b ? std::strong_ordering::less
             : b < a ? std::strong_ordering::greater
                     : std::strong_ordering::equal;
    }
}

template <typename T>
struct _promote {
    using type = T;
};
 
#ifndef PGL_DISABLE_PROMOTION
// Rational is deliberately never promoted: it manages its own overflow by
// reducing to lowest terms, so the storage type stays as the user chose it.
 
// Signed built-in integers promote by width rather than by fixed-width alias.
// Which of `int`, `long` and `long long` each of `int32_t` and `int64_t` names
// is platform-dependent, and the three are distinct types whatever their
// widths: specializing on the aliases alone leaves whichever type they skip —
// `long long` under LP64, `long` under LLP64 — matching nothing but the
// identity primary template. Nothing about that is visible at the use site.
// The coordinate still compiles everywhere it should, and every overflow guard
// built on this trait is simply absent, so an intermediate meant to be
// evaluated 128 bits wide wraps in 64 instead. Worse than a plain wrap, once a
// wrapped value is widened again the modular arithmetic stops agreeing with the
// exact arithmetic around it, which is how a wrong sign reaches a predicate.
template <typename T>
    requires std::signed_integral<T>
struct _promote<T> {
    using type =
        std::conditional_t<sizeof(T) <= 1, int16_t,
        std::conditional_t<sizeof(T) <= 2, int32_t,
        std::conditional_t<sizeof(T) <= 4, int64_t,
        std::conditional_t<sizeof(T) <= 8, pgl::int128, pgl::BigInt>>>>;
};
 
// Kept as an explicit specialization, which outranks the partial one above on
// the platforms where pgl::int128 is the native extension and the standard
// library counts it as an integer.
template <>
struct _promote<pgl::int128> {
    using type = pgl::BigInt;
};
 
template <>
struct _promote<float> {
    using type = double;
};
 
template <>
struct _promote<double> {
    using type = long double;
};
#endif
 
template <typename T>
using promoted_number_t = typename _promote<T>::type;

template <class Target, class Number>
constexpr decltype(auto) asNumber(const Number& value) {
    if constexpr (std::is_same_v<Target, Number>) {
        return (value);
    } else {
        return static_cast<Target>(value);
    }
}

template <class T>
inline constexpr bool arbitraryPrecision = false;
 
template <>
inline constexpr bool arbitraryPrecision<pgl::BigInt> = true;
 
template <class Int>
inline constexpr bool arbitraryPrecision<pgl::Rational<Int>> = arbitraryPrecision<Int>;

template <class ResultNumber>
using floating_result_t =
    std::conditional_t<std::is_floating_point_v<ResultNumber>, ResultNumber, double>;

constexpr auto gcd(std::integral auto a, std::integral auto b) {
    return std::gcd(a, b);
}

inline auto gcd(auto a, auto b) {
    if (a == b) {
        return a;
    }
 
    while (b != 0) {
        const auto remainder = a % b;
        a = b;
        b = remainder;
    }
 
    return a;
}

template <class Integer>
inline std::array<Integer, 3> extendedGcd(Integer a, Integer b) {
    Integer remainder = a, currentRemainder = b;
    Integer coefficientA(1), nextCoefficientA(0);
    Integer coefficientB(0), nextCoefficientB(1);
    while (currentRemainder != Integer(0)) {
        const Integer quotient = remainder / currentRemainder;
        Integer next = remainder - quotient * currentRemainder;
        remainder = currentRemainder;
        currentRemainder = next;
        next = coefficientA - quotient * nextCoefficientA;
        coefficientA = nextCoefficientA;
        nextCoefficientA = next;
        next = coefficientB - quotient * nextCoefficientB;
        coefficientB = nextCoefficientB;
        nextCoefficientB = next;
    }
    if (remainder < Integer(0)) {
        remainder = -remainder;
        coefficientA = -coefficientA;
        coefficientB = -coefficientB;
    }
    return {remainder, coefficientA, coefficientB};
}

template <class Integer>
inline Integer nearestQuotient(const Integer& numerator, const Integer& divisor) {
    Integer quotient = numerator / divisor;          // truncates toward zero
    Integer remainder = numerator - quotient * divisor;
    if (remainder < Integer(0)) {                    // make it a floor division
        quotient -= Integer(1);
        remainder += divisor;
    }
    if (Integer(2) * remainder >= divisor) {
        quotient += Integer(1);
    }
    return quotient;
}

template <class Target, class Integer>
inline bool representableAs(const Integer& value) {
    if constexpr (!numeric_limits<Target>::is_bounded) {
        return true;
    } else {
        return value >= Integer(numeric_limits<Target>::lowest())
            && value <= Integer(numeric_limits<Target>::max());
    }
}

template <class Target, class Integer>
inline Target narrowTo(const Integer& value) {
    if constexpr (std::same_as<Target, Integer>) {
        return value;
    } else {
        return static_cast<Target>(static_cast<pgl::int128>(value));
    }
}

constexpr auto abs(extended_integral auto value) {
    using T = decltype(value);
    if constexpr (std::is_unsigned_v<T>) {
        return value;
    } else {
        return value < 0 ? -value : value;
    }
}

constexpr auto abs(std::floating_point auto value) {
    return std::abs(value);
}

inline constexpr auto abs(auto value) {
    return value >= 0 ? value : -value;
}

template <typename Int>
constexpr pgl::Rational<Int> abs(pgl::Rational<Int> value) {
    // The parameter is a copy, so simplify it in place and read both parts from
    // the reduced form; numerator() and denominator() would each reduce it anew.
    value.simplify();
    return pgl::Rational(pgl::detail::abs(value.numerator()), value.denominator(), true);
}

template <class To, class From>
[[nodiscard]] constexpr To convertCoordinate(From value) {
    return static_cast<To>(value);
}
 
}  // namespace pgl::detail