Number.cpp

#include <xrpl/basics/Number.h>
#include <xrpl/beast/utility/instrumentation.h>
 
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
 
#ifdef _MSC_VER
#pragma message("Using boost::multiprecision::uint128_t")
#include <boost/multiprecision/cpp_int.hpp>
using uint128_t = boost::multiprecision::uint128_t;
#else   // !defined(_MSC_VER)
using uint128_t = __uint128_t;
#endif  // !defined(_MSC_VER)
 
namespace xrpl {
 
thread_local Number::rounding_mode Number::mode_ = Number::to_nearest;
 
Number::rounding_mode
Number::getround()
{
    return mode_;
}
 
Number::rounding_mode
Number::setround(rounding_mode mode)
{
    return std::exchange(mode_, mode);
}
 
// Guard
 
// The Guard class is used to temporarily add extra digits of
// precision to an operation.  This enables the final result
// to be correctly rounded to the internal precision of Number.
 
class Number::Guard
{
    std::uint64_t digits_;   // 16 decimal guard digits
    std::uint8_t xbit_ : 1;  // has a non-zero digit been shifted off the end
    std::uint8_t sbit_ : 1;  // the sign of the guard digits
 
public:
    explicit Guard() : digits_{0}, xbit_{0}, sbit_{0}
    {
    }
 
    // set & test the sign bit
    void
    set_positive() noexcept;
    void
    set_negative() noexcept;
    bool
    is_negative() const noexcept;
 
    // add a digit
    void
    push(unsigned d) noexcept;
 
    // recover a digit
    unsigned
    pop() noexcept;
 
    // Indicate round direction:  1 is up, -1 is down, 0 is even
    // This enables the client to round towards nearest, and on
    // tie, round towards even.
    int
    round() noexcept;
 
    // Modify the result to the correctly rounded value
    void
    doRoundUp(rep& mantissa, int& exponent, std::string location);
 
    // Modify the result to the correctly rounded value
    void
    doRoundDown(rep& mantissa, int& exponent);
 
    // Modify the result to the correctly rounded value
    void
    doRound(rep& drops);
};
 
inline void
Number::Guard::set_positive() noexcept
{
    sbit_ = 0;
}
 
inline void
Number::Guard::set_negative() noexcept
{
    sbit_ = 1;
}
 
inline bool
Number::Guard::is_negative() const noexcept
{
    return sbit_ == 1;
}
 
inline void
Number::Guard::push(unsigned d) noexcept
{
    xbit_ = xbit_ || (digits_ & 0x0000'0000'0000'000F) != 0;
    digits_ >>= 4;
    digits_ |= (d & 0x0000'0000'0000'000FULL) << 60;
}
 
inline unsigned
Number::Guard::pop() noexcept
{
    unsigned d = (digits_ & 0xF000'0000'0000'0000) >> 60;
    digits_ <<= 4;
    return d;
}
 
// Returns:
//     -1 if Guard is less than half
//      0 if Guard is exactly half
//      1 if Guard is greater than half
int
Number::Guard::round() noexcept
{
    auto mode = Number::getround();
 
    if (mode == towards_zero)
        return -1;
 
    if (mode == downward)
    {
        if (sbit_)
        {
            if (digits_ > 0 || xbit_)
                return 1;
        }
        return -1;
    }
 
    if (mode == upward)
    {
        if (sbit_)
            return -1;
        if (digits_ > 0 || xbit_)
            return 1;
        return -1;
    }
 
    // assume round to nearest if mode is not one of the predefined values
    if (digits_ > 0x5000'0000'0000'0000)
        return 1;
    if (digits_ < 0x5000'0000'0000'0000)
        return -1;
    if (xbit_)
        return 1;
    return 0;
}
 
void
Number::Guard::doRoundUp(rep& mantissa, int& exponent, std::string location)
{
    auto r = round();
    if (r == 1 || (r == 0 && (mantissa & 1) == 1))
    {
        ++mantissa;
        if (mantissa > maxMantissa)
        {
            mantissa /= 10;
            ++exponent;
        }
    }
    if (exponent < minExponent)
    {
        mantissa = 0;
        exponent = Number{}.exponent_;
    }
    if (exponent > maxExponent)
        throw std::overflow_error(location);
}
 
void
Number::Guard::doRoundDown(rep& mantissa, int& exponent)
{
    auto r = round();
    if (r == 1 || (r == 0 && (mantissa & 1) == 1))
    {
        --mantissa;
        if (mantissa < minMantissa)
        {
            mantissa *= 10;
            --exponent;
        }
    }
    if (exponent < minExponent)
    {
        mantissa = 0;
        exponent = Number{}.exponent_;
    }
}
 
// Modify the result to the correctly rounded value
void
Number::Guard::doRound(rep& drops)
{
    auto r = round();
    if (r == 1 || (r == 0 && (drops & 1) == 1))
    {
        ++drops;
    }
    if (is_negative())
        drops = -drops;
}
 
// Number
 
constexpr Number one{1000000000000000, -15, Number::unchecked{}};
 
void
Number::normalize()
{
    if (mantissa_ == 0)
    {
        *this = Number{};
        return;
    }
    bool const negative = (mantissa_ < 0);
    auto m = static_cast<std::make_unsigned_t<rep>>(mantissa_);
    if (negative)
        m = -m;
    while ((m < minMantissa) && (exponent_ > minExponent))
    {
        m *= 10;
        --exponent_;
    }
    Guard g;
    if (negative)
        g.set_negative();
    while (m > maxMantissa)
    {
        if (exponent_ >= maxExponent)
            throw std::overflow_error("Number::normalize 1");
        g.push(m % 10);
        m /= 10;
        ++exponent_;
    }
    mantissa_ = m;
    if ((exponent_ < minExponent) || (mantissa_ < minMantissa))
    {
        *this = Number{};
        return;
    }
 
    g.doRoundUp(mantissa_, exponent_, "Number::normalize 2");
 
    if (negative)
        mantissa_ = -mantissa_;
}
 
Number&
Number::operator+=(Number const& y)
{
    if (y == Number{})
        return *this;
    if (*this == Number{})
    {
        *this = y;
        return *this;
    }
    if (*this == -y)
    {
        *this = Number{};
        return *this;
    }
    XRPL_ASSERT(
        isnormal() && y.isnormal(),
        "xrpl::Number::operator+=(Number) : is normal");
    auto xm = mantissa();
    auto xe = exponent();
    int xn = 1;
    if (xm < 0)
    {
        xm = -xm;
        xn = -1;
    }
    auto ym = y.mantissa();
    auto ye = y.exponent();
    int yn = 1;
    if (ym < 0)
    {
        ym = -ym;
        yn = -1;
    }
    Guard g;
    if (xe < ye)
    {
        if (xn == -1)
            g.set_negative();
        do
        {
            g.push(xm % 10);
            xm /= 10;
            ++xe;
        } while (xe < ye);
    }
    else if (xe > ye)
    {
        if (yn == -1)
            g.set_negative();
        do
        {
            g.push(ym % 10);
            ym /= 10;
            ++ye;
        } while (xe > ye);
    }
    if (xn == yn)
    {
        xm += ym;
        if (xm > maxMantissa)
        {
            g.push(xm % 10);
            xm /= 10;
            ++xe;
        }
        g.doRoundUp(xm, xe, "Number::addition overflow");
    }
    else
    {
        if (xm > ym)
        {
            xm = xm - ym;
        }
        else
        {
            xm = ym - xm;
            xe = ye;
            xn = yn;
        }
        while (xm < minMantissa)
        {
            xm *= 10;
            xm -= g.pop();
            --xe;
        }
        g.doRoundDown(xm, xe);
    }
    mantissa_ = xm * xn;
    exponent_ = xe;
    return *this;
}
 
// Optimization equivalent to:
// auto r = static_cast<unsigned>(u % 10);
// u /= 10;
// return r;
// Derived from Hacker's Delight Second Edition Chapter 10
// by Henry S. Warren, Jr.
static inline unsigned
divu10(uint128_t& u)
{
    // q = u * 0.75
    auto q = (u >> 1) + (u >> 2);
    // iterate towards q = u * 0.8
    q += q >> 4;
    q += q >> 8;
    q += q >> 16;
    q += q >> 32;
    q += q >> 64;
    // q /= 8 approximately == u / 10
    q >>= 3;
    // r = u - q * 10  approximately == u % 10
    auto r = static_cast<unsigned>(u - ((q << 3) + (q << 1)));
    // correction c is 1 if r >= 10 else 0
    auto c = (r + 6) >> 4;
    u = q + c;
    r -= c * 10;
    return r;
}
 
Number&
Number::operator*=(Number const& y)
{
    if (*this == Number{})
        return *this;
    if (y == Number{})
    {
        *this = y;
        return *this;
    }
    XRPL_ASSERT(
        isnormal() && y.isnormal(),
        "xrpl::Number::operator*=(Number) : is normal");
    auto xm = mantissa();
    auto xe = exponent();
    int xn = 1;
    if (xm < 0)
    {
        xm = -xm;
        xn = -1;
    }
    auto ym = y.mantissa();
    auto ye = y.exponent();
    int yn = 1;
    if (ym < 0)
    {
        ym = -ym;
        yn = -1;
    }
    auto zm = uint128_t(xm) * uint128_t(ym);
    auto ze = xe + ye;
    auto zn = xn * yn;
    Guard g;
    if (zn == -1)
        g.set_negative();
    while (zm > maxMantissa)
    {
        // The following is optimization for:
        // g.push(static_cast<unsigned>(zm % 10));
        // zm /= 10;
        g.push(divu10(zm));
        ++ze;
    }
    xm = static_cast<rep>(zm);
    xe = ze;
    g.doRoundUp(
        xm,
        xe,
        "Number::multiplication overflow : exponent is " + std::to_string(xe));
    mantissa_ = xm * zn;
    exponent_ = xe;
    XRPL_ASSERT(
        isnormal() || *this == Number{},
        "xrpl::Number::operator*=(Number) : result is normal");
    return *this;
}
 
Number&
Number::operator/=(Number const& y)
{
    if (y == Number{})
        throw std::overflow_error("Number: divide by 0");
    if (*this == Number{})
        return *this;
    int np = 1;
    auto nm = mantissa();
    auto ne = exponent();
    if (nm < 0)
    {
        nm = -nm;
        np = -1;
    }
    int dp = 1;
    auto dm = y.mantissa();
    auto de = y.exponent();
    if (dm < 0)
    {
        dm = -dm;
        dp = -1;
    }
    // Shift by 10^17 gives greatest precision while not overflowing uint128_t
    // or the cast back to int64_t
    uint128_t const f = 100'000'000'000'000'000;
    mantissa_ = static_cast<std::int64_t>(uint128_t(nm) * f / uint128_t(dm));
    exponent_ = ne - de - 17;
    mantissa_ *= np * dp;
    normalize();
    return *this;
}
 
Number::operator rep() const
{
    rep drops = mantissa_;
    int offset = exponent_;
    Guard g;
    if (drops != 0)
    {
        if (drops < 0)
        {
            g.set_negative();
            drops = -drops;
        }
        for (; offset < 0; ++offset)
        {
            g.push(drops % 10);
            drops /= 10;
        }
        for (; offset > 0; --offset)
        {
            if (drops > std::numeric_limits<decltype(drops)>::max() / 10)
                throw std::overflow_error("Number::operator rep() overflow");
            drops *= 10;
        }
        g.doRound(drops);
    }
    return drops;
}
 
Number
Number::truncate() const noexcept
{
    if (exponent_ >= 0 || mantissa_ == 0)
        return *this;
 
    Number ret = *this;
    while (ret.exponent_ < 0 && ret.mantissa_ != 0)
    {
        ret.exponent_ += 1;
        ret.mantissa_ /= rep(10);
    }
    // We are guaranteed that normalize() will never throw an exception
    // because exponent is either negative or zero at this point.
    ret.normalize();
    return ret;
}
 
std::string
to_string(Number const& amount)
{
    // keep full internal accuracy, but make more human friendly if possible
    if (amount == Number{})
        return "0";
 
    auto const exponent = amount.exponent();
    auto mantissa = amount.mantissa();
 
    // Use scientific notation for exponents that are too small or too large
    if (((exponent != 0) && ((exponent < -25) || (exponent > -5))))
    {
        std::string ret = std::to_string(mantissa);
        ret.append(1, 'e');
        ret.append(std::to_string(exponent));
        return ret;
    }
 
    bool negative = false;
 
    if (mantissa < 0)
    {
        mantissa = -mantissa;
        negative = true;
    }
 
    XRPL_ASSERT(
        exponent + 43 > 0, "xrpl::to_string(Number) : minimum exponent");
 
    ptrdiff_t const pad_prefix = 27;
    ptrdiff_t const pad_suffix = 23;
 
    std::string const raw_value(std::to_string(mantissa));
    std::string val;
 
    val.reserve(raw_value.length() + pad_prefix + pad_suffix);
    val.append(pad_prefix, '0');
    val.append(raw_value);
    val.append(pad_suffix, '0');
 
    ptrdiff_t const offset(exponent + 43);
 
    auto pre_from(val.begin());
    auto const pre_to(val.begin() + offset);
 
    auto const post_from(val.begin() + offset);
    auto post_to(val.end());
 
    // Crop leading zeroes. Take advantage of the fact that there's always a
    // fixed amount of leading zeroes and skip them.
    if (std::distance(pre_from, pre_to) > pad_prefix)
        pre_from += pad_prefix;
 
    XRPL_ASSERT(
        post_to >= post_from, "xrpl::to_string(Number) : first distance check");
 
    pre_from = std::find_if(pre_from, pre_to, [](char c) { return c != '0'; });
 
    // Crop trailing zeroes. Take advantage of the fact that there's always a
    // fixed amount of trailing zeroes and skip them.
    if (std::distance(post_from, post_to) > pad_suffix)
        post_to -= pad_suffix;
 
    XRPL_ASSERT(
        post_to >= post_from,
        "xrpl::to_string(Number) : second distance check");
 
    post_to = std::find_if(
                  std::make_reverse_iterator(post_to),
                  std::make_reverse_iterator(post_from),
                  [](char c) { return c != '0'; })
                  .base();
 
    std::string ret;
 
    if (negative)
        ret.append(1, '-');
 
    // Assemble the output:
    if (pre_from == pre_to)
        ret.append(1, '0');
    else
        ret.append(pre_from, pre_to);
 
    if (post_to != post_from)
    {
        ret.append(1, '.');
        ret.append(post_from, post_to);
    }
 
    return ret;
}
 
// Returns f^n
// Uses a log_2(n) number of multiplications
 
Number
power(Number const& f, unsigned n)
{
    if (n == 0)
        return one;
    if (n == 1)
        return f;
    auto r = power(f, n / 2);
    r *= r;
    if (n % 2 != 0)
        r *= f;
    return r;
}
 
// Returns f^(1/d)
// Uses Newton–Raphson iterations until the result stops changing
// to find the non-negative root of the polynomial g(x) = x^d - f
 
// This function, and power(Number f, unsigned n, unsigned d)
// treat corner cases such as 0 roots as advised by Annex F of
// the C standard, which itself is consistent with the IEEE
// floating point standards.
 
Number
root(Number f, unsigned d)
{
    if (f == one || d == 1)
        return f;
    if (d == 0)
    {
        if (f == -one)
            return one;
        if (abs(f) < one)
            return Number{};
        throw std::overflow_error("Number::root infinity");
    }
    if (f < Number{} && d % 2 == 0)
        throw std::overflow_error("Number::root nan");
    if (f == Number{})
        return f;
 
    // Scale f into the range (0, 1) such that f's exponent is a multiple of d
    auto e = f.exponent() + 16;
    auto const di = static_cast<int>(d);
    auto ex = [e = e, di = di]()  // Euclidean remainder of e/d
    {
        int k = (e >= 0 ? e : e - (di - 1)) / di;
        int k2 = e - k * di;
        if (k2 == 0)
            return 0;
        return di - k2;
    }();
    e += ex;
    f = Number{f.mantissa(), f.exponent() - e};  // f /= 10^e;
    bool neg = false;
    if (f < Number{})
    {
        neg = true;
        f = -f;
    }
 
    // Quadratic least squares curve fit of f^(1/d) in the range [0, 1]
    auto const D = ((6 * di + 11) * di + 6) * di + 1;
    auto const a0 = 3 * di * ((2 * di - 3) * di + 1);
    auto const a1 = 24 * di * (2 * di - 1);
    auto const a2 = -30 * (di - 1) * di;
    Number r = ((Number{a2} * f + Number{a1}) * f + Number{a0}) / Number{D};
    if (neg)
    {
        f = -f;
        r = -r;
    }
 
    //  Newton–Raphson iteration of f^(1/d) with initial guess r
    //  halt when r stops changing, checking for bouncing on the last iteration
    Number rm1{};
    Number rm2{};
    do
    {
        rm2 = rm1;
        rm1 = r;
        r = (Number(d - 1) * r + f / power(r, d - 1)) / Number(d);
    } while (r != rm1 && r != rm2);
 
    //  return r * 10^(e/d) to reverse scaling
    return Number{r.mantissa(), r.exponent() + e / di};
}
 
Number
root2(Number f)
{
    if (f == one)
        return f;
    if (f < Number{})
        throw std::overflow_error("Number::root nan");
    if (f == Number{})
        return f;
 
    // Scale f into the range (0, 1) such that f's exponent is a multiple of d
    auto e = f.exponent() + 16;
    if (e % 2 != 0)
        ++e;
    f = Number{f.mantissa(), f.exponent() - e};  // f /= 10^e;
 
    // Quadratic least squares curve fit of f^(1/d) in the range [0, 1]
    auto const D = 105;
    auto const a0 = 18;
    auto const a1 = 144;
    auto const a2 = -60;
    Number r = ((Number{a2} * f + Number{a1}) * f + Number{a0}) / Number{D};
 
    //  Newton–Raphson iteration of f^(1/2) with initial guess r
    //  halt when r stops changing, checking for bouncing on the last iteration
    Number rm1{};
    Number rm2{};
    do
    {
        rm2 = rm1;
        rm1 = r;
        r = (r + f / r) / Number(2);
    } while (r != rm1 && r != rm2);
 
    //  return r * 10^(e/2) to reverse scaling
    return Number{r.mantissa(), r.exponent() + e / 2};
}
 
// Returns f^(n/d)
 
Number
power(Number const& f, unsigned n, unsigned d)
{
    if (f == one)
        return f;
    auto g = std::gcd(n, d);
    if (g == 0)
        throw std::overflow_error("Number::power nan");
    if (d == 0)
    {
        if (f == -one)
            return one;
        if (abs(f) < one)
            return Number{};
        // abs(f) > one
        throw std::overflow_error("Number::power infinity");
    }
    if (n == 0)
        return one;
    n /= g;
    d /= g;
    if ((n % 2) == 1 && (d % 2) == 0 && f < Number{})
        throw std::overflow_error("Number::power nan");
    return root(power(f, n), d);
}
 
}  // namespace xrpl