Number.cpp

//------------------------------------------------------------------------------
/*
    This file is part of rippled: https://github.com/ripple/rippled
    Copyright (c) 2022 Ripple Labs Inc.
 
    Permission to use, copy, modify, and/or distribute this software for any
    purpose  with  or without fee is hereby granted, provided that the above
    copyright notice and this permission notice appear in all copies.
 
    THE  SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
    WITH  REGARD  TO  THIS  SOFTWARE  INCLUDING  ALL  IMPLIED  WARRANTIES  OF
    MERCHANTABILITY  AND  FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
    ANY  SPECIAL ,  DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
    WHATSOEVER  RESULTING  FROM  LOSS  OF USE, DATA OR PROFITS, WHETHER IN AN
    ACTION  OF  CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
    OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
//==============================================================================
 
#include <xrpl/basics/Number.h>
#include <xrpl/beast/utility/instrumentation.h>
#include <boost/predef.h>
#include <algorithm>
#include <numeric>
#include <stdexcept>
#include <type_traits>
#include <utility>
 
#ifdef BOOST_COMP_MSVC
#include <boost/multiprecision/cpp_int.hpp>
using uint128_t = boost::multiprecision::uint128_t;
#else   // !defined(_MSVC_LANG)
using uint128_t = __uint128_t;
#endif  // !defined(_MSVC_LANG)
 
namespace ripple {
 
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 tempoarily add extra digits of
// preicision 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;
};
 
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;
}
 
// 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;
    }
 
    auto r = g.round();
    if (r == 1 || (r == 0 && (mantissa_ & 1) == 1))
    {
        ++mantissa_;
        if (mantissa_ > maxMantissa)
        {
            mantissa_ /= 10;
            ++exponent_;
        }
    }
    if (exponent_ > maxExponent)
        throw std::overflow_error("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(),
        "ripple::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;
        }
        auto r = g.round();
        if (r == 1 || (r == 0 && (xm & 1) == 1))
        {
            ++xm;
            if (xm > maxMantissa)
            {
                xm /= 10;
                ++xe;
            }
        }
        if (xe > maxExponent)
            throw std::overflow_error("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;
        }
        auto r = g.round();
        if (r == 1 || (r == 0 && (xm & 1) == 1))
        {
            --xm;
            if (xm < minMantissa)
            {
                xm *= 10;
                --xe;
            }
        }
        if (xe < minExponent)
        {
            xm = 0;
            xe = Number{}.exponent_;
        }
    }
    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(),
        "ripple::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;
    auto r = g.round();
    if (r == 1 || (r == 0 && (xm & 1) == 1))
    {
        ++xm;
        if (xm > maxMantissa)
        {
            xm /= 10;
            ++xe;
        }
    }
    if (xe < minExponent)
    {
        xm = 0;
        xe = Number{}.exponent_;
    }
    if (xe > maxExponent)
        throw std::overflow_error(
            "Number::multiplication overflow : exponent is " +
            std::to_string(xe));
    mantissa_ = xm * zn;
    exponent_ = xe;
    XRPL_ASSERT(
        isnormal() || *this == Number{},
        "ripple::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
    const uint128_t 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;
        }
        auto r = g.round();
        if (r == 1 || (r == 0 && (drops & 1) == 1))
        {
            ++drops;
        }
        if (g.is_negative())
            drops = -drops;
    }
    return drops;
}
 
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, "ripple::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,
        "ripple::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,
        "ripple::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 ripple