BSNumber.h 44 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335
  1. // David Eberly, Geometric Tools, Redmond WA 98052
  2. // Copyright (c) 1998-2020
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // https://www.boost.org/LICENSE_1_0.txt
  5. // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
  6. // Version: 4.0.2019.11.03
  7. #pragma once
  8. #include <Mathematics/BitHacks.h>
  9. #include <Mathematics/Math.h>
  10. #include <Mathematics/IEEEBinary.h>
  11. // The class BSNumber (binary scientific number) is designed to provide exact
  12. // arithmetic for robust algorithms, typically those for which we need to know
  13. // the exact sign of determinants. The template parameter UInteger must
  14. // have support for at least the following public interface. The fstream
  15. // objects for Write/Read must be created using std::ios::binary. The return
  16. // value of Write/Read is 'true' iff the operation was successful.
  17. //
  18. // class UInteger
  19. // {
  20. // public:
  21. // UInteger();
  22. // UInteger(UInteger const& number);
  23. // UInteger(uint32_t number);
  24. // UInteger(uint64_t number);
  25. // UInteger& operator=(UInteger const& number);
  26. // UInteger(UInteger&& number);
  27. // UInteger& operator=(UInteger&& number);
  28. // void SetNumBits(int numBits);
  29. // int32_t GetNumBits() const;
  30. // bool operator==(UInteger const& number) const;
  31. // bool operator< (UInteger const& number) const;
  32. // void Add(UInteger const& n0, UInteger const& n1);
  33. // void Sub(UInteger const& n0, UInteger const& n1);
  34. // void Mul(UInteger const& n0, UInteger const& n1);
  35. // void ShiftLeft(UInteger const& number, int shift);
  36. // int32_t ShiftRightToOdd(UInteger const& number);
  37. // int32_t RoundUp();
  38. // uint64_t GetPrefix(int numRequested) const;
  39. // bool Write(std::ofstream& output) const;
  40. // bool Read(std::ifstream& input);
  41. // };
  42. //
  43. // GTEngine currently has 32-bits-per-word storage for UInteger. See the
  44. // classes UIntegerAP32 (arbitrary precision), UIntegerFP32<N> (fixed
  45. // precision), and UIntegerALU32 (arithmetic logic unit shared by the previous
  46. // two classes). The document at the following link describes the design,
  47. // implementation, and use of BSNumber and BSRational.
  48. // https://www.geometrictools.com/Documentation/ArbitraryPrecision.pdf
  49. //
  50. // Support for debugging algorithms that use exact rational arithmetic. Each
  51. // BSNumber and BSRational has a double-precision member that is exposed when
  52. // the conditional define is enabled. Be aware that this can be very slow
  53. // because of the conversion to double-precision whenever new objects are
  54. // created by arithmetic operations. As a faster alternative, you can add
  55. // temporary code in your algorithms that explicitly convert specific rational
  56. // numbers to double precision.
  57. //
  58. //#define GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE
  59. //
  60. // Enable this to throw exceptions when NaNs are generated by conversion
  61. // from BSNumber to float or double.
  62. #define GTE_THROW_ON_BSNUMBER_ERRORS
  63. namespace WwiseGTE
  64. {
  65. template <typename UInteger> class BSRational;
  66. template <typename UInteger>
  67. class BSNumber
  68. {
  69. public:
  70. // Construction. The default constructor generates the zero BSNumber.
  71. BSNumber()
  72. :
  73. mSign(0),
  74. mBiasedExponent(0)
  75. {
  76. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  77. mValue = (double)*this;
  78. #endif
  79. }
  80. BSNumber(BSNumber const& number)
  81. {
  82. *this = number;
  83. }
  84. BSNumber(float number)
  85. {
  86. ConvertFrom<IEEEBinary32>(number);
  87. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  88. mValue = (double)*this;
  89. #endif
  90. }
  91. BSNumber(double number)
  92. {
  93. ConvertFrom<IEEEBinary64>(number);
  94. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  95. mValue = (double)*this;
  96. #endif
  97. }
  98. BSNumber(int32_t number)
  99. {
  100. if (number == 0)
  101. {
  102. mSign = 0;
  103. mBiasedExponent = 0;
  104. }
  105. else
  106. {
  107. if (number < 0)
  108. {
  109. mSign = -1;
  110. number = -number;
  111. }
  112. else
  113. {
  114. mSign = 1;
  115. }
  116. mBiasedExponent = BitHacks::GetTrailingBit(number);
  117. mUInteger = (uint32_t)number;
  118. }
  119. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  120. mValue = (double)*this;
  121. #endif
  122. }
  123. BSNumber(uint32_t number)
  124. {
  125. if (number == 0)
  126. {
  127. mSign = 0;
  128. mBiasedExponent = 0;
  129. }
  130. else
  131. {
  132. mSign = 1;
  133. mBiasedExponent = BitHacks::GetTrailingBit(number);
  134. mUInteger = number;
  135. }
  136. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  137. mValue = (double)*this;
  138. #endif
  139. }
  140. BSNumber(int64_t number)
  141. {
  142. if (number == 0)
  143. {
  144. mSign = 0;
  145. mBiasedExponent = 0;
  146. }
  147. else
  148. {
  149. if (number < 0)
  150. {
  151. mSign = -1;
  152. number = -number;
  153. }
  154. else
  155. {
  156. mSign = 1;
  157. }
  158. mBiasedExponent = BitHacks::GetTrailingBit(number);
  159. mUInteger = (uint64_t)number;
  160. }
  161. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  162. mValue = (double)*this;
  163. #endif
  164. }
  165. BSNumber(uint64_t number)
  166. {
  167. if (number == 0)
  168. {
  169. mSign = 0;
  170. mBiasedExponent = 0;
  171. }
  172. else
  173. {
  174. mSign = 1;
  175. mBiasedExponent = BitHacks::GetTrailingBit(number);
  176. mUInteger = number;
  177. }
  178. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  179. mValue = (double)*this;
  180. #endif
  181. }
  182. // The number must be of the form "x" or "+x" or "-x", where x is a
  183. // positive integer with nonzero leading digit.
  184. BSNumber(std::string const& number)
  185. {
  186. LogAssert(number.size() > 0, "A number must be specified.");
  187. // Get the leading '+' or '-' if it exists.
  188. std::string intNumber;
  189. int sign;
  190. if (number[0] == '+')
  191. {
  192. intNumber = number.substr(1);
  193. sign = +1;
  194. LogAssert(intNumber.size() > 1, "Invalid number format.");
  195. }
  196. else if (number[0] == '-')
  197. {
  198. intNumber = number.substr(1);
  199. sign = -1;
  200. LogAssert(intNumber.size() > 1, "Invalid number format.");
  201. }
  202. else
  203. {
  204. intNumber = number;
  205. sign = +1;
  206. }
  207. *this = ConvertToInteger(intNumber);
  208. mSign = sign;
  209. }
  210. BSNumber(char const* number)
  211. :
  212. BSNumber(std::string(number))
  213. {
  214. }
  215. // Implicit conversions. These always use the default rounding mode,
  216. // round-to-nearest-ties-to-even.
  217. inline operator float() const
  218. {
  219. return ConvertTo<IEEEBinary32>();
  220. }
  221. inline operator double() const
  222. {
  223. return ConvertTo<IEEEBinary64>();
  224. }
  225. // Assignment.
  226. BSNumber& operator=(BSNumber const& number)
  227. {
  228. mSign = number.mSign;
  229. mBiasedExponent = number.mBiasedExponent;
  230. mUInteger = number.mUInteger;
  231. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  232. mValue = number.mValue;
  233. #endif
  234. return *this;
  235. }
  236. // Support for std::move.
  237. BSNumber(BSNumber&& number)
  238. {
  239. *this = std::move(number);
  240. }
  241. BSNumber& operator=(BSNumber&& number)
  242. {
  243. mSign = number.mSign;
  244. mBiasedExponent = number.mBiasedExponent;
  245. mUInteger = std::move(number.mUInteger);
  246. number.mSign = 0;
  247. number.mBiasedExponent = 0;
  248. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  249. mValue = number.mValue;
  250. #endif
  251. return *this;
  252. }
  253. // Member access.
  254. inline void SetSign(int32_t sign)
  255. {
  256. mSign = sign;
  257. }
  258. inline int32_t GetSign() const
  259. {
  260. return mSign;
  261. }
  262. inline void SetBiasedExponent(int32_t biasedExponent)
  263. {
  264. mBiasedExponent = biasedExponent;
  265. }
  266. inline int32_t GetBiasedExponent() const
  267. {
  268. return mBiasedExponent;
  269. }
  270. inline void SetExponent(int32_t exponent)
  271. {
  272. mBiasedExponent = exponent - mUInteger.GetNumBits() + 1;
  273. }
  274. inline int32_t GetExponent() const
  275. {
  276. return mBiasedExponent + mUInteger.GetNumBits() - 1;
  277. }
  278. inline UInteger const& GetUInteger() const
  279. {
  280. return mUInteger;
  281. }
  282. inline UInteger& GetUInteger()
  283. {
  284. return mUInteger;
  285. }
  286. // Comparisons.
  287. bool operator==(BSNumber const& number) const
  288. {
  289. return (mSign == number.mSign ? EqualIgnoreSign(*this, number) : false);
  290. }
  291. bool operator!=(BSNumber const& number) const
  292. {
  293. return !operator==(number);
  294. }
  295. bool operator< (BSNumber const& number) const
  296. {
  297. if (mSign > 0)
  298. {
  299. if (number.mSign <= 0)
  300. {
  301. return false;
  302. }
  303. // Both numbers are positive.
  304. return LessThanIgnoreSign(*this, number);
  305. }
  306. else if (mSign < 0)
  307. {
  308. if (number.mSign >= 0)
  309. {
  310. return true;
  311. }
  312. // Both numbers are negative.
  313. return LessThanIgnoreSign(number, *this);
  314. }
  315. else
  316. {
  317. return number.mSign > 0;
  318. }
  319. }
  320. bool operator<=(BSNumber const& number) const
  321. {
  322. return !number.operator<(*this);
  323. }
  324. bool operator> (BSNumber const& number) const
  325. {
  326. return number.operator<(*this);
  327. }
  328. bool operator>=(BSNumber const& number) const
  329. {
  330. return !operator<(number);
  331. }
  332. // Unary operations.
  333. BSNumber operator+() const
  334. {
  335. return *this;
  336. }
  337. BSNumber operator-() const
  338. {
  339. BSNumber result = *this;
  340. result.mSign = -result.mSign;
  341. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  342. result.mValue = (double)result;
  343. #endif
  344. return result;
  345. }
  346. // Arithmetic.
  347. BSNumber operator+(BSNumber const& n1) const
  348. {
  349. BSNumber const& n0 = *this;
  350. if (n0.mSign == 0)
  351. {
  352. return n1;
  353. }
  354. if (n1.mSign == 0)
  355. {
  356. return n0;
  357. }
  358. if (n0.mSign > 0)
  359. {
  360. if (n1.mSign > 0)
  361. {
  362. // n0 + n1 = |n0| + |n1|
  363. return AddIgnoreSign(n0, n1, +1);
  364. }
  365. else // n1.mSign < 0
  366. {
  367. if (!EqualIgnoreSign(n0, n1))
  368. {
  369. if (LessThanIgnoreSign(n1, n0))
  370. {
  371. // n0 + n1 = |n0| - |n1| > 0
  372. return SubIgnoreSign(n0, n1, +1);
  373. }
  374. else
  375. {
  376. // n0 + n1 = -(|n1| - |n0|) < 0
  377. return SubIgnoreSign(n1, n0, -1);
  378. }
  379. }
  380. // else n0 + n1 = 0
  381. }
  382. }
  383. else // n0.mSign < 0
  384. {
  385. if (n1.mSign < 0)
  386. {
  387. // n0 + n1 = -(|n0| + |n1|)
  388. return AddIgnoreSign(n0, n1, -1);
  389. }
  390. else // n1.mSign > 0
  391. {
  392. if (!EqualIgnoreSign(n0, n1))
  393. {
  394. if (LessThanIgnoreSign(n1, n0))
  395. {
  396. // n0 + n1 = -(|n0| - |n1|) < 0
  397. return SubIgnoreSign(n0, n1, -1);
  398. }
  399. else
  400. {
  401. // n0 + n1 = |n1| - |n0| > 0
  402. return SubIgnoreSign(n1, n0, +1);
  403. }
  404. }
  405. // else n0 + n1 = 0
  406. }
  407. }
  408. return BSNumber(); // = 0
  409. }
  410. BSNumber operator-(BSNumber const& n1) const
  411. {
  412. BSNumber const& n0 = *this;
  413. if (n0.mSign == 0)
  414. {
  415. return -n1;
  416. }
  417. if (n1.mSign == 0)
  418. {
  419. return n0;
  420. }
  421. if (n0.mSign > 0)
  422. {
  423. if (n1.mSign < 0)
  424. {
  425. // n0 - n1 = |n0| + |n1|
  426. return AddIgnoreSign(n0, n1, +1);
  427. }
  428. else // n1.mSign > 0
  429. {
  430. if (!EqualIgnoreSign(n0, n1))
  431. {
  432. if (LessThanIgnoreSign(n1, n0))
  433. {
  434. // n0 - n1 = |n0| - |n1| > 0
  435. return SubIgnoreSign(n0, n1, +1);
  436. }
  437. else
  438. {
  439. // n0 - n1 = -(|n1| - |n0|) < 0
  440. return SubIgnoreSign(n1, n0, -1);
  441. }
  442. }
  443. // else n0 - n1 = 0
  444. }
  445. }
  446. else // n0.mSign < 0
  447. {
  448. if (n1.mSign > 0)
  449. {
  450. // n0 - n1 = -(|n0| + |n1|)
  451. return AddIgnoreSign(n0, n1, -1);
  452. }
  453. else // n1.mSign < 0
  454. {
  455. if (!EqualIgnoreSign(n0, n1))
  456. {
  457. if (LessThanIgnoreSign(n1, n0))
  458. {
  459. // n0 - n1 = -(|n0| - |n1|) < 0
  460. return SubIgnoreSign(n0, n1, -1);
  461. }
  462. else
  463. {
  464. // n0 - n1 = |n1| - |n0| > 0
  465. return SubIgnoreSign(n1, n0, +1);
  466. }
  467. }
  468. // else n0 - n1 = 0
  469. }
  470. }
  471. return BSNumber(); // = 0
  472. }
  473. BSNumber operator*(BSNumber const& number) const
  474. {
  475. BSNumber result; // = 0
  476. int sign = mSign * number.mSign;
  477. if (sign != 0)
  478. {
  479. result.mSign = sign;
  480. result.mBiasedExponent = mBiasedExponent + number.mBiasedExponent;
  481. result.mUInteger.Mul(mUInteger, number.mUInteger);
  482. }
  483. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  484. result.mValue = (double)result;
  485. #endif
  486. return result;
  487. }
  488. BSNumber& operator+=(BSNumber const& number)
  489. {
  490. *this = operator+(number);
  491. return *this;
  492. }
  493. BSNumber& operator-=(BSNumber const& number)
  494. {
  495. *this = operator-(number);
  496. return *this;
  497. }
  498. BSNumber& operator*=(BSNumber const& number)
  499. {
  500. *this = operator*(number);
  501. return *this;
  502. }
  503. // Disk input/output. The fstream objects should be created using
  504. // std::ios::binary. The return value is 'true' iff the operation
  505. // was successful.
  506. bool Write(std::ostream& output) const
  507. {
  508. if (output.write((char const*)&mSign, sizeof(mSign)).bad())
  509. {
  510. return false;
  511. }
  512. if (output.write((char const*)&mBiasedExponent,
  513. sizeof(mBiasedExponent)).bad())
  514. {
  515. return false;
  516. }
  517. return mUInteger.Write(output);
  518. }
  519. bool Read(std::istream& input)
  520. {
  521. if (input.read((char*)&mSign, sizeof(mSign)).bad())
  522. {
  523. return false;
  524. }
  525. if (input.read((char*)&mBiasedExponent, sizeof(mBiasedExponent)).bad())
  526. {
  527. return false;
  528. }
  529. return mUInteger.Read(input);
  530. }
  531. private:
  532. // Helper for converting a string to a BSNumber. The string must be
  533. // valid for a nonnegative integer without a leading '+' sign.
  534. static BSNumber ConvertToInteger(std::string const& number)
  535. {
  536. int digit = static_cast<int>(number.back()) - static_cast<int>('0');
  537. BSNumber x(digit);
  538. if (number.size() > 1)
  539. {
  540. LogAssert(number.find_first_of("123456789") == 0, "Invalid number format.");
  541. LogAssert(number.find_first_not_of("0123456789") == std::string::npos, "Invalid number format.");
  542. BSNumber ten(10), pow10(10);
  543. for (size_t i = 1, j = number.size() - 2; i < number.size(); ++i, --j)
  544. {
  545. digit = static_cast<int>(number[j]) - static_cast<int>('0');
  546. if (digit > 0)
  547. {
  548. x += BSNumber(digit) * pow10;
  549. }
  550. pow10 *= ten;
  551. }
  552. }
  553. return x;
  554. }
  555. // Helpers for operator==, operator<, operator+, operator-.
  556. static bool EqualIgnoreSign(BSNumber const& n0, BSNumber const& n1)
  557. {
  558. return n0.mBiasedExponent == n1.mBiasedExponent && n0.mUInteger == n1.mUInteger;
  559. }
  560. static bool LessThanIgnoreSign(BSNumber const& n0, BSNumber const& n1)
  561. {
  562. int32_t e0 = n0.GetExponent(), e1 = n1.GetExponent();
  563. if (e0 < e1)
  564. {
  565. return true;
  566. }
  567. if (e0 > e1)
  568. {
  569. return false;
  570. }
  571. return n0.mUInteger < n1.mUInteger;
  572. }
  573. // Add two positive numbers.
  574. static BSNumber AddIgnoreSign(BSNumber const& n0, BSNumber const& n1, int32_t resultSign)
  575. {
  576. BSNumber result, temp;
  577. int32_t diff = n0.mBiasedExponent - n1.mBiasedExponent;
  578. if (diff > 0)
  579. {
  580. temp.mUInteger.ShiftLeft(n0.mUInteger, diff);
  581. result.mUInteger.Add(temp.mUInteger, n1.mUInteger);
  582. result.mBiasedExponent = n1.mBiasedExponent;
  583. }
  584. else if (diff < 0)
  585. {
  586. temp.mUInteger.ShiftLeft(n1.mUInteger, -diff);
  587. result.mUInteger.Add(n0.mUInteger, temp.mUInteger);
  588. result.mBiasedExponent = n0.mBiasedExponent;
  589. }
  590. else
  591. {
  592. temp.mUInteger.Add(n0.mUInteger, n1.mUInteger);
  593. int32_t shift = result.mUInteger.ShiftRightToOdd(temp.mUInteger);
  594. result.mBiasedExponent = n0.mBiasedExponent + shift;
  595. }
  596. result.mSign = resultSign;
  597. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  598. result.mValue = (double)result;
  599. #endif
  600. return result;
  601. }
  602. // Subtract two positive numbers where n0 > n1.
  603. static BSNumber SubIgnoreSign(BSNumber const& n0, BSNumber const& n1, int32_t resultSign)
  604. {
  605. BSNumber result, temp;
  606. int32_t diff = n0.mBiasedExponent - n1.mBiasedExponent;
  607. if (diff > 0)
  608. {
  609. temp.mUInteger.ShiftLeft(n0.mUInteger, diff);
  610. result.mUInteger.Sub(temp.mUInteger, n1.mUInteger);
  611. result.mBiasedExponent = n1.mBiasedExponent;
  612. }
  613. else if (diff < 0)
  614. {
  615. temp.mUInteger.ShiftLeft(n1.mUInteger, -diff);
  616. result.mUInteger.Sub(n0.mUInteger, temp.mUInteger);
  617. result.mBiasedExponent = n0.mBiasedExponent;
  618. }
  619. else
  620. {
  621. temp.mUInteger.Sub(n0.mUInteger, n1.mUInteger);
  622. int32_t shift = result.mUInteger.ShiftRightToOdd(temp.mUInteger);
  623. result.mBiasedExponent = n0.mBiasedExponent + shift;
  624. }
  625. result.mSign = resultSign;
  626. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  627. result.mValue = (double)result;
  628. #endif
  629. return result;
  630. }
  631. // Support for conversions from floating-point numbers to BSNumber.
  632. template <typename IEEE>
  633. void ConvertFrom(typename IEEE::FloatType number)
  634. {
  635. IEEE x(number);
  636. // Extract sign s, biased exponent e, and trailing significand t.
  637. typename IEEE::UIntType s = x.GetSign();
  638. typename IEEE::UIntType e = x.GetBiased();
  639. typename IEEE::UIntType t = x.GetTrailing();
  640. if (e == 0)
  641. {
  642. if (t == 0) // zeros
  643. {
  644. // x = (-1)^s * 0
  645. mSign = 0;
  646. mBiasedExponent = 0;
  647. }
  648. else // subnormal numbers
  649. {
  650. // x = (-1)^s * 0.t * 2^{1-EXPONENT_BIAS}
  651. int32_t last = BitHacks::GetTrailingBit(t);
  652. int32_t diff = IEEE::NUM_TRAILING_BITS - last;
  653. mSign = (s > 0 ? -1 : 1);
  654. mBiasedExponent = IEEE::MIN_SUB_EXPONENT - diff;
  655. mUInteger = (t >> last);
  656. }
  657. }
  658. else if (e < IEEE::MAX_BIASED_EXPONENT) // normal numbers
  659. {
  660. // x = (-1)^s * 1.t * 2^{e-EXPONENT_BIAS}
  661. if (t > 0)
  662. {
  663. int32_t last = BitHacks::GetTrailingBit(t);
  664. int32_t diff = IEEE::NUM_TRAILING_BITS - last;
  665. mSign = (s > 0 ? -1 : 1);
  666. mBiasedExponent =
  667. static_cast<int32_t>(e) - IEEE::EXPONENT_BIAS - diff;
  668. mUInteger = ((t | IEEE::SUP_TRAILING) >> last);
  669. }
  670. else
  671. {
  672. mSign = (s > 0 ? -1 : 1);
  673. mBiasedExponent = static_cast<int32_t>(e) - IEEE::EXPONENT_BIAS;
  674. mUInteger = (typename IEEE::UIntType)1;
  675. }
  676. }
  677. else // e == MAX_BIASED_EXPONENT, special numbers
  678. {
  679. if (t == 0) // infinities
  680. {
  681. // x = (-1)^s * infinity
  682. #if defined(GTE_THROW_ON_BSNUMBER_ERRORS)
  683. // TODO: This should not be a warning, but BSNumber does
  684. // not have a representation for +infinity or -infinity.
  685. // Consider doing so with mSign in {-2,2} and all other
  686. // members zero-valued.
  687. LogWarning("Input is " + std::string(s > 0 ? "-" : "+") + "infinity.");
  688. #else
  689. // Return (-1)^s * 2^{1+EXPONENT_BIAS} for a graceful
  690. // exit.
  691. mSign = (s > 0 ? -1 : 1);
  692. mBiasedExponent = 1 + IEEE::EXPONENT_BIAS;
  693. mUInteger = (typename IEEE::UIntType)1;
  694. #endif
  695. }
  696. else // not-a-number (NaN)
  697. {
  698. #if defined(GTE_THROW_ON_BSNUMBER_ERRORS)
  699. // TODO: BSNumber does not have a representation for
  700. // NaNs. Consider doing so with mSign in {-3,3} and a
  701. // payload stored in mBits.
  702. LogError("Input is a " +
  703. std::string(t & IEEE::NAN_QUIET_MASK ?
  704. "quiet" : "signaling") + " NaN with payload " +
  705. std::to_string(t & IEEE::NAN_PAYLOAD_MASK) + ".");
  706. #else
  707. // Return 0 for a graceful exit.
  708. mSign = 0;
  709. mBiasedExponent = 0;
  710. #endif
  711. }
  712. }
  713. }
  714. // Support for conversions from BSNumber to floating-point numbers.
  715. template <typename IEEE>
  716. typename IEEE::FloatType ConvertTo() const
  717. {
  718. typename IEEE::UIntType s = (mSign < 0 ? 1 : 0);
  719. typename IEEE::UIntType e, t;
  720. if (mSign != 0)
  721. {
  722. // The conversions use round-to-nearest-ties-to-even
  723. // semantics.
  724. int32_t exponent = GetExponent();
  725. if (exponent < IEEE::MIN_EXPONENT)
  726. {
  727. if (exponent < IEEE::MIN_EXPONENT - 1
  728. || mUInteger.GetNumBits() == 1)
  729. {
  730. // x = 1.0*2^{MIN_EXPONENT-1}, round to zero.
  731. e = 0;
  732. t = 0;
  733. }
  734. else
  735. {
  736. // Round to min subnormal.
  737. e = 0;
  738. t = 1;
  739. }
  740. }
  741. else if (exponent < IEEE::MIN_SUB_EXPONENT)
  742. {
  743. // The second input is in {0, ..., NUM_TRAILING_BITS-1}.
  744. t = GetTrailing<IEEE>(0, IEEE::MIN_SUB_EXPONENT - exponent - 1);
  745. if (t & IEEE::SUP_TRAILING)
  746. {
  747. // Leading NUM_SIGNIFICAND_BITS bits were all 1, so
  748. // round to min normal.
  749. e = 1;
  750. t = 0;
  751. }
  752. else
  753. {
  754. e = 0;
  755. }
  756. }
  757. else if (exponent <= IEEE::EXPONENT_BIAS)
  758. {
  759. e = static_cast<uint32_t>(exponent + IEEE::EXPONENT_BIAS);
  760. t = GetTrailing<IEEE>(1, 0);
  761. if (t & (IEEE::SUP_TRAILING << 1))
  762. {
  763. // Carry-out occurred, so increase exponent by 1 and
  764. // shift right to compensate.
  765. ++e;
  766. t >>= 1;
  767. }
  768. // Eliminate the leading 1 (implied for normals).
  769. t &= ~IEEE::SUP_TRAILING;
  770. }
  771. else
  772. {
  773. // Set to infinity.
  774. e = IEEE::MAX_BIASED_EXPONENT;
  775. t = 0;
  776. }
  777. }
  778. else
  779. {
  780. // The input is zero.
  781. e = 0;
  782. t = 0;
  783. }
  784. IEEE x(s, e, t);
  785. return x.number;
  786. }
  787. template <typename IEEE>
  788. typename IEEE::UIntType GetTrailing(int32_t normal, int32_t sigma) const
  789. {
  790. int32_t const numRequested = IEEE::NUM_SIGNIFICAND_BITS + normal;
  791. // We need numRequested bits to determine rounding direction.
  792. // These are stored in the high-order bits of 'prefix'.
  793. uint64_t prefix = mUInteger.GetPrefix(numRequested);
  794. // The first bit index after the implied binary point for
  795. // rounding.
  796. int32_t diff = numRequested - sigma;
  797. int32_t roundBitIndex = 64 - diff;
  798. // Determine value based on round-to-nearest-ties-to-even.
  799. uint64_t mask = (1ull << roundBitIndex);
  800. uint64_t round;
  801. if (prefix & mask)
  802. {
  803. // The first bit of the remainder is 1.
  804. if (mUInteger.GetNumBits() == diff)
  805. {
  806. // The first bit of the remainder is the lowest-order bit
  807. // of mBits[0]. Apply the ties-to-even rule.
  808. if (prefix & (mask << 1))
  809. {
  810. // The last bit of the trailing significand is odd,
  811. // so round up.
  812. round = 1;
  813. }
  814. else
  815. {
  816. // The last bit of the trailing significand is even,
  817. // so round down.
  818. round = 0;
  819. }
  820. }
  821. else
  822. {
  823. // The first bit of the remainder is not the lowest-order
  824. // bit of mBits[0]. The remainder as a fraction is larger
  825. // than 1/2, so round up.
  826. round = 1;
  827. }
  828. }
  829. else
  830. {
  831. // The first bit of the remainder is 0, so round down.
  832. round = 0;
  833. }
  834. // Get the unrounded trailing significand.
  835. uint64_t trailing = prefix >> (roundBitIndex + 1);
  836. // Apply the rounding.
  837. trailing += round;
  838. return static_cast<typename IEEE::UIntType>(trailing);
  839. }
  840. #if defined(GTE_BINARY_SCIENTIFIC_SHOW_DOUBLE)
  841. public:
  842. // List this first so that it shows up first in the debugger watch
  843. // window.
  844. double mValue;
  845. private:
  846. #endif
  847. // The number 0 is represented by: mSign = 0, mBiasedExponent = 0 and
  848. // mUInteger = 0. For nonzero numbers, mSign != 0 and mUInteger > 0.
  849. int32_t mSign;
  850. int32_t mBiasedExponent;
  851. UInteger mUInteger;
  852. // BSRational depends on the design of BSNumber, so allow it to have
  853. // full access to the implementation.
  854. friend class BSRational<UInteger>;
  855. };
  856. // Explicit conversion to a user-specified precision. The rounding
  857. // mode is one of the flags provided in <cfenv>. The modes are
  858. // FE_TONEAREST: round to nearest ties to even
  859. // FE_DOWNWARD: round towards negative infinity
  860. // FE_TOWARDZERO: round towards zero
  861. // FE_UPWARD: round towards positive infinity
  862. template <typename UInteger>
  863. void Convert(BSNumber<UInteger> const& input, int32_t precision,
  864. int32_t roundingMode, BSNumber<UInteger>& output)
  865. {
  866. if (precision <= 0)
  867. {
  868. LogError("Precision must be positive.");
  869. }
  870. int64_t const maxSize = static_cast<int64_t>(UInteger::GetMaxSize());
  871. int64_t const excess = 32LL * maxSize - static_cast<int64_t>(precision);
  872. if (excess <= 0)
  873. {
  874. LogError("The maximum precision has been exceeded.");
  875. }
  876. if (input.GetSign() == 0)
  877. {
  878. output = BSNumber<UInteger>(0);
  879. return;
  880. }
  881. // Let p = precision and n+1 be the number of bits of the input.
  882. // Compute n+1-p. If it is nonpositive, then the requested precision
  883. // is already satisfied by the input.
  884. int32_t np1mp = input.GetUInteger().GetNumBits() - precision;
  885. if (np1mp <= 0)
  886. {
  887. output = input;
  888. return;
  889. }
  890. // At this point, the requested number of bits is smaller than the
  891. // number of bits in the input. Round the input to the smaller number
  892. // of bits using the specified rounding mode.
  893. UInteger& outW = output.GetUInteger();
  894. outW.SetNumBits(precision);
  895. outW.SetAllBitsToZero();
  896. int32_t const outSize = outW.GetSize();
  897. int32_t const precisionM1 = precision - 1;
  898. int32_t const outLeading = precisionM1 % 32;
  899. uint32_t outMask = (1 << outLeading);
  900. auto& outBits = outW.GetBits();
  901. int32_t outCurrent = outSize - 1;
  902. UInteger const& inW = input.GetUInteger();
  903. int32_t const inSize = inW.GetSize();
  904. int32_t const inLeading = (inW.GetNumBits() - 1) % 32;
  905. uint32_t inMask = (1 << inLeading);
  906. auto const& inBits = inW.GetBits();
  907. int32_t inCurrent = inSize - 1;
  908. int32_t lastBit = -1;
  909. for (int i = precisionM1; i >= 0; --i)
  910. {
  911. if (inBits[inCurrent] & inMask)
  912. {
  913. outBits[outCurrent] |= outMask;
  914. lastBit = 1;
  915. }
  916. else
  917. {
  918. lastBit = 0;
  919. }
  920. if (inMask == 0x00000001u)
  921. {
  922. --inCurrent;
  923. inMask = 0x80000000u;
  924. }
  925. else
  926. {
  927. inMask >>= 1;
  928. }
  929. if (outMask == 0x00000001u)
  930. {
  931. --outCurrent;
  932. outMask = 0x80000000u;
  933. }
  934. else
  935. {
  936. outMask >>= 1;
  937. }
  938. }
  939. // At this point as a sequence of bits, r = u_{n-p} ... u_0.
  940. int32_t sign = input.GetSign();
  941. int32_t outExponent = input.GetExponent();
  942. if (roundingMode == FE_TONEAREST)
  943. {
  944. // Determine whether u_{n-p} is positive.
  945. uint32_t positive = (inBits[inCurrent] & inMask) != 0u;
  946. if (positive && (np1mp > 1 || lastBit == 1))
  947. {
  948. // round up
  949. outExponent += outW.RoundUp();
  950. }
  951. // else round down, equivalent to truncating the r bits
  952. }
  953. else if (roundingMode == FE_UPWARD)
  954. {
  955. // The remainder r must be positive because n-p >= 0 and u_0 = 1.
  956. if (sign > 0)
  957. {
  958. // round up
  959. outExponent += outW.RoundUp();
  960. }
  961. // else round down, equivalent to truncating the r bits
  962. }
  963. else if (roundingMode == FE_DOWNWARD)
  964. {
  965. // The remainder r must be positive because n-p >= 0 and u_0 = 1.
  966. if (sign < 0)
  967. {
  968. // Round down. This is the round-up operation applied to
  969. // w, but the final sign is negative which amounts to
  970. // rounding down.
  971. outExponent += outW.RoundUp();
  972. }
  973. // else round down, equivalent to truncating the r bits
  974. }
  975. else if (roundingMode != FE_TOWARDZERO)
  976. {
  977. // Currently, no additional implementation-dependent modes
  978. // are supported for rounding.
  979. LogError("Implementation-dependent rounding mode not supported.");
  980. }
  981. // else roundingMode == FE_TOWARDZERO. Truncate the r bits, which
  982. // requires no additional work.
  983. output.SetSign(sign);
  984. output.SetBiasedExponent(outExponent - precisionM1);
  985. }
  986. }
  987. namespace std
  988. {
  989. // TODO: Allow for implementations of the math functions in which a
  990. // specified precision is used when computing the result.
  991. template <typename UInteger>
  992. inline WwiseGTE::BSNumber<UInteger> acos(WwiseGTE::BSNumber<UInteger> const& x)
  993. {
  994. return (WwiseGTE::BSNumber<UInteger>)std::acos((double)x);
  995. }
  996. template <typename UInteger>
  997. inline WwiseGTE::BSNumber<UInteger> acosh(WwiseGTE::BSNumber<UInteger> const& x)
  998. {
  999. return (WwiseGTE::BSNumber<UInteger>)std::acosh((double)x);
  1000. }
  1001. template <typename UInteger>
  1002. inline WwiseGTE::BSNumber<UInteger> asin(WwiseGTE::BSNumber<UInteger> const& x)
  1003. {
  1004. return (WwiseGTE::BSNumber<UInteger>)std::asin((double)x);
  1005. }
  1006. template <typename UInteger>
  1007. inline WwiseGTE::BSNumber<UInteger> asinh(WwiseGTE::BSNumber<UInteger> const& x)
  1008. {
  1009. return (WwiseGTE::BSNumber<UInteger>)std::asinh((double)x);
  1010. }
  1011. template <typename UInteger>
  1012. inline WwiseGTE::BSNumber<UInteger> atan(WwiseGTE::BSNumber<UInteger> const& x)
  1013. {
  1014. return (WwiseGTE::BSNumber<UInteger>)std::atan((double)x);
  1015. }
  1016. template <typename UInteger>
  1017. inline WwiseGTE::BSNumber<UInteger> atanh(WwiseGTE::BSNumber<UInteger> const& x)
  1018. {
  1019. return (WwiseGTE::BSNumber<UInteger>)std::atanh((double)x);
  1020. }
  1021. template <typename UInteger>
  1022. inline WwiseGTE::BSNumber<UInteger> atan2(WwiseGTE::BSNumber<UInteger> const& y, WwiseGTE::BSNumber<UInteger> const& x)
  1023. {
  1024. return (WwiseGTE::BSNumber<UInteger>)std::atan2((double)y, (double)x);
  1025. }
  1026. template <typename UInteger>
  1027. inline WwiseGTE::BSNumber<UInteger> ceil(WwiseGTE::BSNumber<UInteger> const& x)
  1028. {
  1029. return (WwiseGTE::BSNumber<UInteger>)std::ceil((double)x);
  1030. }
  1031. template <typename UInteger>
  1032. inline WwiseGTE::BSNumber<UInteger> cos(WwiseGTE::BSNumber<UInteger> const& x)
  1033. {
  1034. return (WwiseGTE::BSNumber<UInteger>)std::cos((double)x);
  1035. }
  1036. template <typename UInteger>
  1037. inline WwiseGTE::BSNumber<UInteger> cosh(WwiseGTE::BSNumber<UInteger> const& x)
  1038. {
  1039. return (WwiseGTE::BSNumber<UInteger>)std::cosh((double)x);
  1040. }
  1041. template <typename UInteger>
  1042. inline WwiseGTE::BSNumber<UInteger> exp(WwiseGTE::BSNumber<UInteger> const& x)
  1043. {
  1044. return (WwiseGTE::BSNumber<UInteger>)std::exp((double)x);
  1045. }
  1046. template <typename UInteger>
  1047. inline WwiseGTE::BSNumber<UInteger> exp2(WwiseGTE::BSNumber<UInteger> const& x)
  1048. {
  1049. return (WwiseGTE::BSNumber<UInteger>)std::exp2((double)x);
  1050. }
  1051. template <typename UInteger>
  1052. inline WwiseGTE::BSNumber<UInteger> fabs(WwiseGTE::BSNumber<UInteger> const& x)
  1053. {
  1054. return (x.GetSign() >= 0 ? x : -x);
  1055. }
  1056. template <typename UInteger>
  1057. inline WwiseGTE::BSNumber<UInteger> floor(WwiseGTE::BSNumber<UInteger> const& x)
  1058. {
  1059. return (WwiseGTE::BSNumber<UInteger>)std::floor((double)x);
  1060. }
  1061. template <typename UInteger>
  1062. inline WwiseGTE::BSNumber<UInteger> fmod(WwiseGTE::BSNumber<UInteger> const& x, WwiseGTE::BSNumber<UInteger> const& y)
  1063. {
  1064. return (WwiseGTE::BSNumber<UInteger>)std::fmod((double)x, (double)y);
  1065. }
  1066. template <typename UInteger>
  1067. inline WwiseGTE::BSNumber<UInteger> frexp(WwiseGTE::BSNumber<UInteger> const& x, int* exponent)
  1068. {
  1069. if (x.GetSign() != 0)
  1070. {
  1071. WwiseGTE::BSNumber<UInteger> result = x;
  1072. *exponent = result.GetExponent() + 1;
  1073. result.SetExponent(-1);
  1074. return result;
  1075. }
  1076. else
  1077. {
  1078. *exponent = 0;
  1079. return WwiseGTE::BSNumber<UInteger>(0);
  1080. }
  1081. }
  1082. template <typename UInteger>
  1083. inline WwiseGTE::BSNumber<UInteger> ldexp(WwiseGTE::BSNumber<UInteger> const& x, int exponent)
  1084. {
  1085. WwiseGTE::BSNumber<UInteger> result = x;
  1086. result.SetBiasedExponent(result.GetBiasedExponent() + exponent);
  1087. return result;
  1088. }
  1089. template <typename UInteger>
  1090. inline WwiseGTE::BSNumber<UInteger> log(WwiseGTE::BSNumber<UInteger> const& x)
  1091. {
  1092. return (WwiseGTE::BSNumber<UInteger>)std::log((double)x);
  1093. }
  1094. template <typename UInteger>
  1095. inline WwiseGTE::BSNumber<UInteger> log2(WwiseGTE::BSNumber<UInteger> const& x)
  1096. {
  1097. return (WwiseGTE::BSNumber<UInteger>)std::log2((double)x);
  1098. }
  1099. template <typename UInteger>
  1100. inline WwiseGTE::BSNumber<UInteger> log10(WwiseGTE::BSNumber<UInteger> const& x)
  1101. {
  1102. return (WwiseGTE::BSNumber<UInteger>)std::log10((double)x);
  1103. }
  1104. template <typename UInteger>
  1105. inline WwiseGTE::BSNumber<UInteger> pow(WwiseGTE::BSNumber<UInteger> const& x, WwiseGTE::BSNumber<UInteger> const& y)
  1106. {
  1107. return (WwiseGTE::BSNumber<UInteger>)std::pow((double)x, (double)y);
  1108. }
  1109. template <typename UInteger>
  1110. inline WwiseGTE::BSNumber<UInteger> sin(WwiseGTE::BSNumber<UInteger> const& x)
  1111. {
  1112. return (WwiseGTE::BSNumber<UInteger>)std::sin((double)x);
  1113. }
  1114. template <typename UInteger>
  1115. inline WwiseGTE::BSNumber<UInteger> sinh(WwiseGTE::BSNumber<UInteger> const& x)
  1116. {
  1117. return (WwiseGTE::BSNumber<UInteger>)std::sinh((double)x);
  1118. }
  1119. template <typename UInteger>
  1120. inline WwiseGTE::BSNumber<UInteger> sqrt(WwiseGTE::BSNumber<UInteger> const& x)
  1121. {
  1122. return (WwiseGTE::BSNumber<UInteger>)std::sqrt((double)x);
  1123. }
  1124. template <typename UInteger>
  1125. inline WwiseGTE::BSNumber<UInteger> tan(WwiseGTE::BSNumber<UInteger> const& x)
  1126. {
  1127. return (WwiseGTE::BSNumber<UInteger>)std::tan((double)x);
  1128. }
  1129. template <typename UInteger>
  1130. inline WwiseGTE::BSNumber<UInteger> tanh(WwiseGTE::BSNumber<UInteger> const& x)
  1131. {
  1132. return (WwiseGTE::BSNumber<UInteger>)std::tanh((double)x);
  1133. }
  1134. // Type trait that says BSNumber is a signed type.
  1135. template <typename UInteger>
  1136. struct is_signed<WwiseGTE::BSNumber<UInteger>> : true_type {};
  1137. }
  1138. namespace WwiseGTE
  1139. {
  1140. template <typename UInteger>
  1141. inline BSNumber<UInteger> atandivpi(BSNumber<UInteger> const& x)
  1142. {
  1143. return (BSNumber<UInteger>)atandivpi((double)x);
  1144. }
  1145. template <typename UInteger>
  1146. inline BSNumber<UInteger> atan2divpi(BSNumber<UInteger> const& y, BSNumber<UInteger> const& x)
  1147. {
  1148. return (BSNumber<UInteger>)atan2divpi((double)y, (double)x);
  1149. }
  1150. template <typename UInteger>
  1151. inline BSNumber<UInteger> clamp(BSNumber<UInteger> const& x, BSNumber<UInteger> const& xmin, BSNumber<UInteger> const& xmax)
  1152. {
  1153. return (BSNumber<UInteger>)clamp((double)x, (double)xmin, (double)xmax);
  1154. }
  1155. template <typename UInteger>
  1156. inline BSNumber<UInteger> cospi(BSNumber<UInteger> const& x)
  1157. {
  1158. return (BSNumber<UInteger>)cospi((double)x);
  1159. }
  1160. template <typename UInteger>
  1161. inline BSNumber<UInteger> exp10(BSNumber<UInteger> const& x)
  1162. {
  1163. return (BSNumber<UInteger>)exp10((double)x);
  1164. }
  1165. template <typename UInteger>
  1166. inline BSNumber<UInteger> invsqrt(BSNumber<UInteger> const& x)
  1167. {
  1168. return (BSNumber<UInteger>)invsqrt((double)x);
  1169. }
  1170. template <typename UInteger>
  1171. inline int isign(BSNumber<UInteger> const& x)
  1172. {
  1173. return isign((double)x);
  1174. }
  1175. template <typename UInteger>
  1176. inline BSNumber<UInteger> saturate(BSNumber<UInteger> const& x)
  1177. {
  1178. return (BSNumber<UInteger>)saturate((double)x);
  1179. }
  1180. template <typename UInteger>
  1181. inline BSNumber<UInteger> sign(BSNumber<UInteger> const& x)
  1182. {
  1183. return (BSNumber<UInteger>)sign((double)x);
  1184. }
  1185. template <typename UInteger>
  1186. inline BSNumber<UInteger> sinpi(BSNumber<UInteger> const& x)
  1187. {
  1188. return (BSNumber<UInteger>)sinpi((double)x);
  1189. }
  1190. template <typename UInteger>
  1191. inline BSNumber<UInteger> sqr(BSNumber<UInteger> const& x)
  1192. {
  1193. return (BSNumber<UInteger>)sqr((double)x);
  1194. }
  1195. // See the comments in GteMath.h about trait is_arbitrary_precision.
  1196. template <typename UInteger>
  1197. struct is_arbitrary_precision_internal<BSNumber<UInteger>> : std::true_type {};
  1198. }