real_polyhedral_homotopy_continuation/shared_module/math/rational.hpp

/* basically copied from boost's */
#pragma once

#include <limits>
#include <numeric>
#include <stdexcept>
#include <type_traits>

namespace detail
{
template <typename __from, typename __to>
static constexpr bool is_compatible_integer_v =
    !std::is_array_v<__from>                                                                                     //
        && std::numeric_limits<__from>::is_specialized                                                           //
        && std::numeric_limits<__from>::is_integer                                                               //
        && (std::numeric_limits<__from>::digits <= std::numeric_limits<__to>::digits)                            //
        && (std::numeric_limits<__from>::radix == std::numeric_limits<__to>::radix)                              //
        && ((std::numeric_limits<__from>::is_signed == false) || (std::numeric_limits<__to>::is_signed == true)) //
        && std::is_convertible_v<__from, __to>                                                                   //
    || std::is_same_v<__from, __to>;

template <typename __from, typename __to>
static constexpr bool is_backward_compatible_integer_v =
    !std::is_array_v<__from>                                                    //
    && std::numeric_limits<__from>::is_specialized                              //
    && std::numeric_limits<__from>::is_integer                                  //
    && !is_compatible_integer_v<__from, __to>                                   //
    && (std::numeric_limits<__from>::radix == std::numeric_limits<__to>::radix) //
    && std::is_convertible_v<__from, __to>;

template <typename T>
static constexpr T abs(T const& x)
{
    return x < T{0} ? -x : x;
}
} // namespace detail

struct bad_rational : public std::domain_error {
    using std::domain_error::domain_error;

    explicit bad_rational() : std::domain_error("bad rational: zero denom") {}
};

template <typename __int_type>
struct rational {
    static_assert(std::numeric_limits<__int_type>::is_specialized);

    constexpr rational() = default;

    template <typename T>
    constexpr rational(__int_type _num)
    {
        if constexpr (detail::is_compatible_integer_v<T, __int_type>) {
            num = _num;
            den = 1;
        } else if constexpr (detail::is_backward_compatible_integer_v<T, __int_type>) {
            assign(_num, static_cast<T>(1));
        }
    }

    template <typename T, typename U>
    constexpr rational(T _num, U _den)
    {
        if constexpr (detail::is_compatible_integer_v<T, __int_type> && detail::is_compatible_integer_v<U, __int_type>) {
            num = _num;
            den = _den;
            normalize();
        } else if constexpr (detail::is_backward_compatible_integer_v<T, __int_type>
                             && detail::is_backward_compatible_integer_v<U, __int_type>) {
            assign(_num, _den);
        }
    }

    template <typename __new_int_type>
    explicit constexpr rational(rational<__new_int_type> const& other) : num(other.numerator())
    {
        if (is_normalized(__int_type(other.numerator()), __int_type(other.denominator()))) {
            if constexpr (!detail::is_compatible_integer_v<__new_int_type, __int_type>) {
                if (is_safe_narrowing_conversion(other.denominator()) && is_safe_narrowing_conversion(other.numerator())) {
                    den = other.denominator();
                }
            } else
                den = other.denominator();
        } else
            throw bad_rational("bad rational: denormalized conversion");
    }

    template <typename T, typename U>
    rational& assign(const T& n, const U& d)
    {
        if constexpr (!detail::is_compatible_integer_v<T, __int_type>
                      || !detail::is_compatible_integer_v<U, __int_type>
                             && detail::is_backward_compatible_integer_v<T, __int_type>
                             && detail::is_backward_compatible_integer_v<U, __int_type>)
            if (!is_safe_narrowing_conversion(n) || !is_safe_narrowing_conversion(d)) throw bad_rational();
        return *this = rational<__int_type>(static_cast<__int_type>(n), static_cast<__int_type>(d));
    }

    constexpr const __int_type& numerator() const noexcept { return num; }

    constexpr const __int_type& denominator() const noexcept { return den; }

    constexpr rational operator+() const noexcept { return *this; }

    constexpr rational operator-() const noexcept { return rational(static_cast<__int_type>(-num), den); }

    template <typename T, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
    constexpr rational& operator+=(const T& n)
    {
        num += n * den;
        return *this;
    }

    template <typename T, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
    constexpr rational& operator-=(const T& n)
    {
        num -= n * den;
        return *this;
    }

    template <typename T, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
    constexpr rational& operator*=(const T& n)
    {
        const auto gcd  = std::gcd(static_cast<__int_type>(n), den);
        num            *= n / gcd;
        den            /= gcd;
        return *this;
    }

    template <typename T, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
    constexpr rational& operator/=(const T& n)
    {
        if (n == zero) throw bad_rational();
        if (num == zero) return *this;

        const auto gcd  = std::gcd(num, static_cast<__int_type>(n));
        num            /= gcd;
        den            *= n / gcd;

        if (den < zero) {
            num = -num;
            den = -den;
        }

        return *this;
    }

    constexpr rational& operator+=(const rational& other)
    {
        const auto other_num = other.numerator();
        const auto other_den = other.denominator();

        auto gcd  = std::gcd(den, other_den);
        den      /= gcd;
        num       = num * (other_den / gcd) + other_num * den;
        gcd       = std::gcd(num, gcd);
        num      /= gcd;
        den      *= other_den / gcd;

        return *this;
    }

    constexpr rational& operator-=(const rational& other)
    {
        const auto other_num = other.numerator();
        const auto other_den = other.denominator();

        auto gcd  = std::gcd(den, other_den);
        den      /= gcd;
        num       = num * (other_den / gcd) - other_num * den;
        gcd       = std::gcd(num, gcd);
        num      /= gcd;
        den      *= other_den / gcd;

        return *this;
    }

    constexpr rational& operator*=(const rational& other)
    {
        const auto other_num = other.numerator();
        const auto other_den = other.denominator();

        const auto gcd1 = std::gcd(num, other_den);
        const auto gcd2 = std::gcd(other_num, den);
        num             = num / gcd1 * other_num / gcd2;
        den             = den / gcd2 * other_den / gcd1;

        return *this;
    }

    constexpr rational& operator/=(const rational& other)
    {
        const auto other_num = other.numerator();
        const auto other_den = other.denominator();

        if (other_num == zero) throw bad_rational();
        if (num == zero) return *this;

        const auto gcd1 = std::gcd(num, other_num);
        const auto gcd2 = std::gcd(den, other_den);
        num             = num / gcd1 * other_den / gcd2;
        den             = den / gcd2 * other_num / gcd1;

        if (den < zero) {
            num = -num;
            den = -den;
        }

        return *this;
    }

    constexpr const rational& operator++()
    {
        num += den;
        return *this;
    }

    constexpr const rational& operator--()
    {
        num -= den;
        return *this;
    }

    constexpr rational operator++(int)
    {
        rational tmp(*this);
        num += den;
        return tmp;
    }

    constexpr rational operator--(int)
    {
        rational tmp(*this);
        num -= den;
        return tmp;
    }

    template <typename T, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
    rational& operator=(const T& n)
    {
        if constexpr (detail::is_compatible_integer_v<T, __int_type>)
            return assign(static_cast<__int_type>(n), __int_type{1});
        else if constexpr (detail::is_backward_compatible_integer_v<T, __int_type>)
            return assign(n, static_cast<T>(1));
    }

    constexpr bool operator<(const rational& other) const noexcept
    {
        assert(den > zero);
        assert(other.den > zero);

        struct {
            __int_type n, d, q, r;
        } ts = {num, den, static_cast<__int_type>(num / den), static_cast<__int_type>(num % den)},
          os = {other.num,
                other.den,
                static_cast<__int_type>(other.num / other.den),
                static_cast<__int_type>(other.num % other.den)};

        unsigned reverse{0u};

        // apply Eculidean GCD algorithm
        while (ts.r < zero) {
            ts.r += ts.d;
            --ts.q;
        }
        while (os.r < zero) {
            os.r += os.d;
            --os.q;
        }

        while (true) {
            if (ts.q != os.q) return reverse ? ts.q > os.q : ts.q < os.q;

            reverse ^= 1u;

            if ((ts.r == zero) || (os.r == zero)) break;

            ts.n = ts.d;
            ts.d = ts.r;
            ts.q = ts.n / ts.d;
            ts.r = ts.n % ts.d;
            os.n = os.d;
            os.d = os.r;
            os.q = os.n / os.d;
            os.r = os.n % os.d;
        }

        if (ts.r == os.r)
            return false;
        else {
#pragma warning(push)
#pragma warning(disable : 4800)

            return (ts.r != zero) != static_cast<bool>(reverse);

#pragma warning(pop)
        }
    }

    constexpr bool operator>(const rational& other) const noexcept { return other < *this; }

    constexpr bool operator==(const rational& other) const noexcept
    {
        return num == other.numerator() && den == other.denominator();
    }

    constexpr operator bool() const noexcept { return num; }

    constexpr bool operator!() const noexcept { return !num; }

    static constexpr __int_type zero{0}, one{1};

private:
    __int_type num{0}, den{1};

    constexpr void normalize()
    {
        if (den == zero) throw bad_rational();
        if (num == zero) {
            den = one;
            return;
        }

        const auto gcd  = std::gcd(num, den);
        num            /= gcd;
        den            /= gcd;

        if (den < -(std::numeric_limits<__int_type>::max())) throw bad_rational("bad rational: non-zero singular denom");

        if (den < zero) {
            num = -num;
            den = -den;
        }

        assert(den > zero && std::gcd(num, den) == one);
    }

    static constexpr bool is_normalized(__int_type const& n, __int_type const& d)
    {
        return d > zero && (n != zero || d == one) && detail::abs(std::gcd(n, d)) == one;
    }

    template <typename T>
    static constexpr bool is_safe_narrowing_conversion(T const& x)
    {
        if constexpr (std::numeric_limits<T>::digits > std::numeric_limits<__int_type>::digits) {
            if constexpr (std::numeric_limits<T>::is_signed == false) {
                return x < (T{1} << std::numeric_limits<T>::digits);
            } else {
                if constexpr (std::numeric_limits<__int_type>::is_signed == false)
                    return x < (T{1} << std::numeric_limits<T>::digits) && x >= 0;
                else
                    return x < (T{1} << (std::numeric_limits<T>::digits - 1))
                           || x >= -(T{1} << (std::numeric_limits<T>::digits - 1));
            }
        } else {
            if constexpr (std::numeric_limits<T>::is_signed == true && std::numeric_limits<__int_type>::is_signed == false)
                return x >= 0;
            else
                return true;
        }
    }
};

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator+(const rational<__int_type>& r, const T& n)
{
    auto tmp  = r;
    tmp      += n;
    return tmp;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator-(const rational<__int_type>& r, const T& n)
{
    auto tmp  = r;
    tmp      -= n;
    return tmp;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator-(const T& n, const rational<__int_type>& r)
{
    auto tmp  = r;
    tmp      -= n;
    return -tmp;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator*(const rational<__int_type>& r, const T& n)
{
    auto tmp  = r;
    tmp      *= n;
    return tmp;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator/(const rational<__int_type>& r, const T& n)
{
    auto tmp  = r;
    tmp      /= n;
    return tmp;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator/(const T& n, const rational<__int_type>& r)
{
    rational<__int_type> tmp(n);
    tmp /= r;
    return tmp;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator<(const rational<__int_type>& r, const T& n)
{
    assert(r.denominator() > r.zero);
    __int_type _q = r.numerator() / r.denominator(), _r = r.numerator() % r.denominator();
    while (_r < r.zero) {
        _r += r.denominator();
        --_q;
    }

    return _q < n;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator<(const T& n, const rational<__int_type>& r)
{
    return r.operator>(n);
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator>(const rational<__int_type>& r, const T& n)
{
    return r == n ? false : !(r < n);
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator>(const T& n, const rational<__int_type>& r)
{
    return r.operator<(n);
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator<=(const rational<__int_type>& r, const T& n)
{
    return !r.operator>(n);
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator<=(const T& n, const rational<__int_type>& r)
{
    return r >= n;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator>=(const rational<__int_type>& r, const T& n)
{
    return !r.operator<(n);
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator>=(const T& n, const rational<__int_type>& r)
{
    return r <= n;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator==(const rational<__int_type>& r, const T& n)
{
    return r.denominator() == r.one && r.numerator() == n;
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator!=(const rational<__int_type>& r, const T& n)
{
    return !r.operator==(n);
}

template <typename T, typename __int_type, typename = std::enable_if_t<detail::is_compatible_integer_v<T, __int_type>>>
inline constexpr rational<__int_type> operator!=(const T& n, const rational<__int_type>& r)
{
    return !(n == r);
}

template <typename T, typename __int_type>
inline constexpr T rational_cast(const rational<__int_type>& r)
{
    return static_cast<T>(r.numerator()) / static_cast<T>(r.denominator());
}

namespace std
{
template <typename __int_type>
inline constexpr rational<__int_type> abs(const rational<__int_type>& r)
{
    return r.numerator() >= __int_type{0} ? r : -r;
}
} // namespace std

using i64_rational = rational<int64_t>;