mirror of
https://github.com/XRPLF/rippled.git
synced 2026-04-29 15:37:57 +00:00
849 lines
20 KiB
C++
849 lines
20 KiB
C++
#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)
|
||
static_assert(std::is_same_v<uint128_t, ripple::numberuint128>);
|
||
|
||
namespace std {
|
||
|
||
template <>
|
||
struct make_unsigned<ripple::numberint128>
|
||
{
|
||
using type = uint128_t;
|
||
};
|
||
|
||
} // namespace std
|
||
|
||
namespace ripple {
|
||
|
||
thread_local Number::rounding_mode Number::mode_ = Number::to_nearest;
|
||
// TODO: Once the Rules switching is implemented, default to largeRange
|
||
thread_local std::reference_wrapper<MantissaRange const> Number::range_ =
|
||
smallRange; // largeRange;
|
||
|
||
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
|
||
template <class T>
|
||
void
|
||
push(T 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;
|
||
|
||
private:
|
||
void
|
||
doPush(unsigned d) 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::doPush(unsigned d) noexcept
|
||
{
|
||
xbit_ = xbit_ || (digits_ & 0x0000'0000'0000'000F) != 0;
|
||
digits_ >>= 4;
|
||
digits_ |= (d & 0x0000'0000'0000'000FULL) << 60;
|
||
}
|
||
|
||
template <class T>
|
||
inline void
|
||
Number::Guard::push(T d) noexcept
|
||
{
|
||
doPush(static_cast<unsigned>(d));
|
||
}
|
||
|
||
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
|
||
Number::oneSmall()
|
||
{
|
||
return Number{
|
||
Number::smallRange.min, -Number::smallRange.log, Number::unchecked{}};
|
||
};
|
||
|
||
constexpr Number oneSml = Number::oneSmall();
|
||
|
||
constexpr Number
|
||
Number::oneLarge()
|
||
{
|
||
return Number{
|
||
Number::largeRange.min, -Number::largeRange.log, Number::unchecked{}};
|
||
};
|
||
|
||
constexpr Number oneLrg = Number::oneLarge();
|
||
|
||
Number
|
||
Number::one()
|
||
{
|
||
if (&range_.get() == &smallRange)
|
||
return oneSml;
|
||
XRPL_ASSERT(&range_.get() == &largeRange, "Number::one() : valid range_");
|
||
return oneLrg;
|
||
}
|
||
|
||
// Use the member names in this static function for now so the diff is cleaner
|
||
void
|
||
Number::normalize(
|
||
internalrep& mantissa_,
|
||
int& exponent_,
|
||
internalrep const& minMantissa,
|
||
internalrep const& maxMantissa)
|
||
{
|
||
constexpr Number zero = Number{};
|
||
if (mantissa_ == 0)
|
||
{
|
||
mantissa_ = zero.mantissa_;
|
||
exponent_ = zero.exponent_;
|
||
return;
|
||
}
|
||
bool const negative = (mantissa_ < 0);
|
||
auto m = static_cast<std::make_unsigned_t<internalrep>>(
|
||
negative ? -mantissa_ : mantissa_);
|
||
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))
|
||
{
|
||
mantissa_ = zero.mantissa_;
|
||
exponent_ = zero.exponent_;
|
||
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_;
|
||
}
|
||
|
||
void
|
||
Number::normalize()
|
||
{
|
||
normalize(
|
||
mantissa_, exponent_, Number::minMantissa(), Number::maxMantissa());
|
||
}
|
||
|
||
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)
|
||
{
|
||
auto const maxMantissa = Number::maxMantissa();
|
||
|
||
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
|
||
{
|
||
auto const minMantissa = Number::minMantissa();
|
||
|
||
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();
|
||
|
||
auto const maxMantissa = Number::maxMantissa();
|
||
|
||
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<internalrep>(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
|
||
// TODO: Can/should this be made bigger for largeRange?
|
||
// log(2^127,10) ~ 38.2
|
||
// largeRange.log = 18
|
||
// f can be up to 10^(37-18) = 10^19 safely
|
||
constexpr uint128_t f = 100'000'000'000'000'000;
|
||
static_assert(f == smallRange.min * 100);
|
||
static_assert(smallRange.log == 15);
|
||
mantissa_ = uint128_t(nm) * f / uint128_t(dm);
|
||
exponent_ = ne - de - 17;
|
||
mantissa_ *= np * dp;
|
||
normalize();
|
||
return *this;
|
||
}
|
||
|
||
Number::operator rep() const
|
||
{
|
||
internalrep 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;
|
||
}
|
||
if (offset == 0 && drops > std::numeric_limits<rep>::max())
|
||
{
|
||
// If offset == 0, then the loop won't run, and the overflow check
|
||
// won't be made, but a int128 can overflow int64 by itself, so
|
||
// check here.
|
||
throw std::overflow_error("Number::operator rep() overflow");
|
||
}
|
||
for (; offset > 0; --offset)
|
||
{
|
||
if (drops > std::numeric_limits<rep>::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 static_cast<rep>(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();
|
||
|
||
bool const negative = amount.mantissa() < 0;
|
||
|
||
auto const mantissa = [&]() {
|
||
auto mantissa = amount.mantissa();
|
||
if (negative)
|
||
{
|
||
mantissa = -mantissa;
|
||
}
|
||
XRPL_ASSERT(
|
||
mantissa < std::numeric_limits<std::uint64_t>::max(),
|
||
"ripple::to_string(Number) : mantissa fits in uin64_t");
|
||
return static_cast<std::uint64_t>(mantissa);
|
||
}();
|
||
|
||
// Use scientific notation for exponents that are too small or too large
|
||
auto const rangeLog = Number::mantissaLog();
|
||
if (((exponent != 0) &&
|
||
((exponent < -(rangeLog + 10)) || (exponent > -(rangeLog - 10)))))
|
||
{
|
||
std::string ret = negative ? "-" : "";
|
||
ret.append(std::to_string(mantissa));
|
||
ret.append(1, 'e');
|
||
ret.append(std::to_string(exponent));
|
||
return ret;
|
||
}
|
||
|
||
// TODO: These numbers are probably wrong for largeRange
|
||
XRPL_ASSERT(
|
||
exponent + 43 > 0, "ripple::to_string(Number) : minimum exponent");
|
||
|
||
ptrdiff_t const pad_prefix = Number::mantissaLog() + 12;
|
||
ptrdiff_t const pad_suffix = Number::mantissaLog() + 8;
|
||
|
||
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 + pad_prefix + 16);
|
||
|
||
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 Number::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)
|
||
{
|
||
auto const one = Number::one();
|
||
|
||
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)
|
||
{
|
||
auto const one = Number::one();
|
||
|
||
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)
|
||
{
|
||
auto const one = Number::one();
|
||
|
||
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
|