ARCHIVE/xahaud
3
0
Fork 0
mirror of https://github.com/Xahau/xahaud.git synced 2025-12-06 17:27:52 +00:00
Code Issues Packages Projects Releases Wiki Activity

Add tests

Browse Source
This commit is contained in:
Howard Hinnant
2022-04-21 16:47:07 -04:00
committed by Elliot Lee
parent c9c54c9799
commit 48e804c40c
3 changed files with 447 additions and 142 deletions
Download Patch File Download Diff File Expand all files Collapse all files

112
src/ripple/basics/impl/Number.cpp
View File

@@ -33,53 +33,66 @@ using uint128_t = __uint128_t;
namespace ripple {
// guard
// Guard
class Number::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_;
std::uint8_t xbit_ : 1;
std::uint8_t sbit_ : 1; // TODO : get rid of
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}
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
Number::Guard::set_positive() noexcept
{
sbit_ = 0;
}
inline void
Number::guard::set_negative() noexcept
Number::Guard::set_negative() noexcept
{
sbit_ = 1;
}
inline bool
Number::guard::is_negative() const noexcept
Number::Guard::is_negative() const noexcept
{
return sbit_ == 1;
}
inline void
Number::guard::push(unsigned d) noexcept
Number::Guard::push(unsigned d) noexcept
{
xbit_ = xbit_ || (digits_ & 0x0000'0000'0000'000F) != 0;
digits_ >>= 4;
@@ -87,7 +100,7 @@ Number::guard::push(unsigned d) noexcept
}
inline unsigned
Number::guard::pop() noexcept
Number::Guard::pop() noexcept
{
unsigned d = (digits_ & 0xF000'0000'0000'0000) >> 60;
digits_ <<= 4;
@@ -95,7 +108,7 @@ Number::guard::pop() noexcept
}
int
Number::guard::round() noexcept
Number::Guard::round() noexcept
{
if (digits_ > 0x5000'0000'0000'0000)
return 1;
@@ -127,10 +140,12 @@ Number::normalize()
m *= 10;
--exponent_;
}
Guard g;
while (m > maxMantissa)
{
if (exponent_ >= maxExponent)
throw std::overflow_error("Number::normalize 1");
g.push(m % 10);
m /= 10;
++exponent_;
}
@@ -141,6 +156,16 @@ Number::normalize()
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");
@@ -180,7 +205,7 @@ Number::operator+=(Number const& y)
ym = -ym;
yn = -1;
}
guard g;
Guard g;
if (xe < ye)
{
if (xn == -1)
@@ -261,7 +286,6 @@ Number::operator+=(Number const& y)
}
mantissa_ = xm * xn;
exponent_ = xe;
assert(isnormal());
return *this;
}
@@ -295,7 +319,7 @@ Number::operator*=(Number const& y)
auto zm = uint128_t(xm) * uint128_t(ym);
auto ze = xe + ye;
auto zn = xn * yn;
guard g;
Guard g;
while (zm > maxMantissa)
{
g.push(static_cast<unsigned>(zm % 10));
@@ -379,7 +403,7 @@ Number::operator rep() const
{
rep drops = mantissa_;
int offset = exponent_;
guard g;
Guard g;
if (drops != 0)
{
if (drops < 0)
@@ -395,7 +419,7 @@ Number::operator rep() const
for (; offset > 0; --offset)
{
if (drops > std::numeric_limits<decltype(drops)>::max() / 10)
throw std::runtime_error("Number::operator rep() overflow");
throw std::overflow_error("Number::operator rep() overflow");
drops *= 10;
}
auto r = g.round();
@@ -505,10 +529,10 @@ to_string(Number const& amount)
}
// Returns f^n
// Uses a log_2(n) number of mulitiplications
// Uses a log_2(n) number of multiplications
Number
power(Number f, unsigned n)
power(Number const& f, unsigned n)
{
if (n == 0)
return one;
@@ -525,6 +549,11 @@ power(Number f, unsigned n)
// 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)
{
@@ -590,10 +619,48 @@ root(Number f, unsigned d)
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 f, unsigned n, unsigned d)
power(Number const& f, unsigned n, unsigned d)
{
if (f == one)
return f;
@@ -606,9 +673,8 @@ power(Number f, unsigned n, unsigned d)
return one;
if (abs(f) < one)
return Number{};
if (abs(f) > one)
throw std::overflow_error("Number::power infinity");
throw std::overflow_error("Number::power nan");
// abs(f) > one
throw std::overflow_error("Number::power infinity");
}
if (n == 0)
return one;