double_integer.tcc

#ifndef DOUBLE_INTEGER_TCC
#define DOUBLE_INTEGER_TCC

/* Double Integer template
// double_integer.tcc
//
// Author: Jan Ringos, http://Tringi.TrimCore.cz, tringi@trimcore.cz
// Version: 2.3
//
// License:
//   This software is provided 'as-is', without any express or implied warranty.
//   In no event will the author be held liable for any damages arising from the
//   use of this software. 
// 
//   Permission is granted to anyone to use this software for any purpose,
//   including commercial applications, and to alter it and redistribute it
//   freely, subject to the following restrictions: 
//
//       1. The origin of this software must not be misrepresented; you must not
//          claim that you wrote the original software. If you use this software
//          in a product, an acknowledgment in the product documentation would
//          be appreciated but is not required.
//
//       2. Altered source versions must be plainly marked as such, and must not
//          be misrepresented as being the original software.
//
//       3. This notice may not be removed or altered from any source
//          distribution.
//
// Changelog:
//      29.01.2011 - initial version
//      28.05.2013 - fixed conversion from negative floating point
//                 - C style fmod, fmodf, fmodl replaced by std::fmod
//      20.03.2015 - completed C++11 release
//      29.09.2023 - removed deprecated numeric_limits denorm_absent/has_denorm_loss
//                 - silenced several comparison warnings
//      30.05.2024 - added explicit conversion to LO and double, tweaked print
//      01.06.2024 - fixed numeric_limits::min, and modernized the numeric_limits
*/

#include "double_integer.hpp"

#include <climits>
#include <cstring>
#include <cctype>
#include <cmath>

#include <limits>
#include <stdexcept>

template <typename LO, typename HI>
constexpr double_integer <LO, HI> ::operator LO () const {
    return this->lo;
};

template <typename LO, typename HI>
constexpr double_integer <LO, HI> ::operator bool () const {
    return !!this->hi || !!this->lo;
};

template <typename LO, typename HI>
bool constexpr
double_integer <LO, HI> ::operator ! () const {
    return !bool(*this);
};

template <typename LO, typename HI>
double_integer <LO, HI> 
double_integer <LO, HI> ::operator - () const {
    double_integer <LO, HI> result = ~*this;
    ++result;
    return result;
};

template <typename LO, typename HI>
double_integer <LO, HI>
double_integer <LO, HI> ::operator ~ () const {
    return double_integer <LO, HI> (~this->lo, ~this->hi);
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator ++ () {
    ++this->lo;
    if (!this->lo) {
        ++this->hi;
    };
    return *this;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator -- () {
    if (!this->lo) {
        --this->hi;
    };
    --this->lo;
    return *this;
};

template <typename LO, typename HI>
double_integer <LO, HI>
double_integer <LO, HI> ::operator ++ (int) {
    double_integer <LO, HI> previous = *this;
    ++*this;
    return previous;
};

template <typename LO, typename HI>
double_integer <LO, HI>
double_integer <LO, HI> ::operator -- (int) {
    double_integer <LO, HI> previous = *this;
    --*this;
    return previous;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator += (const double_integer <LO, HI> & addend) {
    LO previous_lo = this->lo;

    this->lo += addend.lo;
    this->hi += addend.hi;
    this->hi += (this->lo < previous_lo);
    
    return *this;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator -= (const double_integer <LO, HI> & subtrahend) {
    return *this += -subtrahend;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator *= (const double_integer <LO, HI> & _multiplier) {
    if (!_multiplier) {
        *this = double_integer <LO, HI> ();
        return *this;
    };

    double_integer <LO, HI> previous = *this;
    double_integer <LO, HI> multiplier = _multiplier;
    
    this->lo = LO ();
    this->hi = HI ();
    
    auto i = 0u;
    auto n = CHAR_BIT * (sizeof (LO) + sizeof (HI));
    
    while (i < n && !!multiplier) {
        if (!!(multiplier.lo & LO (1u))) {
            *this += previous << i;
        };
        
        multiplier >>= 1u;
        ++i;
    };
    
    return *this;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator /= (const double_integer <LO, HI> & divisor) {
    *this = this->div (divisor);
    return *this;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator %= (const double_integer <LO, HI> & divisor) {
    this->div (divisor, this);
    return *this;
};

template <typename LO, typename HI>
double_integer <LO, HI>
double_integer <LO, HI> ::div (const double_integer <LO, HI> & _ds,
                               double_integer <LO, HI> * remainder) const {
    
//    if (!_ds)
//        return this->lo / _ds.lo;
    
    if (!_ds)
        throw std::logic_error ("div/0");
    
    // reset sign bits
    
    double_integer <LO, HI> zero;
    double_integer <LO, HI> ds = _ds;
    double_integer <LO, HI> dd = *this;
    
    if (ds < zero) {
        ds = -ds;
    };
    if (dd < zero) {
        dd = -dd;
    };

    // remainder only
    //  - shortcut, return zero
    
    if (dd < ds) {
        if (remainder) {
            *remainder = *this;
        };
        return zero;
    };
    
    double_integer <LO, HI> r;
    double_integer <LO, HI> q;
    double_integer <LO, HI> one;
    ++one;
    
    unsigned int b = CHAR_BIT * (sizeof (LO) + sizeof (HI));
    
    while (r < ds) {
        r <<= 1u;
        if (!!(dd & (one << --b))) {
            r.lo |= one.lo;
        };
    };
    
    while (true) {
        if (r < ds) {
            if (!b) break;
            
            --b;
            r <<= 1u;

            if (!!(dd & (one << b))) {
                r.lo |= one.lo;
            };
        } else {
            r -= ds;
            q |= one << b;
        };
    };

    // correct sign bits
    
    if ((_ds < zero) != (*this < zero)) {
        q = -q;
    };
    
    // results

    if (remainder) {
        if (*this < zero) {
            r = -r;
        };
        *remainder = r;
    };
    return q;
    
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator <<= (unsigned int bits) {

    // discard upper bits of the shift count, just like shift for ints does
    bits &= (CHAR_BIT * (sizeof (LO) + sizeof (HI))) - 1u;
    
    if (bits >= CHAR_BIT * sizeof (LO)) {
        bits -= CHAR_BIT * sizeof (LO);
        
        this->hi = this->lo;
        this->lo = LO ();
        this->hi <<= bits;
        
    } else {
        const unsigned int rest = CHAR_BIT * sizeof (LO) - bits;
        
        this->hi = (LO (this->hi) << bits)
                 | ((this->lo & (this->mask <LO> (bits) << rest)) >> rest);
        this->lo <<= bits;
    };
    
    return *this;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator >>= (unsigned int bits) {
    
    // discard upper bits of the shift count, just like shift for ints does
    bits &= (CHAR_BIT * (sizeof (LO) + sizeof (HI))) - 1u;
    
    if (bits >= CHAR_BIT * sizeof (LO)) {
        bits -= CHAR_BIT * sizeof (LO);
        
        this->lo = this->hi;
        this->lo >>= bits;
        
        if (std::numeric_limits <HI> ::is_signed
                && double_integer (LO (), this->hi) < double_integer ()) {
            
            this->hi = ~HI ();
            this->lo |= this->mask <LO> (bits)
                    << (CHAR_BIT * sizeof (LO) - bits);
        } else {
            this->hi = HI ();
        };

    } else {
        this->lo >>= bits;
        this->lo |= (LO (this->hi) & this->mask <LO> (bits))
                 << (CHAR_BIT * sizeof (LO) - bits);
        this->hi >>= bits;
    };
    
    return *this;
};

template <typename LO, typename HI>
template <typename T>
T double_integer <LO, HI> ::mask (unsigned int n) {
    T m = T ();
    
    ++m;
    m <<= n;
    --m;
    
    return m;
};

namespace {
    template <typename LO, typename HI, typename T>
    constexpr HI double_integer_init_HI (const T & a) {
        if constexpr (sizeof (T) > sizeof (LO))
            return a >> (sizeof (LO) * CHAR_BIT);
        else
            return HI ();
    }
    template <typename LO, typename HI, typename T>
    constexpr HI double_integer_init_HI_signed (const T & a) {
        if constexpr (sizeof (T) > sizeof (LO))
            return a >> (sizeof (LO) * CHAR_BIT);
        else
        if (a < 0)
            return HI (-1);
        else
            return HI (0);
    }
}

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (unsigned char a)
    :   lo (LO (a)),
        hi (HI ()) {};
        
template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (unsigned short a)
    :   lo (LO (a)),
        hi (double_integer_init_HI <LO, HI> (a)) {};
        
template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (unsigned int a)
    :   lo (LO (a)),
        hi (double_integer_init_HI <LO, HI> (a)) {};
        
template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (unsigned long a)
    :   lo (LO (a)),
        hi (double_integer_init_HI <LO, HI> (a)) {};
        
template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (unsigned long long a)
    :   lo (LO (a)),
        hi (double_integer_init_HI <LO, HI> (a)) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (signed char a)
    :   lo (LO (a)),
        hi (HI ((a < 0) ? -1 : 0)) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (signed short a)
    :   lo (LO (a)),
        hi (double_integer_init_HI_signed <LO, HI> (a)) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (signed int a)
    :   lo (LO (a)),
        hi (double_integer_init_HI_signed <LO, HI> (a)) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (signed long a)
    :   lo (LO (a)),
        hi (double_integer_init_HI_signed <LO, HI> (a)) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (signed long long a)
    :   lo (LO (a)),
        hi (double_integer_init_HI_signed <LO, HI> (a)) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (float a)
    :   lo (std::fmod (a,      static_cast <float> (LO (~LO ())) + 1.0f)),
        hi ((a < 0.0f) ? (a / (static_cast <float> (LO (~LO ())) + 1.0f) - 1u)
                       : (a / (static_cast <float> (LO (~LO ())) + 1.0f))) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (double a)
    :   lo (std::fmod (a,     static_cast <double> (LO (~LO ())) + 1.0)),
        hi ((a < 0.0) ? (a / (static_cast <double> (LO (~LO ())) + 1.0) - 1u)
                      : (a / (static_cast <double> (LO (~LO ())) + 1.0))) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (long double a)
    :   lo (std::fmod (a,      static_cast <long double> (LO (~LO ())) + 1.0l)),
        hi ((a < 0.0l) ? (a / (static_cast <long double> (LO (~LO ())) + 1.0l) - 1u)
                       : (a / (static_cast <long double> (LO (~LO ())) + 1.0l))) {};

template <typename LO, typename HI>
double_integer <LO, HI> ::double_integer (const char * s, unsigned int _radix)
    :   lo (LO ()),
        hi (HI ()) {
    
    if (!s)
        return;
    
    double_integer <LO, HI> digit;
    double_integer <LO, HI> radix;
    bool minus = false;
    
    // skip blanks to first character
    while (std::isspace (*s)) ++s;

    if (*s == '-') {
        ++s;
        minus = true;
        
        // skip optional spaces after minus
        while (std::isspace (*s)) ++s;
    };
    
    // if no radix was provided, autodetect
    if (_radix) {
        radix.lo = _radix;
    } else {
        if (*s == '0') {
            ++s;
            radix.lo = 8u;
            
            if (*s == 'x' || *s == 'X') {
                ++s;
                radix.lo = 16u;
            };
            if (*s == 'b' || *s == 'B') {
                ++s;
                radix.lo = 2u;
            };
        } else {
            radix.lo = 10u;
        };
    };
    
    while (std::isxdigit (*s)) {
        if (std::isdigit (*s)) {
            digit.lo = static_cast <unsigned int> (*s - '0');
        } else
        if (std::islower (*s)) {
            digit.lo = static_cast <unsigned int> (*s - 'a' + 10u);
        } else {
            digit.lo = static_cast <unsigned int> (*s - 'A' + 10u);
        };
        
        *this *= radix;
        *this += digit;
        ++s;
    };
    
    if (minus) {
        *this = -*this;
    };
    return;
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator = (const double_integer <LO, HI> & x) {
    this->lo = x.lo;
    this->hi = x.hi;
    return *this;
};

template <typename LO, typename HI>
template <class>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator = (const double_integer <LO, LO> & x) {
    this->lo = x.lo;
    this->hi = x.hi;
    return *this;
};

template <typename LO, typename HI>
constexpr
double_integer <LO, HI> ::operator double_integer <LO, LO> () const {
    return double_integer <LO, LO> (this->lo, this->hi);
};

template <typename LO, typename HI>
constexpr double_integer <LO, HI> ::operator double () const {
    return ((double) this->lo)
         + ((double) this->hi * ((double) std::numeric_limits <LO> ::max () + 1.0));
}

template <typename LO, typename HI>
const char * double_integer <LO, HI> ::print (char * buffer, std::size_t length,
                                              unsigned int radix) const {
    if (!*this) {
        if (length >= 2) {
            buffer [0] = '0';
            buffer [1] = 0;
            return buffer;
        } else
            return "0";
    }

    if (radix < 2 || radix > 37)
        return "(inv radix)";

    if (length) {
        std::size_t i = length;
        buffer [--i] = '\0';

        double_integer <LO, HI> zero;
        double_integer <LO, HI> r;
        double_integer <LO, HI> x = *this;
        
        if (x < zero) {
            x = -x;
            if (x < zero) {
                
                buffer [0] = '-';
                double_integer <LO, LO>
                    (this->lo, this->hi) .print (buffer + 1, length - 1, radix);
                
                return buffer;
            };
        };
        
        while (i && !!x) {
            x = x.div (double_integer <LO, HI> (radix), &r);
            
            unsigned char c = *reinterpret_cast <unsigned char *> (&r);
            if (c < 10u)
                buffer [--i] = c + '0';
            else
                buffer [--i] = c + 'A' - 10u;
        };

        if (*this < zero) {
            if (i) {
                buffer [--i] = '-';
            } else
                return NULL;
        };
        
        if (!x) {
            std::memmove (buffer, buffer + i, length - i);
            return buffer;
        };
    };
    
    return NULL;
};

template <typename LO, typename HI>
const char * double_integer <LO, HI> ::c_str (unsigned int radix) const {
    static char sz [CHAR_BIT * (sizeof (LO) + sizeof (HI)) + 2u];
    return this->print (sz, sizeof sz, radix);
};

template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator |= (const double_integer <LO, HI> & param) {
    this->lo |= param.lo;
    this->hi |= param.hi;
    return *this;
};
template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator &= (const double_integer <LO, HI> & param) {
    this->lo &= param.lo;
    this->hi &= param.hi;
    return *this;
};
template <typename LO, typename HI>
double_integer <LO, HI> &
double_integer <LO, HI> ::operator ^= (const double_integer <LO, HI> & param) {
    this->lo ^= param.lo;
    this->hi ^= param.hi;
    return *this;
};

template <typename LO, typename HI>
bool operator < (const double_integer <LO, HI> & a,
                 const double_integer <LO, HI> & b) {
    if (a.hi == b.hi)
        return a.lo < b.lo;
    else
        return a.hi < b.hi;
};
template <typename LO, typename HI>
bool operator == (const double_integer <LO, HI> & a,
                  const double_integer <LO, HI> & b) {
    return a.lo == b.lo
        && a.hi == b.hi;
};
template <typename LO, typename HI>
bool operator != (const double_integer <LO, HI> & a,
                  const double_integer <LO, HI> & b) {
    return a.lo != b.lo
        || a.hi != b.hi;
};
template <typename LO, typename HI>
bool operator || (const double_integer <LO, HI> & a,
                  const double_integer <LO, HI> & b) {
    return (a.hi || a.lo)
        || (b.hi || b.lo);
};
template <typename LO, typename HI>
bool operator && (const double_integer <LO, HI> & a,
                  const double_integer <LO, HI> & b) {
    return (a.hi || a.lo)
        && (b.hi || b.lo);
};

// non-shortcut arithmetic operators implemented as global functions

template <typename LO, typename HI>
double_integer <LO, HI> operator + (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) += b;
};
template <typename LO, typename HI>
double_integer <LO, HI> operator - (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) -= b;
};
template <typename LO, typename HI>
double_integer <LO, HI> operator * (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) *= b;
};
template <typename LO, typename HI>
double_integer <LO, HI> operator / (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) /= b;
};
template <typename LO, typename HI>
double_integer <LO, HI> operator % (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) %= b;
};

template <typename LO, typename HI>
double_integer <LO, HI> operator << (const double_integer <LO, HI> & a,
                                     unsigned int n) {
    return double_integer <LO, HI> (a) <<= n;
};
template <typename LO, typename HI>
double_integer <LO, HI> operator >> (const double_integer <LO, HI> & a,
                                     unsigned int n) {
    return double_integer <LO, HI> (a) >>= n;
};

template <typename LO, typename HI>
double_integer <LO, HI> operator & (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) &= b;
};
template <typename LO, typename HI>
double_integer <LO, HI> operator | (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) |= b;
};
template <typename LO, typename HI>
double_integer <LO, HI> operator ^ (const double_integer <LO, HI> & a,
                                    const double_integer <LO, HI> & b) {
    return double_integer <LO, HI> (a) ^= b;
};

template <typename LO, typename HI>
bool operator >  (const double_integer <LO, HI> & a,
                  const double_integer <LO, HI> & b) { return  (b < a); };
template <typename LO, typename HI>
bool operator <= (const double_integer <LO, HI> & a,
                  const double_integer <LO, HI> & b) { return !(b < a); };
template <typename LO, typename HI>
bool operator >= (const double_integer <LO, HI> & a,
                  const double_integer <LO, HI> & b) { return !(a < b); };

// std::numeric_limits
//  - and some type_traits overrides

namespace std {

    template <typename LO, typename HI>
    struct numeric_limits <double_integer <LO, HI> > {
        static constexpr bool is_specialized = true;
        
        [[nodiscard]] static constexpr double_integer <LO, HI> min () noexcept {
            return double_integer <LO, HI> (is_signed ? std::numeric_limits <LO> ::max ()
                                                      : std::numeric_limits <LO> ::min (),
                                            std::numeric_limits <HI> ::min ());
        };
        
        [[nodiscard]] static constexpr double_integer <LO, HI> max () noexcept {
            return double_integer <LO, HI> (std::numeric_limits <LO> ::max (),
                                            std::numeric_limits <HI> ::max ());
        };

        [[nodiscard]] static constexpr double_integer <LO, HI> lowest () noexcept {
            return min ();
        };

        static constexpr int digits = std::numeric_limits <HI> ::digits
                                    + std::numeric_limits <LO> ::digits;
        static constexpr int digits10 = digits * 643 / 2136;
        
        static constexpr bool is_signed = std::numeric_limits <HI> ::is_signed;
        static constexpr bool is_integer = true;
        static constexpr bool is_exact = true;
        static constexpr int radix = 2;

        [[nodiscard]] static constexpr double_integer <LO, HI> epsilon () noexcept { return {}; }
        [[nodiscard]] static constexpr double_integer <LO, HI> round_error () noexcept { return {}; }
        
        static constexpr int min_exponent = 0;
        static constexpr int min_exponent10 = 0;
        static constexpr int max_exponent = 0;
        static constexpr int max_exponent10 = 0;
        static constexpr int max_digis10 = 0;
        
        static constexpr bool has_infinity = false;
        static constexpr bool has_quiet_NaN = false;
        static constexpr bool has_signaling_NaN = false;
        
        [[nodiscard]] static constexpr double_integer <LO, HI> infinity () noexcept { return {}; }
        [[nodiscard]] static constexpr double_integer <LO, HI> quiet_NaN () noexcept { return {}; }
        [[nodiscard]] static constexpr double_integer <LO, HI> signaling_NaN () noexcept { return {}; }
        [[nodiscard]] static constexpr double_integer <LO, HI> denorm_min () noexcept { return {}; }
        
        static constexpr bool is_iec559 = false;
        static constexpr bool is_bounded = true;
        static constexpr bool is_modulo = true;
        
        static constexpr bool traps = std::numeric_limits <HI> ::traps;
        static constexpr bool tinyness_before = false;
        static constexpr float_round_style round_style = round_indeterminate;
    };
    
    template <typename LO, typename HI>
    struct make_unsigned <double_integer <LO, HI> > {
        typedef double_integer <LO, LO> type;
    };

    template <typename LO, typename HI>
    struct is_integral <double_integer <LO, HI> >
        : std::integral_constant <bool, true> {};
};

#endif