123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335 |
- // David Eberly, Geometric Tools, Redmond WA 98052
- // Copyright (c) 1998-2020
- // Distributed under the Boost Software License, Version 1.0.
- // https://www.boost.org/LICENSE_1_0.txt
- // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
- // Version: 4.0.2019.11.03
- #pragma once
- #include <Mathematics/BitHacks.h>
- #include <Mathematics/Math.h>
- #include <Mathematics/IEEEBinary.h>
- // The class BSNumber (binary scientific number) is designed to provide exact
- // arithmetic for robust algorithms, typically those for which we need to know
- // the exact sign of determinants. The template parameter UInteger must
- // have support for at least the following public interface. The fstream
- // objects for Write/Read must be created using std::ios::binary. The return
- // value of Write/Read is 'true' iff the operation was successful.
- //
- // class UInteger
- // {
- // public:
- // UInteger();
- // UInteger(UInteger const& number);
- // UInteger(uint32_t number);
- // UInteger(uint64_t number);
- // UInteger& operator=(UInteger const& number);
- // UInteger(UInteger&& number);
- // UInteger& operator=(UInteger&& number);
- // void SetNumBits(int numBits);
- // int32_t GetNumBits() const;
- // bool operator==(UInteger const& number) const;
- // bool operator< (UInteger const& number) const;
- // void Add(UInteger const& n0, UInteger const& n1);
- // void Sub(UInteger const& n0, UInteger const& n1);
- // void Mul(UInteger const& n0, UInteger const& n1);
- // void ShiftLeft(UInteger const& number, int shift);
- // int32_t ShiftRightToOdd(UInteger const& number);
- // int32_t RoundUp();
- // uint64_t GetPrefix(int numRequested) const;
- // bool Write(std::ofstream& output) const;
- // bool Read(std::ifstream& input);
- // };
- //
- // GTEngine currently has 32-bits-per-word storage for UInteger. See the
- // classes UIntegerAP32 (arbitrary precision), UIntegerFP32<N> (fixed
- // precision), and UIntegerALU32 (arithmetic logic unit shared by the previous
- // two classes). The document at the following link describes the design,
- // implementation, and use of BSNumber and BSRational.
- // https://www.geometrictools.com/Documentation/ArbitraryPrecision.pdf
- //
- // Support for debugging algorithms that use exact rational arithmetic. Each
- // BSNumber and BSRational has a double-precision member that is exposed when
- // the conditional define is enabled. Be aware that this can be very slow
- // because of the conversion to double-precision whenever new objects are
- // created by arithmetic operations. As a faster alternative, you can add
- // temporary code in your algorithms that explicitly convert specific rational
- // numbers to double precision.
- //
- //#define GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE
- //
- // Enable this to throw exceptions when NaNs are generated by conversion
- // from BSNumber to float or double.
- #define GTE_THROW_ON_BSNUMBER_ERRORS
- namespace WwiseGTE
- {
- template <typename UInteger> class BSRational;
- template <typename UInteger>
- class BSNumber
- {
- public:
- // Construction. The default constructor generates the zero BSNumber.
- BSNumber()
- :
- mSign(0),
- mBiasedExponent(0)
- {
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = (double)*this;
- #endif
- }
- BSNumber(BSNumber const& number)
- {
- *this = number;
- }
- BSNumber(float number)
- {
- ConvertFrom<IEEEBinary32>(number);
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = (double)*this;
- #endif
- }
- BSNumber(double number)
- {
- ConvertFrom<IEEEBinary64>(number);
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = (double)*this;
- #endif
- }
- BSNumber(int32_t number)
- {
- if (number == 0)
- {
- mSign = 0;
- mBiasedExponent = 0;
- }
- else
- {
- if (number < 0)
- {
- mSign = -1;
- number = -number;
- }
- else
- {
- mSign = 1;
- }
- mBiasedExponent = BitHacks::GetTrailingBit(number);
- mUInteger = (uint32_t)number;
- }
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = (double)*this;
- #endif
- }
- BSNumber(uint32_t number)
- {
- if (number == 0)
- {
- mSign = 0;
- mBiasedExponent = 0;
- }
- else
- {
- mSign = 1;
- mBiasedExponent = BitHacks::GetTrailingBit(number);
- mUInteger = number;
- }
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = (double)*this;
- #endif
- }
- BSNumber(int64_t number)
- {
- if (number == 0)
- {
- mSign = 0;
- mBiasedExponent = 0;
- }
- else
- {
- if (number < 0)
- {
- mSign = -1;
- number = -number;
- }
- else
- {
- mSign = 1;
- }
- mBiasedExponent = BitHacks::GetTrailingBit(number);
- mUInteger = (uint64_t)number;
- }
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = (double)*this;
- #endif
- }
- BSNumber(uint64_t number)
- {
- if (number == 0)
- {
- mSign = 0;
- mBiasedExponent = 0;
- }
- else
- {
- mSign = 1;
- mBiasedExponent = BitHacks::GetTrailingBit(number);
- mUInteger = number;
- }
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = (double)*this;
- #endif
- }
- // The number must be of the form "x" or "+x" or "-x", where x is a
- // positive integer with nonzero leading digit.
- BSNumber(std::string const& number)
- {
- LogAssert(number.size() > 0, "A number must be specified.");
- // Get the leading '+' or '-' if it exists.
- std::string intNumber;
- int sign;
- if (number[0] == '+')
- {
- intNumber = number.substr(1);
- sign = +1;
- LogAssert(intNumber.size() > 1, "Invalid number format.");
- }
- else if (number[0] == '-')
- {
- intNumber = number.substr(1);
- sign = -1;
- LogAssert(intNumber.size() > 1, "Invalid number format.");
- }
- else
- {
- intNumber = number;
- sign = +1;
- }
- *this = ConvertToInteger(intNumber);
- mSign = sign;
- }
- BSNumber(char const* number)
- :
- BSNumber(std::string(number))
- {
- }
- // Implicit conversions. These always use the default rounding mode,
- // round-to-nearest-ties-to-even.
- inline operator float() const
- {
- return ConvertTo<IEEEBinary32>();
- }
- inline operator double() const
- {
- return ConvertTo<IEEEBinary64>();
- }
- // Assignment.
- BSNumber& operator=(BSNumber const& number)
- {
- mSign = number.mSign;
- mBiasedExponent = number.mBiasedExponent;
- mUInteger = number.mUInteger;
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = number.mValue;
- #endif
- return *this;
- }
- // Support for std::move.
- BSNumber(BSNumber&& number)
- {
- *this = std::move(number);
- }
- BSNumber& operator=(BSNumber&& number)
- {
- mSign = number.mSign;
- mBiasedExponent = number.mBiasedExponent;
- mUInteger = std::move(number.mUInteger);
- number.mSign = 0;
- number.mBiasedExponent = 0;
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- mValue = number.mValue;
- #endif
- return *this;
- }
- // Member access.
- inline void SetSign(int32_t sign)
- {
- mSign = sign;
- }
- inline int32_t GetSign() const
- {
- return mSign;
- }
- inline void SetBiasedExponent(int32_t biasedExponent)
- {
- mBiasedExponent = biasedExponent;
- }
- inline int32_t GetBiasedExponent() const
- {
- return mBiasedExponent;
- }
- inline void SetExponent(int32_t exponent)
- {
- mBiasedExponent = exponent - mUInteger.GetNumBits() + 1;
- }
- inline int32_t GetExponent() const
- {
- return mBiasedExponent + mUInteger.GetNumBits() - 1;
- }
- inline UInteger const& GetUInteger() const
- {
- return mUInteger;
- }
- inline UInteger& GetUInteger()
- {
- return mUInteger;
- }
- // Comparisons.
- bool operator==(BSNumber const& number) const
- {
- return (mSign == number.mSign ? EqualIgnoreSign(*this, number) : false);
- }
- bool operator!=(BSNumber const& number) const
- {
- return !operator==(number);
- }
- bool operator< (BSNumber const& number) const
- {
- if (mSign > 0)
- {
- if (number.mSign <= 0)
- {
- return false;
- }
- // Both numbers are positive.
- return LessThanIgnoreSign(*this, number);
- }
- else if (mSign < 0)
- {
- if (number.mSign >= 0)
- {
- return true;
- }
- // Both numbers are negative.
- return LessThanIgnoreSign(number, *this);
- }
- else
- {
- return number.mSign > 0;
- }
- }
- bool operator<=(BSNumber const& number) const
- {
- return !number.operator<(*this);
- }
- bool operator> (BSNumber const& number) const
- {
- return number.operator<(*this);
- }
- bool operator>=(BSNumber const& number) const
- {
- return !operator<(number);
- }
- // Unary operations.
- BSNumber operator+() const
- {
- return *this;
- }
- BSNumber operator-() const
- {
- BSNumber result = *this;
- result.mSign = -result.mSign;
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- result.mValue = (double)result;
- #endif
- return result;
- }
- // Arithmetic.
- BSNumber operator+(BSNumber const& n1) const
- {
- BSNumber const& n0 = *this;
- if (n0.mSign == 0)
- {
- return n1;
- }
- if (n1.mSign == 0)
- {
- return n0;
- }
- if (n0.mSign > 0)
- {
- if (n1.mSign > 0)
- {
- // n0 + n1 = |n0| + |n1|
- return AddIgnoreSign(n0, n1, +1);
- }
- else // n1.mSign < 0
- {
- if (!EqualIgnoreSign(n0, n1))
- {
- if (LessThanIgnoreSign(n1, n0))
- {
- // n0 + n1 = |n0| - |n1| > 0
- return SubIgnoreSign(n0, n1, +1);
- }
- else
- {
- // n0 + n1 = -(|n1| - |n0|) < 0
- return SubIgnoreSign(n1, n0, -1);
- }
- }
- // else n0 + n1 = 0
- }
- }
- else // n0.mSign < 0
- {
- if (n1.mSign < 0)
- {
- // n0 + n1 = -(|n0| + |n1|)
- return AddIgnoreSign(n0, n1, -1);
- }
- else // n1.mSign > 0
- {
- if (!EqualIgnoreSign(n0, n1))
- {
- if (LessThanIgnoreSign(n1, n0))
- {
- // n0 + n1 = -(|n0| - |n1|) < 0
- return SubIgnoreSign(n0, n1, -1);
- }
- else
- {
- // n0 + n1 = |n1| - |n0| > 0
- return SubIgnoreSign(n1, n0, +1);
- }
- }
- // else n0 + n1 = 0
- }
- }
- return BSNumber(); // = 0
- }
- BSNumber operator-(BSNumber const& n1) const
- {
- BSNumber const& n0 = *this;
- if (n0.mSign == 0)
- {
- return -n1;
- }
- if (n1.mSign == 0)
- {
- return n0;
- }
- if (n0.mSign > 0)
- {
- if (n1.mSign < 0)
- {
- // n0 - n1 = |n0| + |n1|
- return AddIgnoreSign(n0, n1, +1);
- }
- else // n1.mSign > 0
- {
- if (!EqualIgnoreSign(n0, n1))
- {
- if (LessThanIgnoreSign(n1, n0))
- {
- // n0 - n1 = |n0| - |n1| > 0
- return SubIgnoreSign(n0, n1, +1);
- }
- else
- {
- // n0 - n1 = -(|n1| - |n0|) < 0
- return SubIgnoreSign(n1, n0, -1);
- }
- }
- // else n0 - n1 = 0
- }
- }
- else // n0.mSign < 0
- {
- if (n1.mSign > 0)
- {
- // n0 - n1 = -(|n0| + |n1|)
- return AddIgnoreSign(n0, n1, -1);
- }
- else // n1.mSign < 0
- {
- if (!EqualIgnoreSign(n0, n1))
- {
- if (LessThanIgnoreSign(n1, n0))
- {
- // n0 - n1 = -(|n0| - |n1|) < 0
- return SubIgnoreSign(n0, n1, -1);
- }
- else
- {
- // n0 - n1 = |n1| - |n0| > 0
- return SubIgnoreSign(n1, n0, +1);
- }
- }
- // else n0 - n1 = 0
- }
- }
- return BSNumber(); // = 0
- }
- BSNumber operator*(BSNumber const& number) const
- {
- BSNumber result; // = 0
- int sign = mSign * number.mSign;
- if (sign != 0)
- {
- result.mSign = sign;
- result.mBiasedExponent = mBiasedExponent + number.mBiasedExponent;
- result.mUInteger.Mul(mUInteger, number.mUInteger);
- }
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- result.mValue = (double)result;
- #endif
- return result;
- }
- BSNumber& operator+=(BSNumber const& number)
- {
- *this = operator+(number);
- return *this;
- }
- BSNumber& operator-=(BSNumber const& number)
- {
- *this = operator-(number);
- return *this;
- }
- BSNumber& operator*=(BSNumber const& number)
- {
- *this = operator*(number);
- return *this;
- }
- // Disk input/output. The fstream objects should be created using
- // std::ios::binary. The return value is 'true' iff the operation
- // was successful.
- bool Write(std::ostream& output) const
- {
- if (output.write((char const*)&mSign, sizeof(mSign)).bad())
- {
- return false;
- }
- if (output.write((char const*)&mBiasedExponent,
- sizeof(mBiasedExponent)).bad())
- {
- return false;
- }
- return mUInteger.Write(output);
- }
- bool Read(std::istream& input)
- {
- if (input.read((char*)&mSign, sizeof(mSign)).bad())
- {
- return false;
- }
- if (input.read((char*)&mBiasedExponent, sizeof(mBiasedExponent)).bad())
- {
- return false;
- }
- return mUInteger.Read(input);
- }
- private:
- // Helper for converting a string to a BSNumber. The string must be
- // valid for a nonnegative integer without a leading '+' sign.
- static BSNumber ConvertToInteger(std::string const& number)
- {
- int digit = static_cast<int>(number.back()) - static_cast<int>('0');
- BSNumber x(digit);
- if (number.size() > 1)
- {
- LogAssert(number.find_first_of("123456789") == 0, "Invalid number format.");
- LogAssert(number.find_first_not_of("0123456789") == std::string::npos, "Invalid number format.");
- BSNumber ten(10), pow10(10);
- for (size_t i = 1, j = number.size() - 2; i < number.size(); ++i, --j)
- {
- digit = static_cast<int>(number[j]) - static_cast<int>('0');
- if (digit > 0)
- {
- x += BSNumber(digit) * pow10;
- }
- pow10 *= ten;
- }
- }
- return x;
- }
- // Helpers for operator==, operator<, operator+, operator-.
- static bool EqualIgnoreSign(BSNumber const& n0, BSNumber const& n1)
- {
- return n0.mBiasedExponent == n1.mBiasedExponent && n0.mUInteger == n1.mUInteger;
- }
- static bool LessThanIgnoreSign(BSNumber const& n0, BSNumber const& n1)
- {
- int32_t e0 = n0.GetExponent(), e1 = n1.GetExponent();
- if (e0 < e1)
- {
- return true;
- }
- if (e0 > e1)
- {
- return false;
- }
- return n0.mUInteger < n1.mUInteger;
- }
- // Add two positive numbers.
- static BSNumber AddIgnoreSign(BSNumber const& n0, BSNumber const& n1, int32_t resultSign)
- {
- BSNumber result, temp;
- int32_t diff = n0.mBiasedExponent - n1.mBiasedExponent;
- if (diff > 0)
- {
- temp.mUInteger.ShiftLeft(n0.mUInteger, diff);
- result.mUInteger.Add(temp.mUInteger, n1.mUInteger);
- result.mBiasedExponent = n1.mBiasedExponent;
- }
- else if (diff < 0)
- {
- temp.mUInteger.ShiftLeft(n1.mUInteger, -diff);
- result.mUInteger.Add(n0.mUInteger, temp.mUInteger);
- result.mBiasedExponent = n0.mBiasedExponent;
- }
- else
- {
- temp.mUInteger.Add(n0.mUInteger, n1.mUInteger);
- int32_t shift = result.mUInteger.ShiftRightToOdd(temp.mUInteger);
- result.mBiasedExponent = n0.mBiasedExponent + shift;
- }
- result.mSign = resultSign;
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- result.mValue = (double)result;
- #endif
- return result;
- }
- // Subtract two positive numbers where n0 > n1.
- static BSNumber SubIgnoreSign(BSNumber const& n0, BSNumber const& n1, int32_t resultSign)
- {
- BSNumber result, temp;
- int32_t diff = n0.mBiasedExponent - n1.mBiasedExponent;
- if (diff > 0)
- {
- temp.mUInteger.ShiftLeft(n0.mUInteger, diff);
- result.mUInteger.Sub(temp.mUInteger, n1.mUInteger);
- result.mBiasedExponent = n1.mBiasedExponent;
- }
- else if (diff < 0)
- {
- temp.mUInteger.ShiftLeft(n1.mUInteger, -diff);
- result.mUInteger.Sub(n0.mUInteger, temp.mUInteger);
- result.mBiasedExponent = n0.mBiasedExponent;
- }
- else
- {
- temp.mUInteger.Sub(n0.mUInteger, n1.mUInteger);
- int32_t shift = result.mUInteger.ShiftRightToOdd(temp.mUInteger);
- result.mBiasedExponent = n0.mBiasedExponent + shift;
- }
- result.mSign = resultSign;
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- result.mValue = (double)result;
- #endif
- return result;
- }
- // Support for conversions from floating-point numbers to BSNumber.
- template <typename IEEE>
- void ConvertFrom(typename IEEE::FloatType number)
- {
- IEEE x(number);
- // Extract sign s, biased exponent e, and trailing significand t.
- typename IEEE::UIntType s = x.GetSign();
- typename IEEE::UIntType e = x.GetBiased();
- typename IEEE::UIntType t = x.GetTrailing();
- if (e == 0)
- {
- if (t == 0) // zeros
- {
- // x = (-1)^s * 0
- mSign = 0;
- mBiasedExponent = 0;
- }
- else // subnormal numbers
- {
- // x = (-1)^s * 0.t * 2^{1-EXPONENT_BIAS}
- int32_t last = BitHacks::GetTrailingBit(t);
- int32_t diff = IEEE::NUM_TRAILING_BITS - last;
- mSign = (s > 0 ? -1 : 1);
- mBiasedExponent = IEEE::MIN_SUB_EXPONENT - diff;
- mUInteger = (t >> last);
- }
- }
- else if (e < IEEE::MAX_BIASED_EXPONENT) // normal numbers
- {
- // x = (-1)^s * 1.t * 2^{e-EXPONENT_BIAS}
- if (t > 0)
- {
- int32_t last = BitHacks::GetTrailingBit(t);
- int32_t diff = IEEE::NUM_TRAILING_BITS - last;
- mSign = (s > 0 ? -1 : 1);
- mBiasedExponent =
- static_cast<int32_t>(e) - IEEE::EXPONENT_BIAS - diff;
- mUInteger = ((t | IEEE::SUP_TRAILING) >> last);
- }
- else
- {
- mSign = (s > 0 ? -1 : 1);
- mBiasedExponent = static_cast<int32_t>(e) - IEEE::EXPONENT_BIAS;
- mUInteger = (typename IEEE::UIntType)1;
- }
- }
- else // e == MAX_BIASED_EXPONENT, special numbers
- {
- if (t == 0) // infinities
- {
- // x = (-1)^s * infinity
- #if defined(GTE_THROW_ON_BSNUMBER_ERRORS)
- // TODO: This should not be a warning, but BSNumber does
- // not have a representation for +infinity or -infinity.
- // Consider doing so with mSign in {-2,2} and all other
- // members zero-valued.
- LogWarning("Input is " + std::string(s > 0 ? "-" : "+") + "infinity.");
- #else
- // Return (-1)^s * 2^{1+EXPONENT_BIAS} for a graceful
- // exit.
- mSign = (s > 0 ? -1 : 1);
- mBiasedExponent = 1 + IEEE::EXPONENT_BIAS;
- mUInteger = (typename IEEE::UIntType)1;
- #endif
- }
- else // not-a-number (NaN)
- {
- #if defined(GTE_THROW_ON_BSNUMBER_ERRORS)
- // TODO: BSNumber does not have a representation for
- // NaNs. Consider doing so with mSign in {-3,3} and a
- // payload stored in mBits.
- LogError("Input is a " +
- std::string(t & IEEE::NAN_QUIET_MASK ?
- "quiet" : "signaling") + " NaN with payload " +
- std::to_string(t & IEEE::NAN_PAYLOAD_MASK) + ".");
- #else
- // Return 0 for a graceful exit.
- mSign = 0;
- mBiasedExponent = 0;
- #endif
- }
- }
- }
- // Support for conversions from BSNumber to floating-point numbers.
- template <typename IEEE>
- typename IEEE::FloatType ConvertTo() const
- {
- typename IEEE::UIntType s = (mSign < 0 ? 1 : 0);
- typename IEEE::UIntType e, t;
- if (mSign != 0)
- {
- // The conversions use round-to-nearest-ties-to-even
- // semantics.
- int32_t exponent = GetExponent();
- if (exponent < IEEE::MIN_EXPONENT)
- {
- if (exponent < IEEE::MIN_EXPONENT - 1
- || mUInteger.GetNumBits() == 1)
- {
- // x = 1.0*2^{MIN_EXPONENT-1}, round to zero.
- e = 0;
- t = 0;
- }
- else
- {
- // Round to min subnormal.
- e = 0;
- t = 1;
- }
- }
- else if (exponent < IEEE::MIN_SUB_EXPONENT)
- {
- // The second input is in {0, ..., NUM_TRAILING_BITS-1}.
- t = GetTrailing<IEEE>(0, IEEE::MIN_SUB_EXPONENT - exponent - 1);
- if (t & IEEE::SUP_TRAILING)
- {
- // Leading NUM_SIGNIFICAND_BITS bits were all 1, so
- // round to min normal.
- e = 1;
- t = 0;
- }
- else
- {
- e = 0;
- }
- }
- else if (exponent <= IEEE::EXPONENT_BIAS)
- {
- e = static_cast<uint32_t>(exponent + IEEE::EXPONENT_BIAS);
- t = GetTrailing<IEEE>(1, 0);
- if (t & (IEEE::SUP_TRAILING << 1))
- {
- // Carry-out occurred, so increase exponent by 1 and
- // shift right to compensate.
- ++e;
- t >>= 1;
- }
- // Eliminate the leading 1 (implied for normals).
- t &= ~IEEE::SUP_TRAILING;
- }
- else
- {
- // Set to infinity.
- e = IEEE::MAX_BIASED_EXPONENT;
- t = 0;
- }
- }
- else
- {
- // The input is zero.
- e = 0;
- t = 0;
- }
- IEEE x(s, e, t);
- return x.number;
- }
- template <typename IEEE>
- typename IEEE::UIntType GetTrailing(int32_t normal, int32_t sigma) const
- {
- int32_t const numRequested = IEEE::NUM_SIGNIFICAND_BITS + normal;
- // We need numRequested bits to determine rounding direction.
- // These are stored in the high-order bits of 'prefix'.
- uint64_t prefix = mUInteger.GetPrefix(numRequested);
- // The first bit index after the implied binary point for
- // rounding.
- int32_t diff = numRequested - sigma;
- int32_t roundBitIndex = 64 - diff;
- // Determine value based on round-to-nearest-ties-to-even.
- uint64_t mask = (1ull << roundBitIndex);
- uint64_t round;
- if (prefix & mask)
- {
- // The first bit of the remainder is 1.
- if (mUInteger.GetNumBits() == diff)
- {
- // The first bit of the remainder is the lowest-order bit
- // of mBits[0]. Apply the ties-to-even rule.
- if (prefix & (mask << 1))
- {
- // The last bit of the trailing significand is odd,
- // so round up.
- round = 1;
- }
- else
- {
- // The last bit of the trailing significand is even,
- // so round down.
- round = 0;
- }
- }
- else
- {
- // The first bit of the remainder is not the lowest-order
- // bit of mBits[0]. The remainder as a fraction is larger
- // than 1/2, so round up.
- round = 1;
- }
- }
- else
- {
- // The first bit of the remainder is 0, so round down.
- round = 0;
- }
- // Get the unrounded trailing significand.
- uint64_t trailing = prefix >> (roundBitIndex + 1);
- // Apply the rounding.
- trailing += round;
- return static_cast<typename IEEE::UIntType>(trailing);
- }
- #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
- public:
- // List this first so that it shows up first in the debugger watch
- // window.
- double mValue;
- private:
- #endif
- // The number 0 is represented by: mSign = 0, mBiasedExponent = 0 and
- // mUInteger = 0. For nonzero numbers, mSign != 0 and mUInteger > 0.
- int32_t mSign;
- int32_t mBiasedExponent;
- UInteger mUInteger;
- // BSRational depends on the design of BSNumber, so allow it to have
- // full access to the implementation.
- friend class BSRational<UInteger>;
- };
- // Explicit conversion to a user-specified precision. The rounding
- // mode is one of the flags provided in <cfenv>. The modes are
- // FE_TONEAREST: round to nearest ties to even
- // FE_DOWNWARD: round towards negative infinity
- // FE_TOWARDZERO: round towards zero
- // FE_UPWARD: round towards positive infinity
- template <typename UInteger>
- void Convert(BSNumber<UInteger> const& input, int32_t precision,
- int32_t roundingMode, BSNumber<UInteger>& output)
- {
- if (precision <= 0)
- {
- LogError("Precision must be positive.");
- }
- int64_t const maxSize = static_cast<int64_t>(UInteger::GetMaxSize());
- int64_t const excess = 32LL * maxSize - static_cast<int64_t>(precision);
- if (excess <= 0)
- {
- LogError("The maximum precision has been exceeded.");
- }
- if (input.GetSign() == 0)
- {
- output = BSNumber<UInteger>(0);
- return;
- }
- // Let p = precision and n+1 be the number of bits of the input.
- // Compute n+1-p. If it is nonpositive, then the requested precision
- // is already satisfied by the input.
- int32_t np1mp = input.GetUInteger().GetNumBits() - precision;
- if (np1mp <= 0)
- {
- output = input;
- return;
- }
- // At this point, the requested number of bits is smaller than the
- // number of bits in the input. Round the input to the smaller number
- // of bits using the specified rounding mode.
- UInteger& outW = output.GetUInteger();
- outW.SetNumBits(precision);
- outW.SetAllBitsToZero();
- int32_t const outSize = outW.GetSize();
- int32_t const precisionM1 = precision - 1;
- int32_t const outLeading = precisionM1 % 32;
- uint32_t outMask = (1 << outLeading);
- auto& outBits = outW.GetBits();
- int32_t outCurrent = outSize - 1;
- UInteger const& inW = input.GetUInteger();
- int32_t const inSize = inW.GetSize();
- int32_t const inLeading = (inW.GetNumBits() - 1) % 32;
- uint32_t inMask = (1 << inLeading);
- auto const& inBits = inW.GetBits();
- int32_t inCurrent = inSize - 1;
- int32_t lastBit = -1;
- for (int i = precisionM1; i >= 0; --i)
- {
- if (inBits[inCurrent] & inMask)
- {
- outBits[outCurrent] |= outMask;
- lastBit = 1;
- }
- else
- {
- lastBit = 0;
- }
- if (inMask == 0x00000001u)
- {
- --inCurrent;
- inMask = 0x80000000u;
- }
- else
- {
- inMask >>= 1;
- }
- if (outMask == 0x00000001u)
- {
- --outCurrent;
- outMask = 0x80000000u;
- }
- else
- {
- outMask >>= 1;
- }
- }
- // At this point as a sequence of bits, r = u_{n-p} ... u_0.
- int32_t sign = input.GetSign();
- int32_t outExponent = input.GetExponent();
- if (roundingMode == FE_TONEAREST)
- {
- // Determine whether u_{n-p} is positive.
- uint32_t positive = (inBits[inCurrent] & inMask) != 0u;
- if (positive && (np1mp > 1 || lastBit == 1))
- {
- // round up
- outExponent += outW.RoundUp();
- }
- // else round down, equivalent to truncating the r bits
- }
- else if (roundingMode == FE_UPWARD)
- {
- // The remainder r must be positive because n-p >= 0 and u_0 = 1.
- if (sign > 0)
- {
- // round up
- outExponent += outW.RoundUp();
- }
- // else round down, equivalent to truncating the r bits
- }
- else if (roundingMode == FE_DOWNWARD)
- {
- // The remainder r must be positive because n-p >= 0 and u_0 = 1.
- if (sign < 0)
- {
- // Round down. This is the round-up operation applied to
- // w, but the final sign is negative which amounts to
- // rounding down.
- outExponent += outW.RoundUp();
- }
- // else round down, equivalent to truncating the r bits
- }
- else if (roundingMode != FE_TOWARDZERO)
- {
- // Currently, no additional implementation-dependent modes
- // are supported for rounding.
- LogError("Implementation-dependent rounding mode not supported.");
- }
- // else roundingMode == FE_TOWARDZERO. Truncate the r bits, which
- // requires no additional work.
- output.SetSign(sign);
- output.SetBiasedExponent(outExponent - precisionM1);
- }
- }
- namespace std
- {
- // TODO: Allow for implementations of the math functions in which a
- // specified precision is used when computing the result.
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> acos(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::acos((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> acosh(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::acosh((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> asin(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::asin((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> asinh(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::asinh((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> atan(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::atan((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> atanh(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::atanh((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> atan2(WwiseGTE::BSNumber<UInteger> const& y, WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::atan2((double)y, (double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> ceil(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::ceil((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> cos(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::cos((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> cosh(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::cosh((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> exp(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::exp((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> exp2(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::exp2((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> fabs(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (x.GetSign() >= 0 ? x : -x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> floor(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::floor((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> fmod(WwiseGTE::BSNumber<UInteger> const& x, WwiseGTE::BSNumber<UInteger> const& y)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::fmod((double)x, (double)y);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> frexp(WwiseGTE::BSNumber<UInteger> const& x, int* exponent)
- {
- if (x.GetSign() != 0)
- {
- WwiseGTE::BSNumber<UInteger> result = x;
- *exponent = result.GetExponent() + 1;
- result.SetExponent(-1);
- return result;
- }
- else
- {
- *exponent = 0;
- return WwiseGTE::BSNumber<UInteger>(0);
- }
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> ldexp(WwiseGTE::BSNumber<UInteger> const& x, int exponent)
- {
- WwiseGTE::BSNumber<UInteger> result = x;
- result.SetBiasedExponent(result.GetBiasedExponent() + exponent);
- return result;
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> log(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::log((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> log2(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::log2((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> log10(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::log10((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> pow(WwiseGTE::BSNumber<UInteger> const& x, WwiseGTE::BSNumber<UInteger> const& y)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::pow((double)x, (double)y);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> sin(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::sin((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> sinh(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::sinh((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> sqrt(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::sqrt((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> tan(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::tan((double)x);
- }
- template <typename UInteger>
- inline WwiseGTE::BSNumber<UInteger> tanh(WwiseGTE::BSNumber<UInteger> const& x)
- {
- return (WwiseGTE::BSNumber<UInteger>)std::tanh((double)x);
- }
- // Type trait that says BSNumber is a signed type.
- template <typename UInteger>
- struct is_signed<WwiseGTE::BSNumber<UInteger>> : true_type {};
- }
- namespace WwiseGTE
- {
- template <typename UInteger>
- inline BSNumber<UInteger> atandivpi(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)atandivpi((double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> atan2divpi(BSNumber<UInteger> const& y, BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)atan2divpi((double)y, (double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> clamp(BSNumber<UInteger> const& x, BSNumber<UInteger> const& xmin, BSNumber<UInteger> const& xmax)
- {
- return (BSNumber<UInteger>)clamp((double)x, (double)xmin, (double)xmax);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> cospi(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)cospi((double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> exp10(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)exp10((double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> invsqrt(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)invsqrt((double)x);
- }
- template <typename UInteger>
- inline int isign(BSNumber<UInteger> const& x)
- {
- return isign((double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> saturate(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)saturate((double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> sign(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)sign((double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> sinpi(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)sinpi((double)x);
- }
- template <typename UInteger>
- inline BSNumber<UInteger> sqr(BSNumber<UInteger> const& x)
- {
- return (BSNumber<UInteger>)sqr((double)x);
- }
- // See the comments in GteMath.h about trait is_arbitrary_precision.
- template <typename UInteger>
- struct is_arbitrary_precision_internal<BSNumber<UInteger>> : std::true_type {};
- }
|