dragonbox.h 116 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285
  1. // Copyright 2020 Junekey Jeon
  2. //
  3. // The contents of this file may be used under the terms of
  4. // the Apache License v2.0 with LLVM Exceptions.
  5. //
  6. // (See accompanying file LICENSE-Apache or copy at
  7. // https://llvm.org/foundation/relicensing/LICENSE.txt)
  8. //
  9. // Alternatively, the contents of this file may be used under the terms of
  10. // the Boost Software License, Version 1.0.
  11. // (See accompanying file LICENSE-Boost or copy at
  12. // https://www.boost.org/LICENSE_1_0.txt)
  13. //
  14. // Unless required by applicable law or agreed to in writing, this software
  15. // is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
  16. // KIND, either express or implied.
  17. #ifndef JKJ_DRAGONBOX
  18. #define JKJ_DRAGONBOX
  19. #include <cassert>
  20. #include <cstdint>
  21. #include <cstring>
  22. #include <limits>
  23. #include <type_traits>
  24. // Suppress additional buffer overrun check
  25. // I have no idea why MSVC thinks some functions here are vulnerable to the buffer overrun attacks
  26. // No, they aren't.
  27. #if defined(__GNUC__) || defined(__clang__)
  28. #define JKJ_SAFEBUFFERS
  29. #define JKJ_FORCEINLINE inline __attribute__((always_inline))
  30. #elif defined(_MSC_VER)
  31. #define JKJ_SAFEBUFFERS __declspec(safebuffers)
  32. #define JKJ_FORCEINLINE __forceinline
  33. #else
  34. #define JKJ_SAFEBUFFERS
  35. #define JKJ_FORCEINLINE inline
  36. #endif
  37. #if defined(_MSC_VER)
  38. #include <intrin.h> // this includes immintrin.h as well
  39. #elif (defined(__GNUC__) || defined(__clang__)) && defined(__x86_64__)
  40. #include <immintrin.h>
  41. #endif
  42. namespace jkj::dragonbox {
  43. namespace detail {
  44. template <class T>
  45. constexpr std::size_t physical_bits = sizeof(T) * std::numeric_limits<unsigned char>::digits;
  46. template <class T>
  47. constexpr std::size_t value_bits =
  48. std::numeric_limits<std::enable_if_t<std::is_unsigned_v<T>, T>>::digits;
  49. }
  50. enum class ieee754_format {
  51. binary32,
  52. binary64
  53. };
  54. template <ieee754_format format_>
  55. struct ieee754_format_info;
  56. template <>
  57. struct ieee754_format_info<ieee754_format::binary32> {
  58. static constexpr auto format = ieee754_format::binary32;
  59. static constexpr int significand_bits = 23;
  60. static constexpr int exponent_bits = 8;
  61. static constexpr int min_exponent = -126;
  62. static constexpr int max_exponent = 127;
  63. static constexpr int exponent_bias = -127;
  64. static constexpr int decimal_digits = 9;
  65. };
  66. template <>
  67. struct ieee754_format_info<ieee754_format::binary64> {
  68. static constexpr auto format = ieee754_format::binary64;
  69. static constexpr int significand_bits = 52;
  70. static constexpr int exponent_bits = 11;
  71. static constexpr int min_exponent = -1022;
  72. static constexpr int max_exponent = 1023;
  73. static constexpr int exponent_bias = -1023;
  74. static constexpr int decimal_digits = 17;
  75. };
  76. // To reduce boilerplates
  77. template <class T>
  78. struct default_ieee754_traits {
  79. static_assert(detail::physical_bits<T> == 32 || detail::physical_bits<T> == 64);
  80. using type = T;
  81. static constexpr ieee754_format format =
  82. detail::physical_bits<T> == 32 ? ieee754_format::binary32 : ieee754_format::binary64;
  83. using carrier_uint = std::conditional_t<
  84. detail::physical_bits<T> == 32,
  85. std::uint32_t,
  86. std::uint64_t>;
  87. static_assert(sizeof(carrier_uint) == sizeof(T));
  88. static constexpr int carrier_bits = int(detail::physical_bits<carrier_uint>);
  89. static T carrier_to_float(carrier_uint u) noexcept {
  90. T x;
  91. std::memcpy(&x, &u, sizeof(carrier_uint));
  92. return x;
  93. }
  94. static carrier_uint float_to_carrier(T x) noexcept {
  95. carrier_uint u;
  96. std::memcpy(&u, &x, sizeof(carrier_uint));
  97. return u;
  98. }
  99. static constexpr unsigned int extract_exponent_bits(carrier_uint u) noexcept {
  100. constexpr int significand_bits = ieee754_format_info<format>::significand_bits;
  101. constexpr int exponent_bits = ieee754_format_info<format>::exponent_bits;
  102. static_assert(detail::value_bits<unsigned int> > exponent_bits);
  103. constexpr auto exponent_bits_mask = (unsigned int)(((unsigned int)(1) << exponent_bits) - 1);
  104. return (unsigned int)((u >> significand_bits) & exponent_bits_mask);
  105. }
  106. static constexpr carrier_uint extract_significand_bits(carrier_uint u) noexcept {
  107. constexpr int significand_bits = ieee754_format_info<format>::significand_bits;
  108. constexpr auto significand_bits_mask = carrier_uint((carrier_uint(1) << significand_bits) - 1);
  109. return carrier_uint(u & significand_bits_mask);
  110. }
  111. // Allows positive zero and positive NaN's, but not allow negative zero
  112. static constexpr bool is_positive(carrier_uint u) noexcept {
  113. return (u >> (carrier_bits - 1)) == 0;
  114. }
  115. // Allows negative zero and negative NaN's, but not allow positive zero
  116. static constexpr bool is_negative(carrier_uint u) noexcept {
  117. return (u >> (carrier_bits - 1)) != 0;
  118. }
  119. static constexpr int exponent_bias = 1 - (1 << (carrier_bits - ieee754_format_info<format>::significand_bits - 2));
  120. static constexpr bool is_finite(carrier_uint u) noexcept {
  121. constexpr int significand_bits = ieee754_format_info<format>::significand_bits;
  122. constexpr int exponent_bits = ieee754_format_info<format>::exponent_bits;
  123. constexpr auto exponent_bits_mask =
  124. carrier_uint(((carrier_uint(1) << exponent_bits) - 1) << significand_bits);
  125. return (u & exponent_bits_mask) != exponent_bits_mask;
  126. }
  127. static constexpr bool is_nonzero(carrier_uint u) noexcept {
  128. return (u << 1) != 0;
  129. }
  130. // Allows positive and negative zeros
  131. static constexpr bool is_subnormal(carrier_uint u) noexcept {
  132. constexpr int significand_bits = ieee754_format_info<format>::significand_bits;
  133. constexpr int exponent_bits = ieee754_format_info<format>::exponent_bits;
  134. constexpr auto exponent_bits_mask =
  135. carrier_uint(((carrier_uint(1) << exponent_bits) - 1) << significand_bits);
  136. return (u & exponent_bits_mask) == 0;
  137. }
  138. static constexpr bool is_positive_infinity(carrier_uint u) noexcept {
  139. constexpr int significand_bits = ieee754_format_info<format>::significand_bits;
  140. constexpr int exponent_bits = ieee754_format_info<format>::exponent_bits;
  141. constexpr auto positive_infinity =
  142. carrier_uint((carrier_uint(1) << exponent_bits) - 1) << significand_bits;
  143. return u == positive_infinity;
  144. }
  145. static constexpr bool is_negative_infinity(carrier_uint u) noexcept {
  146. constexpr int significand_bits = ieee754_format_info<format>::significand_bits;
  147. constexpr int exponent_bits = ieee754_format_info<format>::exponent_bits;
  148. constexpr auto negative_infinity =
  149. (carrier_uint((carrier_uint(1) << exponent_bits) - 1) << significand_bits)
  150. | (carrier_uint(1) << (carrier_bits - 1));
  151. return u == negative_infinity;
  152. }
  153. static constexpr bool is_infinity(carrier_uint u) noexcept {
  154. return is_positive_infinity(u) || is_negative_infinity(u);
  155. }
  156. static constexpr bool is_nan(carrier_uint u) noexcept {
  157. return !is_finite(u) && (extract_significand_bits(u) != 0);
  158. }
  159. };
  160. // Speciailze this class template for possible extensions
  161. template <class T>
  162. struct ieee754_traits : default_ieee754_traits<T> {
  163. // I don't know if there is a truly reliable way of detecting
  164. // IEEE-754 binary32/binary64 formats; I just did my best here
  165. static_assert(std::numeric_limits<T>::is_iec559 &&
  166. std::numeric_limits<T>::radix == 2 &&
  167. (detail::physical_bits<T> == 32 || detail::physical_bits<T> == 64),
  168. "default_ieee754_traits only worsk for 32-bits or 64-bits types "
  169. "supporting binary32 or binary64 formats!");
  170. };
  171. // Convenient wrapper for ieee754_traits
  172. // In order to reduce the argument passing overhead,
  173. // this class should be as simple as possible
  174. // (e.g., no inheritance, no private non-static data member, etc.;
  175. // this is an unfortunate fact about x64 calling convention)
  176. template <class T>
  177. struct ieee754_bits {
  178. using carrier_uint = typename ieee754_traits<T>::carrier_uint;
  179. carrier_uint u;
  180. ieee754_bits() = default;
  181. constexpr explicit ieee754_bits(carrier_uint u) noexcept : u{ u } {}
  182. constexpr explicit ieee754_bits(T x) noexcept : u{ ieee754_traits<T>::float_to_carrier(x) } {}
  183. constexpr T to_float() const noexcept {
  184. return ieee754_traits<T>::carrier_to_float(u);
  185. }
  186. constexpr carrier_uint extract_significand_bits() const noexcept {
  187. return ieee754_traits<T>::extract_significand_bits(u);
  188. }
  189. constexpr unsigned int extract_exponent_bits() const noexcept {
  190. return ieee754_traits<T>::extract_exponent_bits(u);
  191. }
  192. constexpr carrier_uint binary_significand() const noexcept {
  193. using format_info = ieee754_format_info<ieee754_traits<T>::format>;
  194. auto s = extract_significand_bits();
  195. if (extract_exponent_bits() == 0) {
  196. return s;
  197. }
  198. else {
  199. return s | (carrier_uint(1) << format_info::significand_bits);
  200. }
  201. }
  202. constexpr int binary_exponent() const noexcept {
  203. using format_info = ieee754_format_info<ieee754_traits<T>::format>;
  204. auto e = extract_exponent_bits();
  205. if (e == 0) {
  206. return format_info::min_exponent;
  207. }
  208. else {
  209. return e + format_info::exponent_bias;
  210. }
  211. }
  212. constexpr bool is_finite() const noexcept {
  213. return ieee754_traits<T>::is_finite(u);
  214. }
  215. constexpr bool is_nonzero() const noexcept {
  216. return ieee754_traits<T>::is_nonzero(u);
  217. }
  218. // Allows positive and negative zeros
  219. constexpr bool is_subnormal() const noexcept {
  220. return ieee754_traits<T>::is_subnormal(u);
  221. }
  222. // Allows positive zero and positive NaN's, but not allow negative zero
  223. constexpr bool is_positive() const noexcept {
  224. return ieee754_traits<T>::is_positive(u);
  225. }
  226. // Allows negative zero and negative NaN's, but not allow positive zero
  227. constexpr bool is_negative() const noexcept {
  228. return ieee754_traits<T>::is_negative(u);
  229. }
  230. constexpr bool is_positive_infinity() const noexcept {
  231. return ieee754_traits<T>::is_positive_infinity(u);
  232. }
  233. constexpr bool is_negative_infinity() const noexcept {
  234. return ieee754_traits<T>::is_negative_infinity(u);
  235. }
  236. // Allows both plus and minus infinities
  237. constexpr bool is_infinity() const noexcept {
  238. return ieee754_traits<T>::is_infinity(u);
  239. }
  240. constexpr bool is_nan() const noexcept {
  241. return ieee754_traits<T>::is_nan(u);
  242. }
  243. };
  244. namespace detail {
  245. ////////////////////////////////////////////////////////////////////////////////////////
  246. // Bit operation intrinsics
  247. ////////////////////////////////////////////////////////////////////////////////////////
  248. namespace bits {
  249. template <class UInt>
  250. inline int countr_zero(UInt n) noexcept {
  251. static_assert(std::is_unsigned_v<UInt> && value_bits<UInt> <= 64);
  252. #if (defined(__GNUC__) || defined(__clang__)) && defined(__x86_64__)
  253. #define JKJ_HAS_COUNTR_ZERO_INTRINSIC 1
  254. if constexpr (std::is_same_v<UInt, unsigned long>) {
  255. return __builtin_ctzl(n);
  256. }
  257. else if constexpr (std::is_same_v<UInt, unsigned long long>) {
  258. return __builtin_ctzll(n);
  259. }
  260. else {
  261. static_assert(sizeof(UInt) <= sizeof(unsigned int));
  262. return __builtin_ctz((unsigned int)n);
  263. }
  264. #elif defined(_MSC_VER)
  265. #define JKJ_HAS_COUNTR_ZERO_INTRINSIC 1
  266. if constexpr (std::is_same_v<UInt, unsigned __int64>) {
  267. #if defined(_M_X64)
  268. return int(_tzcnt_u64(n));
  269. #else
  270. return ((unsigned int)(n) == 0) ?
  271. (32 + (_tzcnt_u32((unsigned int)(n >> 32)))) :
  272. (_tzcnt_u32((unsigned int)n));
  273. #endif
  274. }
  275. else {
  276. static_assert(sizeof(UInt) <= sizeof(unsigned int));
  277. return int(_tzcnt_u32((unsigned int)n));
  278. }
  279. #else
  280. #define JKJ_HAS_COUNTR_ZERO_INTRINSIC 0
  281. int count = int(value_bits<UInt>);
  282. auto n32 = std::uint32_t(n);
  283. if constexpr (value_bits<UInt> > 32) {
  284. if (n32 != 0) {
  285. count = 31;
  286. }
  287. else {
  288. n32 = std::uint32_t(n >> 32);
  289. if (n32 != 0) {
  290. count -= 1;
  291. }
  292. }
  293. }
  294. if constexpr (value_bits<UInt> > 16) {
  295. if ((n32 & 0x0000ffff) != 0) count -= 16;
  296. }
  297. if constexpr (value_bits<UInt> > 8) {
  298. if ((n32 & 0x00ff00ff) != 0) count -= 8;
  299. }
  300. if ((n32 & 0x0f0f0f0f) != 0) count -= 4;
  301. if ((n32 & 0x33333333) != 0) count -= 2;
  302. if ((n32 & 0x55555555) != 0) count -= 1;
  303. return count;
  304. #endif
  305. }
  306. }
  307. ////////////////////////////////////////////////////////////////////////////////////////
  308. // Utilities for wide unsigned integer arithmetic
  309. ////////////////////////////////////////////////////////////////////////////////////////
  310. namespace wuint {
  311. struct uint128 {
  312. uint128() = default;
  313. #if (defined(__GNUC__) || defined(__clang__)) && defined(__SIZEOF_INT128__) && defined(__x86_64__)
  314. unsigned __int128 internal_;
  315. constexpr uint128(std::uint64_t high, std::uint64_t low) noexcept :
  316. internal_{ ((unsigned __int128)low) | (((unsigned __int128)high) << 64) } {}
  317. constexpr uint128(unsigned __int128 u) noexcept : internal_{ u } {}
  318. constexpr std::uint64_t high() const noexcept {
  319. return std::uint64_t(internal_ >> 64);
  320. }
  321. constexpr std::uint64_t low() const noexcept {
  322. return std::uint64_t(internal_);
  323. }
  324. uint128& operator+=(std::uint64_t n) & noexcept {
  325. internal_ += n;
  326. return *this;
  327. }
  328. #else
  329. std::uint64_t high_;
  330. std::uint64_t low_;
  331. constexpr uint128(std::uint64_t high, std::uint64_t low) noexcept :
  332. high_{ high }, low_{ low } {}
  333. constexpr std::uint64_t high() const noexcept {
  334. return high_;
  335. }
  336. constexpr std::uint64_t low() const noexcept {
  337. return low_;
  338. }
  339. uint128& operator+=(std::uint64_t n) & noexcept {
  340. #if defined(_MSC_VER) && defined(_M_X64)
  341. auto carry = _addcarry_u64(0, low_, n, &low_);
  342. _addcarry_u64(carry, high_, 0, &high_);
  343. return *this;
  344. #else
  345. auto sum = low_ + n;
  346. high_ += (sum < low_ ? 1 : 0);
  347. low_ = sum;
  348. return *this;
  349. #endif
  350. }
  351. #endif
  352. };
  353. static inline std::uint64_t umul64(std::uint32_t x, std::uint32_t y) noexcept {
  354. #if defined(_MSC_VER) && defined(_M_IX86)
  355. return __emulu(x, y);
  356. #else
  357. return x * std::uint64_t(y);
  358. #endif
  359. }
  360. // Get 128-bit result of multiplication of two 64-bit unsigned integers
  361. JKJ_SAFEBUFFERS inline uint128 umul128(std::uint64_t x, std::uint64_t y) noexcept {
  362. #if (defined(__GNUC__) || defined(__clang__)) && defined(__SIZEOF_INT128__) && defined(__x86_64__)
  363. return (unsigned __int128)(x) * (unsigned __int128)(y);
  364. #elif defined(_MSC_VER) && defined(_M_X64)
  365. uint128 result;
  366. result.low_ = _umul128(x, y, &result.high_);
  367. return result;
  368. #else
  369. auto a = std::uint32_t(x >> 32);
  370. auto b = std::uint32_t(x);
  371. auto c = std::uint32_t(y >> 32);
  372. auto d = std::uint32_t(y);
  373. auto ac = umul64(a, c);
  374. auto bc = umul64(b, c);
  375. auto ad = umul64(a, d);
  376. auto bd = umul64(b, d);
  377. auto intermediate = (bd >> 32) + std::uint32_t(ad) + std::uint32_t(bc);
  378. return{ ac + (intermediate >> 32) + (ad >> 32) + (bc >> 32),
  379. (intermediate << 32) + std::uint32_t(bd) };
  380. #endif
  381. }
  382. JKJ_SAFEBUFFERS inline std::uint64_t umul128_upper64(std::uint64_t x, std::uint64_t y) noexcept {
  383. #if (defined(__GNUC__) || defined(__clang__)) && defined(__SIZEOF_INT128__) && defined(__x86_64__)
  384. auto p = (unsigned __int128)(x) * (unsigned __int128)(y);
  385. return std::uint64_t(p >> 64);
  386. #elif defined(_MSC_VER) && defined(_M_X64)
  387. return __umulh(x, y);
  388. #else
  389. auto a = std::uint32_t(x >> 32);
  390. auto b = std::uint32_t(x);
  391. auto c = std::uint32_t(y >> 32);
  392. auto d = std::uint32_t(y);
  393. auto ac = umul64(a, c);
  394. auto bc = umul64(b, c);
  395. auto ad = umul64(a, d);
  396. auto bd = umul64(b, d);
  397. auto intermediate = (bd >> 32) + std::uint32_t(ad) + std::uint32_t(bc);
  398. return ac + (intermediate >> 32) + (ad >> 32) + (bc >> 32);
  399. #endif
  400. }
  401. // Get upper 64-bits of multiplication of a 64-bit unsigned integer and a 128-bit unsigned integer
  402. JKJ_SAFEBUFFERS inline std::uint64_t umul192_upper64(std::uint64_t x, uint128 y) noexcept {
  403. auto g0 = umul128(x, y.high());
  404. g0 += umul128_upper64(x, y.low());
  405. return g0.high();
  406. }
  407. // Get upper 32-bits of multiplication of a 32-bit unsigned integer and a 64-bit unsigned integer
  408. inline std::uint32_t umul96_upper32(std::uint32_t x, std::uint64_t y) noexcept {
  409. #if defined(__x86_64__) || defined(_M_X64)
  410. return std::uint32_t(umul128_upper64(x, y));
  411. #else
  412. //std::uint32_t a = 0;
  413. auto b = x;
  414. auto c = std::uint32_t(y >> 32);
  415. auto d = std::uint32_t(y);
  416. //std::uint64_t ac = 0;
  417. auto bc = umul64(b, c);
  418. //std::uint64_t ad = 0;
  419. auto bd = umul64(b, d);
  420. auto intermediate = (bd >> 32) + bc;
  421. return std::uint32_t(intermediate >> 32);
  422. #endif
  423. }
  424. // Get middle 64-bits of multiplication of a 64-bit unsigned integer and a 128-bit unsigned integer
  425. JKJ_SAFEBUFFERS inline std::uint64_t umul192_middle64(std::uint64_t x, uint128 y) noexcept {
  426. auto g01 = x * y.high();
  427. auto g10 = umul128_upper64(x, y.low());
  428. return g01 + g10;
  429. }
  430. // Get middle 32-bits of multiplication of a 32-bit unsigned integer and a 64-bit unsigned integer
  431. inline std::uint64_t umul96_lower64(std::uint32_t x, std::uint64_t y) noexcept {
  432. return x * y;
  433. }
  434. }
  435. template <int k, class Int>
  436. constexpr Int compute_power(Int a) noexcept {
  437. static_assert(k >= 0);
  438. Int p = 1;
  439. for (int i = 0; i < k; ++i) {
  440. p *= a;
  441. }
  442. return p;
  443. }
  444. template <int a, class UInt>
  445. constexpr int count_factors(UInt n) noexcept {
  446. static_assert(a > 1);
  447. int c = 0;
  448. while (n % a == 0) {
  449. n /= a;
  450. ++c;
  451. }
  452. return c;
  453. }
  454. ////////////////////////////////////////////////////////////////////////////////////////
  455. // Utilities for fast/constexpr log computation
  456. ////////////////////////////////////////////////////////////////////////////////////////
  457. namespace log {
  458. constexpr std::int32_t floor_shift(
  459. std::uint32_t integer_part,
  460. std::uint64_t fractional_digits,
  461. std::size_t shift_amount) noexcept
  462. {
  463. assert(shift_amount < 32);
  464. // Ensure no overflow
  465. assert(shift_amount == 0 || integer_part < (std::uint32_t(1) << (32 - shift_amount)));
  466. return shift_amount == 0 ? std::int32_t(integer_part) :
  467. std::int32_t(
  468. (integer_part << shift_amount) |
  469. (fractional_digits >> (64 - shift_amount)));
  470. }
  471. // Compute floor(e * c - s)
  472. template <
  473. std::uint32_t c_integer_part,
  474. std::uint64_t c_fractional_digits,
  475. std::size_t shift_amount,
  476. std::int32_t max_exponent,
  477. std::uint32_t s_integer_part = 0,
  478. std::uint64_t s_fractional_digits = 0
  479. >
  480. constexpr int compute(int e) noexcept {
  481. assert(e <= max_exponent && e >= -max_exponent);
  482. constexpr auto c = floor_shift(c_integer_part, c_fractional_digits, shift_amount);
  483. constexpr auto s = floor_shift(s_integer_part, s_fractional_digits, shift_amount);
  484. return int((std::int32_t(e) * c - s) >> shift_amount);
  485. }
  486. static constexpr std::uint64_t log10_2_fractional_digits{ 0x4d10'4d42'7de7'fbcc };
  487. static constexpr std::uint64_t log10_4_over_3_fractional_digits{ 0x1ffb'fc2b'bc78'0375 };
  488. static constexpr std::size_t floor_log10_pow2_shift_amount = 22;
  489. static constexpr int floor_log10_pow2_input_limit = 1700;
  490. static constexpr int floor_log10_pow2_minus_log10_4_over_3_input_limit = 1700;
  491. static constexpr std::uint64_t log2_10_fractional_digits{ 0x5269'e12f'346e'2bf9 };
  492. static constexpr std::size_t floor_log2_pow10_shift_amount = 19;
  493. static constexpr int floor_log2_pow10_input_limit = 1233;
  494. static constexpr std::uint64_t log5_2_fractional_digits{ 0x6e40'd1a4'143d'cb94 };
  495. static constexpr std::uint64_t log5_3_fractional_digits{ 0xaebf'4791'5d44'3b24 };
  496. static constexpr std::size_t floor_log5_pow2_shift_amount = 20;
  497. static constexpr int floor_log5_pow2_input_limit = 1492;
  498. static constexpr int floor_log5_pow2_minus_log5_3_input_limit = 2427;
  499. // For constexpr computation
  500. // Returns -1 when n = 0
  501. template <class UInt>
  502. constexpr int floor_log2(UInt n) noexcept {
  503. int count = -1;
  504. while (n != 0) {
  505. ++count;
  506. n >>= 1;
  507. }
  508. return count;
  509. }
  510. constexpr int floor_log10_pow2(int e) noexcept {
  511. using namespace log;
  512. return compute<
  513. 0, log10_2_fractional_digits,
  514. floor_log10_pow2_shift_amount,
  515. floor_log10_pow2_input_limit>(e);
  516. }
  517. constexpr int floor_log2_pow10(int e) noexcept {
  518. using namespace log;
  519. return compute<
  520. 3, log2_10_fractional_digits,
  521. floor_log2_pow10_shift_amount,
  522. floor_log2_pow10_input_limit>(e);
  523. }
  524. constexpr int floor_log5_pow2(int e) noexcept {
  525. using namespace log;
  526. return compute<
  527. 0, log5_2_fractional_digits,
  528. floor_log5_pow2_shift_amount,
  529. floor_log5_pow2_input_limit>(e);
  530. }
  531. constexpr int floor_log5_pow2_minus_log5_3(int e) noexcept {
  532. using namespace log;
  533. return compute<
  534. 0, log5_2_fractional_digits,
  535. floor_log5_pow2_shift_amount,
  536. floor_log5_pow2_minus_log5_3_input_limit,
  537. 0, log5_3_fractional_digits>(e);
  538. }
  539. constexpr int floor_log10_pow2_minus_log10_4_over_3(int e) noexcept {
  540. using namespace log;
  541. return compute<
  542. 0, log10_2_fractional_digits,
  543. floor_log10_pow2_shift_amount,
  544. floor_log10_pow2_minus_log10_4_over_3_input_limit,
  545. 0, log10_4_over_3_fractional_digits>(e);
  546. }
  547. }
  548. ////////////////////////////////////////////////////////////////////////////////////////
  549. // Utilities for fast divisibility test
  550. ////////////////////////////////////////////////////////////////////////////////////////
  551. namespace div {
  552. template <class UInt, UInt a>
  553. constexpr UInt modular_inverse(int bit_width = int(value_bits<UInt>)) noexcept {
  554. // By Euler's theorem, a^phi(2^n) == 1 (mod 2^n),
  555. // where phi(2^n) = 2^(n-1), so the modular inverse of a is
  556. // a^(2^(n-1) - 1) = a^(1 + 2 + 2^2 + ... + 2^(n-2))
  557. std::common_type_t<UInt, unsigned int> mod_inverse = 1;
  558. for (int i = 1; i < bit_width; ++i) {
  559. mod_inverse = mod_inverse * mod_inverse * a;
  560. }
  561. if (bit_width < value_bits<UInt>) {
  562. auto mask = UInt((UInt(1) << bit_width) - 1);
  563. return UInt(mod_inverse & mask);
  564. }
  565. else {
  566. return UInt(mod_inverse);
  567. }
  568. }
  569. template <class UInt, UInt a, int N>
  570. struct table_t {
  571. static_assert(std::is_unsigned_v<UInt>);
  572. static_assert(a % 2 != 0);
  573. static_assert(N > 0);
  574. static constexpr int size = N;
  575. UInt mod_inv[N];
  576. UInt max_quotients[N];
  577. };
  578. template <class UInt, UInt a, int N>
  579. struct table_holder {
  580. static constexpr table_t<UInt, a, N> table = [] {
  581. constexpr auto mod_inverse = modular_inverse<UInt, a>();
  582. table_t<UInt, a, N> table{};
  583. std::common_type_t<UInt, unsigned int> pow_of_mod_inverse = 1;
  584. UInt pow_of_a = 1;
  585. for (int i = 0; i < N; ++i) {
  586. table.mod_inv[i] = UInt(pow_of_mod_inverse);
  587. table.max_quotients[i] = UInt(std::numeric_limits<UInt>::max() / pow_of_a);
  588. pow_of_mod_inverse *= mod_inverse;
  589. pow_of_a *= a;
  590. }
  591. return table;
  592. }();
  593. };
  594. template <std::size_t table_size, class UInt>
  595. constexpr bool divisible_by_power_of_5(UInt x, unsigned int exp) noexcept {
  596. auto const& table = table_holder<UInt, 5, table_size>::table;
  597. assert(exp < (unsigned int)(table.size));
  598. return (x * table.mod_inv[exp]) <= table.max_quotients[exp];
  599. }
  600. template <class UInt>
  601. constexpr bool divisible_by_power_of_2(UInt x, unsigned int exp) noexcept {
  602. assert(exp >= 1);
  603. assert(x != 0);
  604. #if JKJ_HAS_COUNTR_ZERO_INTRINSIC
  605. return bits::countr_zero(x) >= int(exp);
  606. #else
  607. if (exp >= int(value_bits<UInt>)) {
  608. return false;
  609. }
  610. auto mask = UInt((UInt(1) << exp) - 1);
  611. return (x & mask) == 0;
  612. #endif
  613. }
  614. // Replace n by floor(n / 5^N)
  615. // Returns true if and only if n is divisible by 5^N
  616. // Precondition: n <= 2 * 5^(N+1)
  617. template <int N>
  618. struct check_divisibility_and_divide_by_pow5_info;
  619. template <>
  620. struct check_divisibility_and_divide_by_pow5_info<1> {
  621. static constexpr std::uint32_t magic_number = 0xcccd;
  622. static constexpr int bits_for_comparison = 16;
  623. static constexpr std::uint32_t threshold = 0x3333;
  624. static constexpr int shift_amount = 18;
  625. };
  626. template <>
  627. struct check_divisibility_and_divide_by_pow5_info<2> {
  628. static constexpr std::uint32_t magic_number = 0xa429;
  629. static constexpr int bits_for_comparison = 8;
  630. static constexpr std::uint32_t threshold = 0x0a;
  631. static constexpr int shift_amount = 20;
  632. };
  633. template <int N>
  634. constexpr bool check_divisibility_and_divide_by_pow5(std::uint32_t& n) noexcept
  635. {
  636. // Make sure the computation for max_n does not overflow
  637. static_assert(N + 1 <= log::floor_log5_pow2(31));
  638. assert(n <= compute_power<N + 1>(std::uint32_t(5)) * 2);
  639. using info = check_divisibility_and_divide_by_pow5_info<N>;
  640. n *= info::magic_number;
  641. constexpr std::uint32_t comparison_mask =
  642. info::bits_for_comparison >= 32 ? std::numeric_limits<std::uint32_t>::max() :
  643. std::uint32_t((std::uint32_t(1) << info::bits_for_comparison) - 1);
  644. if ((n & comparison_mask) <= info::threshold) {
  645. n >>= info::shift_amount;
  646. return true;
  647. }
  648. else {
  649. n >>= info::shift_amount;
  650. return false;
  651. }
  652. }
  653. // Compute floor(n / 10^N) for small n and N
  654. // Precondition: n <= 10^(N+1)
  655. template <int N>
  656. struct small_division_by_pow10_info;
  657. template <>
  658. struct small_division_by_pow10_info<1> {
  659. static constexpr std::uint32_t magic_number = 0xcccd;
  660. static constexpr int shift_amount = 19;
  661. };
  662. template <>
  663. struct small_division_by_pow10_info<2> {
  664. static constexpr std::uint32_t magic_number = 0xa3d8;
  665. static constexpr int shift_amount = 22;
  666. };
  667. template <int N>
  668. constexpr std::uint32_t small_division_by_pow10(std::uint32_t n) noexcept
  669. {
  670. assert(n <= compute_power<N + 1>(std::uint32_t(10)));
  671. return (n * small_division_by_pow10_info<N>::magic_number)
  672. >> small_division_by_pow10_info<N>::shift_amount;
  673. }
  674. // Compute floor(n / 10^N) for small N
  675. // Precondition: n <= 2^a * 5^b (a = max_pow2, b = max_pow5)
  676. template <int N, int max_pow2, int max_pow5, class UInt>
  677. constexpr UInt divide_by_pow10(UInt n) noexcept
  678. {
  679. static_assert(N >= 0);
  680. // Ensure no overflow
  681. static_assert(max_pow2 + (log::floor_log2_pow10(max_pow5) - max_pow5) < value_bits<UInt>);
  682. // Specialize for 64bit division by 1000
  683. // Ensure that the correctness condition is met
  684. if constexpr (std::is_same_v<UInt, std::uint64_t> && N == 3 &&
  685. max_pow2 + (log::floor_log2_pow10(N + max_pow5) - (N + max_pow5)) < 70)
  686. {
  687. return wuint::umul128_upper64(n, 0x8312'6e97'8d4f'df3c) >> 9;
  688. }
  689. else {
  690. constexpr auto divisor = compute_power<N>(UInt(10));
  691. return n / divisor;
  692. }
  693. }
  694. }
  695. }
  696. ////////////////////////////////////////////////////////////////////////////////////////
  697. // DIY floating-point data type
  698. ////////////////////////////////////////////////////////////////////////////////////////
  699. template <class Float, bool is_signed, bool trailing_zero_flag>
  700. struct fp_t;
  701. template <class Float>
  702. struct fp_t<Float, false, false> {
  703. using float_type = Float;
  704. using carrier_uint = typename ieee754_traits<Float>::carrier_uint;
  705. carrier_uint significand;
  706. int exponent;
  707. };
  708. template <class Float>
  709. struct fp_t<Float, true, false> {
  710. using float_type = Float;
  711. using carrier_uint = typename ieee754_traits<Float>::carrier_uint;
  712. carrier_uint significand;
  713. int exponent;
  714. bool is_negative;
  715. };
  716. template <class Float>
  717. struct fp_t<Float, false, true> {
  718. using float_type = Float;
  719. using carrier_uint = typename ieee754_traits<Float>::carrier_uint;
  720. carrier_uint significand;
  721. int exponent;
  722. bool may_have_trailing_zeros;
  723. };
  724. template <class Float>
  725. struct fp_t<Float, true, true> {
  726. using float_type = Float;
  727. using carrier_uint = typename ieee754_traits<Float>::carrier_uint;
  728. carrier_uint significand;
  729. int exponent;
  730. bool is_negative;
  731. bool may_have_trailing_zeros;
  732. };
  733. template <class Float>
  734. using unsigned_fp_t = fp_t<Float, false, false>;
  735. template <class Float>
  736. using signed_fp_t = fp_t<Float, true, false>;
  737. ////////////////////////////////////////////////////////////////////////////////////////
  738. // Computed cache entries
  739. ////////////////////////////////////////////////////////////////////////////////////////
  740. namespace detail {
  741. template <ieee754_format format>
  742. struct cache_holder;
  743. template <>
  744. struct cache_holder<ieee754_format::binary32> {
  745. using cache_entry_type = std::uint64_t;
  746. static constexpr int cache_bits = 64;
  747. static constexpr int min_k = -31;
  748. static constexpr int max_k = 46;
  749. static constexpr cache_entry_type cache[] = {
  750. 0x81ceb32c4b43fcf5,
  751. 0xa2425ff75e14fc32,
  752. 0xcad2f7f5359a3b3f,
  753. 0xfd87b5f28300ca0e,
  754. 0x9e74d1b791e07e49,
  755. 0xc612062576589ddb,
  756. 0xf79687aed3eec552,
  757. 0x9abe14cd44753b53,
  758. 0xc16d9a0095928a28,
  759. 0xf1c90080baf72cb2,
  760. 0x971da05074da7bef,
  761. 0xbce5086492111aeb,
  762. 0xec1e4a7db69561a6,
  763. 0x9392ee8e921d5d08,
  764. 0xb877aa3236a4b44a,
  765. 0xe69594bec44de15c,
  766. 0x901d7cf73ab0acda,
  767. 0xb424dc35095cd810,
  768. 0xe12e13424bb40e14,
  769. 0x8cbccc096f5088cc,
  770. 0xafebff0bcb24aaff,
  771. 0xdbe6fecebdedd5bf,
  772. 0x89705f4136b4a598,
  773. 0xabcc77118461cefd,
  774. 0xd6bf94d5e57a42bd,
  775. 0x8637bd05af6c69b6,
  776. 0xa7c5ac471b478424,
  777. 0xd1b71758e219652c,
  778. 0x83126e978d4fdf3c,
  779. 0xa3d70a3d70a3d70b,
  780. 0xcccccccccccccccd,
  781. 0x8000000000000000,
  782. 0xa000000000000000,
  783. 0xc800000000000000,
  784. 0xfa00000000000000,
  785. 0x9c40000000000000,
  786. 0xc350000000000000,
  787. 0xf424000000000000,
  788. 0x9896800000000000,
  789. 0xbebc200000000000,
  790. 0xee6b280000000000,
  791. 0x9502f90000000000,
  792. 0xba43b74000000000,
  793. 0xe8d4a51000000000,
  794. 0x9184e72a00000000,
  795. 0xb5e620f480000000,
  796. 0xe35fa931a0000000,
  797. 0x8e1bc9bf04000000,
  798. 0xb1a2bc2ec5000000,
  799. 0xde0b6b3a76400000,
  800. 0x8ac7230489e80000,
  801. 0xad78ebc5ac620000,
  802. 0xd8d726b7177a8000,
  803. 0x878678326eac9000,
  804. 0xa968163f0a57b400,
  805. 0xd3c21bcecceda100,
  806. 0x84595161401484a0,
  807. 0xa56fa5b99019a5c8,
  808. 0xcecb8f27f4200f3a,
  809. 0x813f3978f8940984,
  810. 0xa18f07d736b90be5,
  811. 0xc9f2c9cd04674ede,
  812. 0xfc6f7c4045812296,
  813. 0x9dc5ada82b70b59d,
  814. 0xc5371912364ce305,
  815. 0xf684df56c3e01bc6,
  816. 0x9a130b963a6c115c,
  817. 0xc097ce7bc90715b3,
  818. 0xf0bdc21abb48db20,
  819. 0x96769950b50d88f4,
  820. 0xbc143fa4e250eb31,
  821. 0xeb194f8e1ae525fd,
  822. 0x92efd1b8d0cf37be,
  823. 0xb7abc627050305ad,
  824. 0xe596b7b0c643c719,
  825. 0x8f7e32ce7bea5c6f,
  826. 0xb35dbf821ae4f38b,
  827. 0xe0352f62a19e306e
  828. };
  829. };
  830. template <>
  831. struct cache_holder<ieee754_format::binary64> {
  832. using cache_entry_type = wuint::uint128;
  833. static constexpr int cache_bits = 128;
  834. static constexpr int min_k = -292;
  835. static constexpr int max_k = 326;
  836. static constexpr cache_entry_type cache[] = {
  837. { 0xff77b1fcbebcdc4f, 0x25e8e89c13bb0f7b },
  838. { 0x9faacf3df73609b1, 0x77b191618c54e9ad },
  839. { 0xc795830d75038c1d, 0xd59df5b9ef6a2418 },
  840. { 0xf97ae3d0d2446f25, 0x4b0573286b44ad1e },
  841. { 0x9becce62836ac577, 0x4ee367f9430aec33 },
  842. { 0xc2e801fb244576d5, 0x229c41f793cda740 },
  843. { 0xf3a20279ed56d48a, 0x6b43527578c11110 },
  844. { 0x9845418c345644d6, 0x830a13896b78aaaa },
  845. { 0xbe5691ef416bd60c, 0x23cc986bc656d554 },
  846. { 0xedec366b11c6cb8f, 0x2cbfbe86b7ec8aa9 },
  847. { 0x94b3a202eb1c3f39, 0x7bf7d71432f3d6aa },
  848. { 0xb9e08a83a5e34f07, 0xdaf5ccd93fb0cc54 },
  849. { 0xe858ad248f5c22c9, 0xd1b3400f8f9cff69 },
  850. { 0x91376c36d99995be, 0x23100809b9c21fa2 },
  851. { 0xb58547448ffffb2d, 0xabd40a0c2832a78b },
  852. { 0xe2e69915b3fff9f9, 0x16c90c8f323f516d },
  853. { 0x8dd01fad907ffc3b, 0xae3da7d97f6792e4 },
  854. { 0xb1442798f49ffb4a, 0x99cd11cfdf41779d },
  855. { 0xdd95317f31c7fa1d, 0x40405643d711d584 },
  856. { 0x8a7d3eef7f1cfc52, 0x482835ea666b2573 },
  857. { 0xad1c8eab5ee43b66, 0xda3243650005eed0 },
  858. { 0xd863b256369d4a40, 0x90bed43e40076a83 },
  859. { 0x873e4f75e2224e68, 0x5a7744a6e804a292 },
  860. { 0xa90de3535aaae202, 0x711515d0a205cb37 },
  861. { 0xd3515c2831559a83, 0x0d5a5b44ca873e04 },
  862. { 0x8412d9991ed58091, 0xe858790afe9486c3 },
  863. { 0xa5178fff668ae0b6, 0x626e974dbe39a873 },
  864. { 0xce5d73ff402d98e3, 0xfb0a3d212dc81290 },
  865. { 0x80fa687f881c7f8e, 0x7ce66634bc9d0b9a },
  866. { 0xa139029f6a239f72, 0x1c1fffc1ebc44e81 },
  867. { 0xc987434744ac874e, 0xa327ffb266b56221 },
  868. { 0xfbe9141915d7a922, 0x4bf1ff9f0062baa9 },
  869. { 0x9d71ac8fada6c9b5, 0x6f773fc3603db4aa },
  870. { 0xc4ce17b399107c22, 0xcb550fb4384d21d4 },
  871. { 0xf6019da07f549b2b, 0x7e2a53a146606a49 },
  872. { 0x99c102844f94e0fb, 0x2eda7444cbfc426e },
  873. { 0xc0314325637a1939, 0xfa911155fefb5309 },
  874. { 0xf03d93eebc589f88, 0x793555ab7eba27cb },
  875. { 0x96267c7535b763b5, 0x4bc1558b2f3458df },
  876. { 0xbbb01b9283253ca2, 0x9eb1aaedfb016f17 },
  877. { 0xea9c227723ee8bcb, 0x465e15a979c1cadd },
  878. { 0x92a1958a7675175f, 0x0bfacd89ec191eca },
  879. { 0xb749faed14125d36, 0xcef980ec671f667c },
  880. { 0xe51c79a85916f484, 0x82b7e12780e7401b },
  881. { 0x8f31cc0937ae58d2, 0xd1b2ecb8b0908811 },
  882. { 0xb2fe3f0b8599ef07, 0x861fa7e6dcb4aa16 },
  883. { 0xdfbdcece67006ac9, 0x67a791e093e1d49b },
  884. { 0x8bd6a141006042bd, 0xe0c8bb2c5c6d24e1 },
  885. { 0xaecc49914078536d, 0x58fae9f773886e19 },
  886. { 0xda7f5bf590966848, 0xaf39a475506a899f },
  887. { 0x888f99797a5e012d, 0x6d8406c952429604 },
  888. { 0xaab37fd7d8f58178, 0xc8e5087ba6d33b84 },
  889. { 0xd5605fcdcf32e1d6, 0xfb1e4a9a90880a65 },
  890. { 0x855c3be0a17fcd26, 0x5cf2eea09a550680 },
  891. { 0xa6b34ad8c9dfc06f, 0xf42faa48c0ea481f },
  892. { 0xd0601d8efc57b08b, 0xf13b94daf124da27 },
  893. { 0x823c12795db6ce57, 0x76c53d08d6b70859 },
  894. { 0xa2cb1717b52481ed, 0x54768c4b0c64ca6f },
  895. { 0xcb7ddcdda26da268, 0xa9942f5dcf7dfd0a },
  896. { 0xfe5d54150b090b02, 0xd3f93b35435d7c4d },
  897. { 0x9efa548d26e5a6e1, 0xc47bc5014a1a6db0 },
  898. { 0xc6b8e9b0709f109a, 0x359ab6419ca1091c },
  899. { 0xf867241c8cc6d4c0, 0xc30163d203c94b63 },
  900. { 0x9b407691d7fc44f8, 0x79e0de63425dcf1e },
  901. { 0xc21094364dfb5636, 0x985915fc12f542e5 },
  902. { 0xf294b943e17a2bc4, 0x3e6f5b7b17b2939e },
  903. { 0x979cf3ca6cec5b5a, 0xa705992ceecf9c43 },
  904. { 0xbd8430bd08277231, 0x50c6ff782a838354 },
  905. { 0xece53cec4a314ebd, 0xa4f8bf5635246429 },
  906. { 0x940f4613ae5ed136, 0x871b7795e136be9a },
  907. { 0xb913179899f68584, 0x28e2557b59846e40 },
  908. { 0xe757dd7ec07426e5, 0x331aeada2fe589d0 },
  909. { 0x9096ea6f3848984f, 0x3ff0d2c85def7622 },
  910. { 0xb4bca50b065abe63, 0x0fed077a756b53aa },
  911. { 0xe1ebce4dc7f16dfb, 0xd3e8495912c62895 },
  912. { 0x8d3360f09cf6e4bd, 0x64712dd7abbbd95d },
  913. { 0xb080392cc4349dec, 0xbd8d794d96aacfb4 },
  914. { 0xdca04777f541c567, 0xecf0d7a0fc5583a1 },
  915. { 0x89e42caaf9491b60, 0xf41686c49db57245 },
  916. { 0xac5d37d5b79b6239, 0x311c2875c522ced6 },
  917. { 0xd77485cb25823ac7, 0x7d633293366b828c },
  918. { 0x86a8d39ef77164bc, 0xae5dff9c02033198 },
  919. { 0xa8530886b54dbdeb, 0xd9f57f830283fdfd },
  920. { 0xd267caa862a12d66, 0xd072df63c324fd7c },
  921. { 0x8380dea93da4bc60, 0x4247cb9e59f71e6e },
  922. { 0xa46116538d0deb78, 0x52d9be85f074e609 },
  923. { 0xcd795be870516656, 0x67902e276c921f8c },
  924. { 0x806bd9714632dff6, 0x00ba1cd8a3db53b7 },
  925. { 0xa086cfcd97bf97f3, 0x80e8a40eccd228a5 },
  926. { 0xc8a883c0fdaf7df0, 0x6122cd128006b2ce },
  927. { 0xfad2a4b13d1b5d6c, 0x796b805720085f82 },
  928. { 0x9cc3a6eec6311a63, 0xcbe3303674053bb1 },
  929. { 0xc3f490aa77bd60fc, 0xbedbfc4411068a9d },
  930. { 0xf4f1b4d515acb93b, 0xee92fb5515482d45 },
  931. { 0x991711052d8bf3c5, 0x751bdd152d4d1c4b },
  932. { 0xbf5cd54678eef0b6, 0xd262d45a78a0635e },
  933. { 0xef340a98172aace4, 0x86fb897116c87c35 },
  934. { 0x9580869f0e7aac0e, 0xd45d35e6ae3d4da1 },
  935. { 0xbae0a846d2195712, 0x8974836059cca10a },
  936. { 0xe998d258869facd7, 0x2bd1a438703fc94c },
  937. { 0x91ff83775423cc06, 0x7b6306a34627ddd0 },
  938. { 0xb67f6455292cbf08, 0x1a3bc84c17b1d543 },
  939. { 0xe41f3d6a7377eeca, 0x20caba5f1d9e4a94 },
  940. { 0x8e938662882af53e, 0x547eb47b7282ee9d },
  941. { 0xb23867fb2a35b28d, 0xe99e619a4f23aa44 },
  942. { 0xdec681f9f4c31f31, 0x6405fa00e2ec94d5 },
  943. { 0x8b3c113c38f9f37e, 0xde83bc408dd3dd05 },
  944. { 0xae0b158b4738705e, 0x9624ab50b148d446 },
  945. { 0xd98ddaee19068c76, 0x3badd624dd9b0958 },
  946. { 0x87f8a8d4cfa417c9, 0xe54ca5d70a80e5d7 },
  947. { 0xa9f6d30a038d1dbc, 0x5e9fcf4ccd211f4d },
  948. { 0xd47487cc8470652b, 0x7647c32000696720 },
  949. { 0x84c8d4dfd2c63f3b, 0x29ecd9f40041e074 },
  950. { 0xa5fb0a17c777cf09, 0xf468107100525891 },
  951. { 0xcf79cc9db955c2cc, 0x7182148d4066eeb5 },
  952. { 0x81ac1fe293d599bf, 0xc6f14cd848405531 },
  953. { 0xa21727db38cb002f, 0xb8ada00e5a506a7d },
  954. { 0xca9cf1d206fdc03b, 0xa6d90811f0e4851d },
  955. { 0xfd442e4688bd304a, 0x908f4a166d1da664 },
  956. { 0x9e4a9cec15763e2e, 0x9a598e4e043287ff },
  957. { 0xc5dd44271ad3cdba, 0x40eff1e1853f29fe },
  958. { 0xf7549530e188c128, 0xd12bee59e68ef47d },
  959. { 0x9a94dd3e8cf578b9, 0x82bb74f8301958cf },
  960. { 0xc13a148e3032d6e7, 0xe36a52363c1faf02 },
  961. { 0xf18899b1bc3f8ca1, 0xdc44e6c3cb279ac2 },
  962. { 0x96f5600f15a7b7e5, 0x29ab103a5ef8c0ba },
  963. { 0xbcb2b812db11a5de, 0x7415d448f6b6f0e8 },
  964. { 0xebdf661791d60f56, 0x111b495b3464ad22 },
  965. { 0x936b9fcebb25c995, 0xcab10dd900beec35 },
  966. { 0xb84687c269ef3bfb, 0x3d5d514f40eea743 },
  967. { 0xe65829b3046b0afa, 0x0cb4a5a3112a5113 },
  968. { 0x8ff71a0fe2c2e6dc, 0x47f0e785eaba72ac },
  969. { 0xb3f4e093db73a093, 0x59ed216765690f57 },
  970. { 0xe0f218b8d25088b8, 0x306869c13ec3532d },
  971. { 0x8c974f7383725573, 0x1e414218c73a13fc },
  972. { 0xafbd2350644eeacf, 0xe5d1929ef90898fb },
  973. { 0xdbac6c247d62a583, 0xdf45f746b74abf3a },
  974. { 0x894bc396ce5da772, 0x6b8bba8c328eb784 },
  975. { 0xab9eb47c81f5114f, 0x066ea92f3f326565 },
  976. { 0xd686619ba27255a2, 0xc80a537b0efefebe },
  977. { 0x8613fd0145877585, 0xbd06742ce95f5f37 },
  978. { 0xa798fc4196e952e7, 0x2c48113823b73705 },
  979. { 0xd17f3b51fca3a7a0, 0xf75a15862ca504c6 },
  980. { 0x82ef85133de648c4, 0x9a984d73dbe722fc },
  981. { 0xa3ab66580d5fdaf5, 0xc13e60d0d2e0ebbb },
  982. { 0xcc963fee10b7d1b3, 0x318df905079926a9 },
  983. { 0xffbbcfe994e5c61f, 0xfdf17746497f7053 },
  984. { 0x9fd561f1fd0f9bd3, 0xfeb6ea8bedefa634 },
  985. { 0xc7caba6e7c5382c8, 0xfe64a52ee96b8fc1 },
  986. { 0xf9bd690a1b68637b, 0x3dfdce7aa3c673b1 },
  987. { 0x9c1661a651213e2d, 0x06bea10ca65c084f },
  988. { 0xc31bfa0fe5698db8, 0x486e494fcff30a63 },
  989. { 0xf3e2f893dec3f126, 0x5a89dba3c3efccfb },
  990. { 0x986ddb5c6b3a76b7, 0xf89629465a75e01d },
  991. { 0xbe89523386091465, 0xf6bbb397f1135824 },
  992. { 0xee2ba6c0678b597f, 0x746aa07ded582e2d },
  993. { 0x94db483840b717ef, 0xa8c2a44eb4571cdd },
  994. { 0xba121a4650e4ddeb, 0x92f34d62616ce414 },
  995. { 0xe896a0d7e51e1566, 0x77b020baf9c81d18 },
  996. { 0x915e2486ef32cd60, 0x0ace1474dc1d122f },
  997. { 0xb5b5ada8aaff80b8, 0x0d819992132456bb },
  998. { 0xe3231912d5bf60e6, 0x10e1fff697ed6c6a },
  999. { 0x8df5efabc5979c8f, 0xca8d3ffa1ef463c2 },
  1000. { 0xb1736b96b6fd83b3, 0xbd308ff8a6b17cb3 },
  1001. { 0xddd0467c64bce4a0, 0xac7cb3f6d05ddbdf },
  1002. { 0x8aa22c0dbef60ee4, 0x6bcdf07a423aa96c },
  1003. { 0xad4ab7112eb3929d, 0x86c16c98d2c953c7 },
  1004. { 0xd89d64d57a607744, 0xe871c7bf077ba8b8 },
  1005. { 0x87625f056c7c4a8b, 0x11471cd764ad4973 },
  1006. { 0xa93af6c6c79b5d2d, 0xd598e40d3dd89bd0 },
  1007. { 0xd389b47879823479, 0x4aff1d108d4ec2c4 },
  1008. { 0x843610cb4bf160cb, 0xcedf722a585139bb },
  1009. { 0xa54394fe1eedb8fe, 0xc2974eb4ee658829 },
  1010. { 0xce947a3da6a9273e, 0x733d226229feea33 },
  1011. { 0x811ccc668829b887, 0x0806357d5a3f5260 },
  1012. { 0xa163ff802a3426a8, 0xca07c2dcb0cf26f8 },
  1013. { 0xc9bcff6034c13052, 0xfc89b393dd02f0b6 },
  1014. { 0xfc2c3f3841f17c67, 0xbbac2078d443ace3 },
  1015. { 0x9d9ba7832936edc0, 0xd54b944b84aa4c0e },
  1016. { 0xc5029163f384a931, 0x0a9e795e65d4df12 },
  1017. { 0xf64335bcf065d37d, 0x4d4617b5ff4a16d6 },
  1018. { 0x99ea0196163fa42e, 0x504bced1bf8e4e46 },
  1019. { 0xc06481fb9bcf8d39, 0xe45ec2862f71e1d7 },
  1020. { 0xf07da27a82c37088, 0x5d767327bb4e5a4d },
  1021. { 0x964e858c91ba2655, 0x3a6a07f8d510f870 },
  1022. { 0xbbe226efb628afea, 0x890489f70a55368c },
  1023. { 0xeadab0aba3b2dbe5, 0x2b45ac74ccea842f },
  1024. { 0x92c8ae6b464fc96f, 0x3b0b8bc90012929e },
  1025. { 0xb77ada0617e3bbcb, 0x09ce6ebb40173745 },
  1026. { 0xe55990879ddcaabd, 0xcc420a6a101d0516 },
  1027. { 0x8f57fa54c2a9eab6, 0x9fa946824a12232e },
  1028. { 0xb32df8e9f3546564, 0x47939822dc96abfa },
  1029. { 0xdff9772470297ebd, 0x59787e2b93bc56f8 },
  1030. { 0x8bfbea76c619ef36, 0x57eb4edb3c55b65b },
  1031. { 0xaefae51477a06b03, 0xede622920b6b23f2 },
  1032. { 0xdab99e59958885c4, 0xe95fab368e45ecee },
  1033. { 0x88b402f7fd75539b, 0x11dbcb0218ebb415 },
  1034. { 0xaae103b5fcd2a881, 0xd652bdc29f26a11a },
  1035. { 0xd59944a37c0752a2, 0x4be76d3346f04960 },
  1036. { 0x857fcae62d8493a5, 0x6f70a4400c562ddc },
  1037. { 0xa6dfbd9fb8e5b88e, 0xcb4ccd500f6bb953 },
  1038. { 0xd097ad07a71f26b2, 0x7e2000a41346a7a8 },
  1039. { 0x825ecc24c873782f, 0x8ed400668c0c28c9 },
  1040. { 0xa2f67f2dfa90563b, 0x728900802f0f32fb },
  1041. { 0xcbb41ef979346bca, 0x4f2b40a03ad2ffba },
  1042. { 0xfea126b7d78186bc, 0xe2f610c84987bfa9 },
  1043. { 0x9f24b832e6b0f436, 0x0dd9ca7d2df4d7ca },
  1044. { 0xc6ede63fa05d3143, 0x91503d1c79720dbc },
  1045. { 0xf8a95fcf88747d94, 0x75a44c6397ce912b },
  1046. { 0x9b69dbe1b548ce7c, 0xc986afbe3ee11abb },
  1047. { 0xc24452da229b021b, 0xfbe85badce996169 },
  1048. { 0xf2d56790ab41c2a2, 0xfae27299423fb9c4 },
  1049. { 0x97c560ba6b0919a5, 0xdccd879fc967d41b },
  1050. { 0xbdb6b8e905cb600f, 0x5400e987bbc1c921 },
  1051. { 0xed246723473e3813, 0x290123e9aab23b69 },
  1052. { 0x9436c0760c86e30b, 0xf9a0b6720aaf6522 },
  1053. { 0xb94470938fa89bce, 0xf808e40e8d5b3e6a },
  1054. { 0xe7958cb87392c2c2, 0xb60b1d1230b20e05 },
  1055. { 0x90bd77f3483bb9b9, 0xb1c6f22b5e6f48c3 },
  1056. { 0xb4ecd5f01a4aa828, 0x1e38aeb6360b1af4 },
  1057. { 0xe2280b6c20dd5232, 0x25c6da63c38de1b1 },
  1058. { 0x8d590723948a535f, 0x579c487e5a38ad0f },
  1059. { 0xb0af48ec79ace837, 0x2d835a9df0c6d852 },
  1060. { 0xdcdb1b2798182244, 0xf8e431456cf88e66 },
  1061. { 0x8a08f0f8bf0f156b, 0x1b8e9ecb641b5900 },
  1062. { 0xac8b2d36eed2dac5, 0xe272467e3d222f40 },
  1063. { 0xd7adf884aa879177, 0x5b0ed81dcc6abb10 },
  1064. { 0x86ccbb52ea94baea, 0x98e947129fc2b4ea },
  1065. { 0xa87fea27a539e9a5, 0x3f2398d747b36225 },
  1066. { 0xd29fe4b18e88640e, 0x8eec7f0d19a03aae },
  1067. { 0x83a3eeeef9153e89, 0x1953cf68300424ad },
  1068. { 0xa48ceaaab75a8e2b, 0x5fa8c3423c052dd8 },
  1069. { 0xcdb02555653131b6, 0x3792f412cb06794e },
  1070. { 0x808e17555f3ebf11, 0xe2bbd88bbee40bd1 },
  1071. { 0xa0b19d2ab70e6ed6, 0x5b6aceaeae9d0ec5 },
  1072. { 0xc8de047564d20a8b, 0xf245825a5a445276 },
  1073. { 0xfb158592be068d2e, 0xeed6e2f0f0d56713 },
  1074. { 0x9ced737bb6c4183d, 0x55464dd69685606c },
  1075. { 0xc428d05aa4751e4c, 0xaa97e14c3c26b887 },
  1076. { 0xf53304714d9265df, 0xd53dd99f4b3066a9 },
  1077. { 0x993fe2c6d07b7fab, 0xe546a8038efe402a },
  1078. { 0xbf8fdb78849a5f96, 0xde98520472bdd034 },
  1079. { 0xef73d256a5c0f77c, 0x963e66858f6d4441 },
  1080. { 0x95a8637627989aad, 0xdde7001379a44aa9 },
  1081. { 0xbb127c53b17ec159, 0x5560c018580d5d53 },
  1082. { 0xe9d71b689dde71af, 0xaab8f01e6e10b4a7 },
  1083. { 0x9226712162ab070d, 0xcab3961304ca70e9 },
  1084. { 0xb6b00d69bb55c8d1, 0x3d607b97c5fd0d23 },
  1085. { 0xe45c10c42a2b3b05, 0x8cb89a7db77c506b },
  1086. { 0x8eb98a7a9a5b04e3, 0x77f3608e92adb243 },
  1087. { 0xb267ed1940f1c61c, 0x55f038b237591ed4 },
  1088. { 0xdf01e85f912e37a3, 0x6b6c46dec52f6689 },
  1089. { 0x8b61313bbabce2c6, 0x2323ac4b3b3da016 },
  1090. { 0xae397d8aa96c1b77, 0xabec975e0a0d081b },
  1091. { 0xd9c7dced53c72255, 0x96e7bd358c904a22 },
  1092. { 0x881cea14545c7575, 0x7e50d64177da2e55 },
  1093. { 0xaa242499697392d2, 0xdde50bd1d5d0b9ea },
  1094. { 0xd4ad2dbfc3d07787, 0x955e4ec64b44e865 },
  1095. { 0x84ec3c97da624ab4, 0xbd5af13bef0b113f },
  1096. { 0xa6274bbdd0fadd61, 0xecb1ad8aeacdd58f },
  1097. { 0xcfb11ead453994ba, 0x67de18eda5814af3 },
  1098. { 0x81ceb32c4b43fcf4, 0x80eacf948770ced8 },
  1099. { 0xa2425ff75e14fc31, 0xa1258379a94d028e },
  1100. { 0xcad2f7f5359a3b3e, 0x096ee45813a04331 },
  1101. { 0xfd87b5f28300ca0d, 0x8bca9d6e188853fd },
  1102. { 0x9e74d1b791e07e48, 0x775ea264cf55347e },
  1103. { 0xc612062576589dda, 0x95364afe032a819e },
  1104. { 0xf79687aed3eec551, 0x3a83ddbd83f52205 },
  1105. { 0x9abe14cd44753b52, 0xc4926a9672793543 },
  1106. { 0xc16d9a0095928a27, 0x75b7053c0f178294 },
  1107. { 0xf1c90080baf72cb1, 0x5324c68b12dd6339 },
  1108. { 0x971da05074da7bee, 0xd3f6fc16ebca5e04 },
  1109. { 0xbce5086492111aea, 0x88f4bb1ca6bcf585 },
  1110. { 0xec1e4a7db69561a5, 0x2b31e9e3d06c32e6 },
  1111. { 0x9392ee8e921d5d07, 0x3aff322e62439fd0 },
  1112. { 0xb877aa3236a4b449, 0x09befeb9fad487c3 },
  1113. { 0xe69594bec44de15b, 0x4c2ebe687989a9b4 },
  1114. { 0x901d7cf73ab0acd9, 0x0f9d37014bf60a11 },
  1115. { 0xb424dc35095cd80f, 0x538484c19ef38c95 },
  1116. { 0xe12e13424bb40e13, 0x2865a5f206b06fba },
  1117. { 0x8cbccc096f5088cb, 0xf93f87b7442e45d4 },
  1118. { 0xafebff0bcb24aafe, 0xf78f69a51539d749 },
  1119. { 0xdbe6fecebdedd5be, 0xb573440e5a884d1c },
  1120. { 0x89705f4136b4a597, 0x31680a88f8953031 },
  1121. { 0xabcc77118461cefc, 0xfdc20d2b36ba7c3e },
  1122. { 0xd6bf94d5e57a42bc, 0x3d32907604691b4d },
  1123. { 0x8637bd05af6c69b5, 0xa63f9a49c2c1b110 },
  1124. { 0xa7c5ac471b478423, 0x0fcf80dc33721d54 },
  1125. { 0xd1b71758e219652b, 0xd3c36113404ea4a9 },
  1126. { 0x83126e978d4fdf3b, 0x645a1cac083126ea },
  1127. { 0xa3d70a3d70a3d70a, 0x3d70a3d70a3d70a4 },
  1128. { 0xcccccccccccccccc, 0xcccccccccccccccd },
  1129. { 0x8000000000000000, 0x0000000000000000 },
  1130. { 0xa000000000000000, 0x0000000000000000 },
  1131. { 0xc800000000000000, 0x0000000000000000 },
  1132. { 0xfa00000000000000, 0x0000000000000000 },
  1133. { 0x9c40000000000000, 0x0000000000000000 },
  1134. { 0xc350000000000000, 0x0000000000000000 },
  1135. { 0xf424000000000000, 0x0000000000000000 },
  1136. { 0x9896800000000000, 0x0000000000000000 },
  1137. { 0xbebc200000000000, 0x0000000000000000 },
  1138. { 0xee6b280000000000, 0x0000000000000000 },
  1139. { 0x9502f90000000000, 0x0000000000000000 },
  1140. { 0xba43b74000000000, 0x0000000000000000 },
  1141. { 0xe8d4a51000000000, 0x0000000000000000 },
  1142. { 0x9184e72a00000000, 0x0000000000000000 },
  1143. { 0xb5e620f480000000, 0x0000000000000000 },
  1144. { 0xe35fa931a0000000, 0x0000000000000000 },
  1145. { 0x8e1bc9bf04000000, 0x0000000000000000 },
  1146. { 0xb1a2bc2ec5000000, 0x0000000000000000 },
  1147. { 0xde0b6b3a76400000, 0x0000000000000000 },
  1148. { 0x8ac7230489e80000, 0x0000000000000000 },
  1149. { 0xad78ebc5ac620000, 0x0000000000000000 },
  1150. { 0xd8d726b7177a8000, 0x0000000000000000 },
  1151. { 0x878678326eac9000, 0x0000000000000000 },
  1152. { 0xa968163f0a57b400, 0x0000000000000000 },
  1153. { 0xd3c21bcecceda100, 0x0000000000000000 },
  1154. { 0x84595161401484a0, 0x0000000000000000 },
  1155. { 0xa56fa5b99019a5c8, 0x0000000000000000 },
  1156. { 0xcecb8f27f4200f3a, 0x0000000000000000 },
  1157. { 0x813f3978f8940984, 0x4000000000000000 },
  1158. { 0xa18f07d736b90be5, 0x5000000000000000 },
  1159. { 0xc9f2c9cd04674ede, 0xa400000000000000 },
  1160. { 0xfc6f7c4045812296, 0x4d00000000000000 },
  1161. { 0x9dc5ada82b70b59d, 0xf020000000000000 },
  1162. { 0xc5371912364ce305, 0x6c28000000000000 },
  1163. { 0xf684df56c3e01bc6, 0xc732000000000000 },
  1164. { 0x9a130b963a6c115c, 0x3c7f400000000000 },
  1165. { 0xc097ce7bc90715b3, 0x4b9f100000000000 },
  1166. { 0xf0bdc21abb48db20, 0x1e86d40000000000 },
  1167. { 0x96769950b50d88f4, 0x1314448000000000 },
  1168. { 0xbc143fa4e250eb31, 0x17d955a000000000 },
  1169. { 0xeb194f8e1ae525fd, 0x5dcfab0800000000 },
  1170. { 0x92efd1b8d0cf37be, 0x5aa1cae500000000 },
  1171. { 0xb7abc627050305ad, 0xf14a3d9e40000000 },
  1172. { 0xe596b7b0c643c719, 0x6d9ccd05d0000000 },
  1173. { 0x8f7e32ce7bea5c6f, 0xe4820023a2000000 },
  1174. { 0xb35dbf821ae4f38b, 0xdda2802c8a800000 },
  1175. { 0xe0352f62a19e306e, 0xd50b2037ad200000 },
  1176. { 0x8c213d9da502de45, 0x4526f422cc340000 },
  1177. { 0xaf298d050e4395d6, 0x9670b12b7f410000 },
  1178. { 0xdaf3f04651d47b4c, 0x3c0cdd765f114000 },
  1179. { 0x88d8762bf324cd0f, 0xa5880a69fb6ac800 },
  1180. { 0xab0e93b6efee0053, 0x8eea0d047a457a00 },
  1181. { 0xd5d238a4abe98068, 0x72a4904598d6d880 },
  1182. { 0x85a36366eb71f041, 0x47a6da2b7f864750 },
  1183. { 0xa70c3c40a64e6c51, 0x999090b65f67d924 },
  1184. { 0xd0cf4b50cfe20765, 0xfff4b4e3f741cf6d },
  1185. { 0x82818f1281ed449f, 0xbff8f10e7a8921a4 },
  1186. { 0xa321f2d7226895c7, 0xaff72d52192b6a0d },
  1187. { 0xcbea6f8ceb02bb39, 0x9bf4f8a69f764490 },
  1188. { 0xfee50b7025c36a08, 0x02f236d04753d5b4 },
  1189. { 0x9f4f2726179a2245, 0x01d762422c946590 },
  1190. { 0xc722f0ef9d80aad6, 0x424d3ad2b7b97ef5 },
  1191. { 0xf8ebad2b84e0d58b, 0xd2e0898765a7deb2 },
  1192. { 0x9b934c3b330c8577, 0x63cc55f49f88eb2f },
  1193. { 0xc2781f49ffcfa6d5, 0x3cbf6b71c76b25fb },
  1194. { 0xf316271c7fc3908a, 0x8bef464e3945ef7a },
  1195. { 0x97edd871cfda3a56, 0x97758bf0e3cbb5ac },
  1196. { 0xbde94e8e43d0c8ec, 0x3d52eeed1cbea317 },
  1197. { 0xed63a231d4c4fb27, 0x4ca7aaa863ee4bdd },
  1198. { 0x945e455f24fb1cf8, 0x8fe8caa93e74ef6a },
  1199. { 0xb975d6b6ee39e436, 0xb3e2fd538e122b44 },
  1200. { 0xe7d34c64a9c85d44, 0x60dbbca87196b616 },
  1201. { 0x90e40fbeea1d3a4a, 0xbc8955e946fe31cd },
  1202. { 0xb51d13aea4a488dd, 0x6babab6398bdbe41 },
  1203. { 0xe264589a4dcdab14, 0xc696963c7eed2dd1 },
  1204. { 0x8d7eb76070a08aec, 0xfc1e1de5cf543ca2 },
  1205. { 0xb0de65388cc8ada8, 0x3b25a55f43294bcb },
  1206. { 0xdd15fe86affad912, 0x49ef0eb713f39ebe },
  1207. { 0x8a2dbf142dfcc7ab, 0x6e3569326c784337 },
  1208. { 0xacb92ed9397bf996, 0x49c2c37f07965404 },
  1209. { 0xd7e77a8f87daf7fb, 0xdc33745ec97be906 },
  1210. { 0x86f0ac99b4e8dafd, 0x69a028bb3ded71a3 },
  1211. { 0xa8acd7c0222311bc, 0xc40832ea0d68ce0c },
  1212. { 0xd2d80db02aabd62b, 0xf50a3fa490c30190 },
  1213. { 0x83c7088e1aab65db, 0x792667c6da79e0fa },
  1214. { 0xa4b8cab1a1563f52, 0x577001b891185938 },
  1215. { 0xcde6fd5e09abcf26, 0xed4c0226b55e6f86 },
  1216. { 0x80b05e5ac60b6178, 0x544f8158315b05b4 },
  1217. { 0xa0dc75f1778e39d6, 0x696361ae3db1c721 },
  1218. { 0xc913936dd571c84c, 0x03bc3a19cd1e38e9 },
  1219. { 0xfb5878494ace3a5f, 0x04ab48a04065c723 },
  1220. { 0x9d174b2dcec0e47b, 0x62eb0d64283f9c76 },
  1221. { 0xc45d1df942711d9a, 0x3ba5d0bd324f8394 },
  1222. { 0xf5746577930d6500, 0xca8f44ec7ee36479 },
  1223. { 0x9968bf6abbe85f20, 0x7e998b13cf4e1ecb },
  1224. { 0xbfc2ef456ae276e8, 0x9e3fedd8c321a67e },
  1225. { 0xefb3ab16c59b14a2, 0xc5cfe94ef3ea101e },
  1226. { 0x95d04aee3b80ece5, 0xbba1f1d158724a12 },
  1227. { 0xbb445da9ca61281f, 0x2a8a6e45ae8edc97 },
  1228. { 0xea1575143cf97226, 0xf52d09d71a3293bd },
  1229. { 0x924d692ca61be758, 0x593c2626705f9c56 },
  1230. { 0xb6e0c377cfa2e12e, 0x6f8b2fb00c77836c },
  1231. { 0xe498f455c38b997a, 0x0b6dfb9c0f956447 },
  1232. { 0x8edf98b59a373fec, 0x4724bd4189bd5eac },
  1233. { 0xb2977ee300c50fe7, 0x58edec91ec2cb657 },
  1234. { 0xdf3d5e9bc0f653e1, 0x2f2967b66737e3ed },
  1235. { 0x8b865b215899f46c, 0xbd79e0d20082ee74 },
  1236. { 0xae67f1e9aec07187, 0xecd8590680a3aa11 },
  1237. { 0xda01ee641a708de9, 0xe80e6f4820cc9495 },
  1238. { 0x884134fe908658b2, 0x3109058d147fdcdd },
  1239. { 0xaa51823e34a7eede, 0xbd4b46f0599fd415 },
  1240. { 0xd4e5e2cdc1d1ea96, 0x6c9e18ac7007c91a },
  1241. { 0x850fadc09923329e, 0x03e2cf6bc604ddb0 },
  1242. { 0xa6539930bf6bff45, 0x84db8346b786151c },
  1243. { 0xcfe87f7cef46ff16, 0xe612641865679a63 },
  1244. { 0x81f14fae158c5f6e, 0x4fcb7e8f3f60c07e },
  1245. { 0xa26da3999aef7749, 0xe3be5e330f38f09d },
  1246. { 0xcb090c8001ab551c, 0x5cadf5bfd3072cc5 },
  1247. { 0xfdcb4fa002162a63, 0x73d9732fc7c8f7f6 },
  1248. { 0x9e9f11c4014dda7e, 0x2867e7fddcdd9afa },
  1249. { 0xc646d63501a1511d, 0xb281e1fd541501b8 },
  1250. { 0xf7d88bc24209a565, 0x1f225a7ca91a4226 },
  1251. { 0x9ae757596946075f, 0x3375788de9b06958 },
  1252. { 0xc1a12d2fc3978937, 0x0052d6b1641c83ae },
  1253. { 0xf209787bb47d6b84, 0xc0678c5dbd23a49a },
  1254. { 0x9745eb4d50ce6332, 0xf840b7ba963646e0 },
  1255. { 0xbd176620a501fbff, 0xb650e5a93bc3d898 },
  1256. { 0xec5d3fa8ce427aff, 0xa3e51f138ab4cebe },
  1257. { 0x93ba47c980e98cdf, 0xc66f336c36b10137 },
  1258. { 0xb8a8d9bbe123f017, 0xb80b0047445d4184 },
  1259. { 0xe6d3102ad96cec1d, 0xa60dc059157491e5 },
  1260. { 0x9043ea1ac7e41392, 0x87c89837ad68db2f },
  1261. { 0xb454e4a179dd1877, 0x29babe4598c311fb },
  1262. { 0xe16a1dc9d8545e94, 0xf4296dd6fef3d67a },
  1263. { 0x8ce2529e2734bb1d, 0x1899e4a65f58660c },
  1264. { 0xb01ae745b101e9e4, 0x5ec05dcff72e7f8f },
  1265. { 0xdc21a1171d42645d, 0x76707543f4fa1f73 },
  1266. { 0x899504ae72497eba, 0x6a06494a791c53a8 },
  1267. { 0xabfa45da0edbde69, 0x0487db9d17636892 },
  1268. { 0xd6f8d7509292d603, 0x45a9d2845d3c42b6 },
  1269. { 0x865b86925b9bc5c2, 0x0b8a2392ba45a9b2 },
  1270. { 0xa7f26836f282b732, 0x8e6cac7768d7141e },
  1271. { 0xd1ef0244af2364ff, 0x3207d795430cd926 },
  1272. { 0x8335616aed761f1f, 0x7f44e6bd49e807b8 },
  1273. { 0xa402b9c5a8d3a6e7, 0x5f16206c9c6209a6 },
  1274. { 0xcd036837130890a1, 0x36dba887c37a8c0f },
  1275. { 0x802221226be55a64, 0xc2494954da2c9789 },
  1276. { 0xa02aa96b06deb0fd, 0xf2db9baa10b7bd6c },
  1277. { 0xc83553c5c8965d3d, 0x6f92829494e5acc7 },
  1278. { 0xfa42a8b73abbf48c, 0xcb772339ba1f17f9 },
  1279. { 0x9c69a97284b578d7, 0xff2a760414536efb },
  1280. { 0xc38413cf25e2d70d, 0xfef5138519684aba },
  1281. { 0xf46518c2ef5b8cd1, 0x7eb258665fc25d69 },
  1282. { 0x98bf2f79d5993802, 0xef2f773ffbd97a61 },
  1283. { 0xbeeefb584aff8603, 0xaafb550ffacfd8fa },
  1284. { 0xeeaaba2e5dbf6784, 0x95ba2a53f983cf38 },
  1285. { 0x952ab45cfa97a0b2, 0xdd945a747bf26183 },
  1286. { 0xba756174393d88df, 0x94f971119aeef9e4 },
  1287. { 0xe912b9d1478ceb17, 0x7a37cd5601aab85d },
  1288. { 0x91abb422ccb812ee, 0xac62e055c10ab33a },
  1289. { 0xb616a12b7fe617aa, 0x577b986b314d6009 },
  1290. { 0xe39c49765fdf9d94, 0xed5a7e85fda0b80b },
  1291. { 0x8e41ade9fbebc27d, 0x14588f13be847307 },
  1292. { 0xb1d219647ae6b31c, 0x596eb2d8ae258fc8 },
  1293. { 0xde469fbd99a05fe3, 0x6fca5f8ed9aef3bb },
  1294. { 0x8aec23d680043bee, 0x25de7bb9480d5854 },
  1295. { 0xada72ccc20054ae9, 0xaf561aa79a10ae6a },
  1296. { 0xd910f7ff28069da4, 0x1b2ba1518094da04 },
  1297. { 0x87aa9aff79042286, 0x90fb44d2f05d0842 },
  1298. { 0xa99541bf57452b28, 0x353a1607ac744a53 },
  1299. { 0xd3fa922f2d1675f2, 0x42889b8997915ce8 },
  1300. { 0x847c9b5d7c2e09b7, 0x69956135febada11 },
  1301. { 0xa59bc234db398c25, 0x43fab9837e699095 },
  1302. { 0xcf02b2c21207ef2e, 0x94f967e45e03f4bb },
  1303. { 0x8161afb94b44f57d, 0x1d1be0eebac278f5 },
  1304. { 0xa1ba1ba79e1632dc, 0x6462d92a69731732 },
  1305. { 0xca28a291859bbf93, 0x7d7b8f7503cfdcfe },
  1306. { 0xfcb2cb35e702af78, 0x5cda735244c3d43e },
  1307. { 0x9defbf01b061adab, 0x3a0888136afa64a7 },
  1308. { 0xc56baec21c7a1916, 0x088aaa1845b8fdd0 },
  1309. { 0xf6c69a72a3989f5b, 0x8aad549e57273d45 },
  1310. { 0x9a3c2087a63f6399, 0x36ac54e2f678864b },
  1311. { 0xc0cb28a98fcf3c7f, 0x84576a1bb416a7dd },
  1312. { 0xf0fdf2d3f3c30b9f, 0x656d44a2a11c51d5 },
  1313. { 0x969eb7c47859e743, 0x9f644ae5a4b1b325 },
  1314. { 0xbc4665b596706114, 0x873d5d9f0dde1fee },
  1315. { 0xeb57ff22fc0c7959, 0xa90cb506d155a7ea },
  1316. { 0x9316ff75dd87cbd8, 0x09a7f12442d588f2 },
  1317. { 0xb7dcbf5354e9bece, 0x0c11ed6d538aeb2f },
  1318. { 0xe5d3ef282a242e81, 0x8f1668c8a86da5fa },
  1319. { 0x8fa475791a569d10, 0xf96e017d694487bc },
  1320. { 0xb38d92d760ec4455, 0x37c981dcc395a9ac },
  1321. { 0xe070f78d3927556a, 0x85bbe253f47b1417 },
  1322. { 0x8c469ab843b89562, 0x93956d7478ccec8e },
  1323. { 0xaf58416654a6babb, 0x387ac8d1970027b2 },
  1324. { 0xdb2e51bfe9d0696a, 0x06997b05fcc0319e },
  1325. { 0x88fcf317f22241e2, 0x441fece3bdf81f03 },
  1326. { 0xab3c2fddeeaad25a, 0xd527e81cad7626c3 },
  1327. { 0xd60b3bd56a5586f1, 0x8a71e223d8d3b074 },
  1328. { 0x85c7056562757456, 0xf6872d5667844e49 },
  1329. { 0xa738c6bebb12d16c, 0xb428f8ac016561db },
  1330. { 0xd106f86e69d785c7, 0xe13336d701beba52 },
  1331. { 0x82a45b450226b39c, 0xecc0024661173473 },
  1332. { 0xa34d721642b06084, 0x27f002d7f95d0190 },
  1333. { 0xcc20ce9bd35c78a5, 0x31ec038df7b441f4 },
  1334. { 0xff290242c83396ce, 0x7e67047175a15271 },
  1335. { 0x9f79a169bd203e41, 0x0f0062c6e984d386 },
  1336. { 0xc75809c42c684dd1, 0x52c07b78a3e60868 },
  1337. { 0xf92e0c3537826145, 0xa7709a56ccdf8a82 },
  1338. { 0x9bbcc7a142b17ccb, 0x88a66076400bb691 },
  1339. { 0xc2abf989935ddbfe, 0x6acff893d00ea435 },
  1340. { 0xf356f7ebf83552fe, 0x0583f6b8c4124d43 },
  1341. { 0x98165af37b2153de, 0xc3727a337a8b704a },
  1342. { 0xbe1bf1b059e9a8d6, 0x744f18c0592e4c5c },
  1343. { 0xeda2ee1c7064130c, 0x1162def06f79df73 },
  1344. { 0x9485d4d1c63e8be7, 0x8addcb5645ac2ba8 },
  1345. { 0xb9a74a0637ce2ee1, 0x6d953e2bd7173692 },
  1346. { 0xe8111c87c5c1ba99, 0xc8fa8db6ccdd0437 },
  1347. { 0x910ab1d4db9914a0, 0x1d9c9892400a22a2 },
  1348. { 0xb54d5e4a127f59c8, 0x2503beb6d00cab4b },
  1349. { 0xe2a0b5dc971f303a, 0x2e44ae64840fd61d },
  1350. { 0x8da471a9de737e24, 0x5ceaecfed289e5d2 },
  1351. { 0xb10d8e1456105dad, 0x7425a83e872c5f47 },
  1352. { 0xdd50f1996b947518, 0xd12f124e28f77719 },
  1353. { 0x8a5296ffe33cc92f, 0x82bd6b70d99aaa6f },
  1354. { 0xace73cbfdc0bfb7b, 0x636cc64d1001550b },
  1355. { 0xd8210befd30efa5a, 0x3c47f7e05401aa4e },
  1356. { 0x8714a775e3e95c78, 0x65acfaec34810a71 },
  1357. { 0xa8d9d1535ce3b396, 0x7f1839a741a14d0d },
  1358. { 0xd31045a8341ca07c, 0x1ede48111209a050 },
  1359. { 0x83ea2b892091e44d, 0x934aed0aab460432 },
  1360. { 0xa4e4b66b68b65d60, 0xf81da84d5617853f },
  1361. { 0xce1de40642e3f4b9, 0x36251260ab9d668e },
  1362. { 0x80d2ae83e9ce78f3, 0xc1d72b7c6b426019 },
  1363. { 0xa1075a24e4421730, 0xb24cf65b8612f81f },
  1364. { 0xc94930ae1d529cfc, 0xdee033f26797b627 },
  1365. { 0xfb9b7cd9a4a7443c, 0x169840ef017da3b1 },
  1366. { 0x9d412e0806e88aa5, 0x8e1f289560ee864e },
  1367. { 0xc491798a08a2ad4e, 0xf1a6f2bab92a27e2 },
  1368. { 0xf5b5d7ec8acb58a2, 0xae10af696774b1db },
  1369. { 0x9991a6f3d6bf1765, 0xacca6da1e0a8ef29 },
  1370. { 0xbff610b0cc6edd3f, 0x17fd090a58d32af3 },
  1371. { 0xeff394dcff8a948e, 0xddfc4b4cef07f5b0 },
  1372. { 0x95f83d0a1fb69cd9, 0x4abdaf101564f98e },
  1373. { 0xbb764c4ca7a4440f, 0x9d6d1ad41abe37f1 },
  1374. { 0xea53df5fd18d5513, 0x84c86189216dc5ed },
  1375. { 0x92746b9be2f8552c, 0x32fd3cf5b4e49bb4 },
  1376. { 0xb7118682dbb66a77, 0x3fbc8c33221dc2a1 },
  1377. { 0xe4d5e82392a40515, 0x0fabaf3feaa5334a },
  1378. { 0x8f05b1163ba6832d, 0x29cb4d87f2a7400e },
  1379. { 0xb2c71d5bca9023f8, 0x743e20e9ef511012 },
  1380. { 0xdf78e4b2bd342cf6, 0x914da9246b255416 },
  1381. { 0x8bab8eefb6409c1a, 0x1ad089b6c2f7548e },
  1382. { 0xae9672aba3d0c320, 0xa184ac2473b529b1 },
  1383. { 0xda3c0f568cc4f3e8, 0xc9e5d72d90a2741e },
  1384. { 0x8865899617fb1871, 0x7e2fa67c7a658892 },
  1385. { 0xaa7eebfb9df9de8d, 0xddbb901b98feeab7 },
  1386. { 0xd51ea6fa85785631, 0x552a74227f3ea565 },
  1387. { 0x8533285c936b35de, 0xd53a88958f87275f },
  1388. { 0xa67ff273b8460356, 0x8a892abaf368f137 },
  1389. { 0xd01fef10a657842c, 0x2d2b7569b0432d85 },
  1390. { 0x8213f56a67f6b29b, 0x9c3b29620e29fc73 },
  1391. { 0xa298f2c501f45f42, 0x8349f3ba91b47b8f },
  1392. { 0xcb3f2f7642717713, 0x241c70a936219a73 },
  1393. { 0xfe0efb53d30dd4d7, 0xed238cd383aa0110 },
  1394. { 0x9ec95d1463e8a506, 0xf4363804324a40aa },
  1395. { 0xc67bb4597ce2ce48, 0xb143c6053edcd0d5 },
  1396. { 0xf81aa16fdc1b81da, 0xdd94b7868e94050a },
  1397. { 0x9b10a4e5e9913128, 0xca7cf2b4191c8326 },
  1398. { 0xc1d4ce1f63f57d72, 0xfd1c2f611f63a3f0 },
  1399. { 0xf24a01a73cf2dccf, 0xbc633b39673c8cec },
  1400. { 0x976e41088617ca01, 0xd5be0503e085d813 },
  1401. { 0xbd49d14aa79dbc82, 0x4b2d8644d8a74e18 },
  1402. { 0xec9c459d51852ba2, 0xddf8e7d60ed1219e },
  1403. { 0x93e1ab8252f33b45, 0xcabb90e5c942b503 },
  1404. { 0xb8da1662e7b00a17, 0x3d6a751f3b936243 },
  1405. { 0xe7109bfba19c0c9d, 0x0cc512670a783ad4 },
  1406. { 0x906a617d450187e2, 0x27fb2b80668b24c5 },
  1407. { 0xb484f9dc9641e9da, 0xb1f9f660802dedf6 },
  1408. { 0xe1a63853bbd26451, 0x5e7873f8a0396973 },
  1409. { 0x8d07e33455637eb2, 0xdb0b487b6423e1e8 },
  1410. { 0xb049dc016abc5e5f, 0x91ce1a9a3d2cda62 },
  1411. { 0xdc5c5301c56b75f7, 0x7641a140cc7810fb },
  1412. { 0x89b9b3e11b6329ba, 0xa9e904c87fcb0a9d },
  1413. { 0xac2820d9623bf429, 0x546345fa9fbdcd44 },
  1414. { 0xd732290fbacaf133, 0xa97c177947ad4095 },
  1415. { 0x867f59a9d4bed6c0, 0x49ed8eabcccc485d },
  1416. { 0xa81f301449ee8c70, 0x5c68f256bfff5a74 },
  1417. { 0xd226fc195c6a2f8c, 0x73832eec6fff3111 },
  1418. { 0x83585d8fd9c25db7, 0xc831fd53c5ff7eab },
  1419. { 0xa42e74f3d032f525, 0xba3e7ca8b77f5e55 },
  1420. { 0xcd3a1230c43fb26f, 0x28ce1bd2e55f35eb },
  1421. { 0x80444b5e7aa7cf85, 0x7980d163cf5b81b3 },
  1422. { 0xa0555e361951c366, 0xd7e105bcc332621f },
  1423. { 0xc86ab5c39fa63440, 0x8dd9472bf3fefaa7 },
  1424. { 0xfa856334878fc150, 0xb14f98f6f0feb951 },
  1425. { 0x9c935e00d4b9d8d2, 0x6ed1bf9a569f33d3 },
  1426. { 0xc3b8358109e84f07, 0x0a862f80ec4700c8 },
  1427. { 0xf4a642e14c6262c8, 0xcd27bb612758c0fa },
  1428. { 0x98e7e9cccfbd7dbd, 0x8038d51cb897789c },
  1429. { 0xbf21e44003acdd2c, 0xe0470a63e6bd56c3 },
  1430. { 0xeeea5d5004981478, 0x1858ccfce06cac74 },
  1431. { 0x95527a5202df0ccb, 0x0f37801e0c43ebc8 },
  1432. { 0xbaa718e68396cffd, 0xd30560258f54e6ba },
  1433. { 0xe950df20247c83fd, 0x47c6b82ef32a2069 },
  1434. { 0x91d28b7416cdd27e, 0x4cdc331d57fa5441 },
  1435. { 0xb6472e511c81471d, 0xe0133fe4adf8e952 },
  1436. { 0xe3d8f9e563a198e5, 0x58180fddd97723a6 },
  1437. { 0x8e679c2f5e44ff8f, 0x570f09eaa7ea7648 },
  1438. { 0xb201833b35d63f73, 0x2cd2cc6551e513da },
  1439. { 0xde81e40a034bcf4f, 0xf8077f7ea65e58d1 },
  1440. { 0x8b112e86420f6191, 0xfb04afaf27faf782 },
  1441. { 0xadd57a27d29339f6, 0x79c5db9af1f9b563 },
  1442. { 0xd94ad8b1c7380874, 0x18375281ae7822bc },
  1443. { 0x87cec76f1c830548, 0x8f2293910d0b15b5 },
  1444. { 0xa9c2794ae3a3c69a, 0xb2eb3875504ddb22 },
  1445. { 0xd433179d9c8cb841, 0x5fa60692a46151eb },
  1446. { 0x849feec281d7f328, 0xdbc7c41ba6bcd333 },
  1447. { 0xa5c7ea73224deff3, 0x12b9b522906c0800 },
  1448. { 0xcf39e50feae16bef, 0xd768226b34870a00 },
  1449. { 0x81842f29f2cce375, 0xe6a1158300d46640 },
  1450. { 0xa1e53af46f801c53, 0x60495ae3c1097fd0 },
  1451. { 0xca5e89b18b602368, 0x385bb19cb14bdfc4 },
  1452. { 0xfcf62c1dee382c42, 0x46729e03dd9ed7b5 },
  1453. { 0x9e19db92b4e31ba9, 0x6c07a2c26a8346d1 },
  1454. { 0xc5a05277621be293, 0xc7098b7305241885 },
  1455. { 0xf70867153aa2db38, 0xb8cbee4fc66d1ea7 }
  1456. };
  1457. };
  1458. // Compressed cache for double
  1459. struct compressed_cache_detail {
  1460. static constexpr int compression_ratio = 27;
  1461. static constexpr std::size_t compressed_table_size =
  1462. (cache_holder<ieee754_format::binary64>::max_k -
  1463. cache_holder<ieee754_format::binary64>::min_k + compression_ratio) / compression_ratio;
  1464. struct cache_holder_t {
  1465. wuint::uint128 table[compressed_table_size];
  1466. };
  1467. static constexpr cache_holder_t cache = [] {
  1468. cache_holder_t res{};
  1469. for (std::size_t i = 0; i < compressed_table_size; ++i) {
  1470. res.table[i] = cache_holder<ieee754_format::binary64>::cache[i * compression_ratio];
  1471. }
  1472. return res;
  1473. }();
  1474. struct pow5_holder_t {
  1475. std::uint64_t table[compression_ratio];
  1476. };
  1477. static constexpr pow5_holder_t pow5 = [] {
  1478. pow5_holder_t res{};
  1479. std::uint64_t p = 1;
  1480. for (std::size_t i = 0; i < compression_ratio; ++i) {
  1481. res.table[i] = p;
  1482. p *= 5;
  1483. }
  1484. return res;
  1485. }();
  1486. static constexpr std::uint32_t errors[] = {
  1487. 0x50001400, 0x54044100, 0x54014555, 0x55954415, 0x54115555,
  1488. 0x00000001, 0x50000000, 0x00104000, 0x54010004, 0x05004001,
  1489. 0x55555544, 0x41545555, 0x54040551, 0x15445545, 0x51555514,
  1490. 0x10000015, 0x00101100, 0x01100015, 0x00000000, 0x00000000,
  1491. 0x00000000, 0x00000000, 0x04450514, 0x45414110, 0x55555145,
  1492. 0x50544050, 0x15040155, 0x11054140, 0x50111514, 0x11451454,
  1493. 0x00400541, 0x00000000, 0x55555450, 0x10056551, 0x10054011,
  1494. 0x55551014, 0x69514555, 0x05151109, 0x00155555
  1495. };
  1496. };
  1497. }
  1498. ////////////////////////////////////////////////////////////////////////////////////////
  1499. // Policies
  1500. ////////////////////////////////////////////////////////////////////////////////////////
  1501. namespace detail {
  1502. // Forward declare the implementation class
  1503. template <class Float>
  1504. struct impl;
  1505. namespace policy_impl {
  1506. // Sign policy
  1507. namespace sign {
  1508. struct base {};
  1509. struct ignore : base {
  1510. using sign_policy = ignore;
  1511. static constexpr bool return_has_sign = false;
  1512. template <class Float, class Fp>
  1513. static constexpr void handle_sign(ieee754_bits<Float>, Fp&) noexcept {}
  1514. };
  1515. struct return_sign : base {
  1516. using sign_policy = return_sign;
  1517. static constexpr bool return_has_sign = true;
  1518. template <class Float, class Fp>
  1519. static constexpr void handle_sign(ieee754_bits<Float> br, Fp& fp) noexcept {
  1520. fp.is_negative = br.is_negative();
  1521. }
  1522. };
  1523. }
  1524. // Trailing zero policy
  1525. namespace trailing_zero {
  1526. struct base {};
  1527. struct ignore : base {
  1528. using trailing_zero_policy = ignore;
  1529. static constexpr bool report_trailing_zeros = false;
  1530. template <class Fp>
  1531. static constexpr void on_trailing_zeros(Fp&) noexcept {}
  1532. template <class Fp>
  1533. static constexpr void no_trailing_zeros(Fp&) noexcept {}
  1534. };
  1535. struct remove : base {
  1536. using trailing_zero_policy = remove;
  1537. static constexpr bool report_trailing_zeros = false;
  1538. template <class Fp>
  1539. static constexpr void on_trailing_zeros(Fp& fp) noexcept {
  1540. fp.exponent +=
  1541. impl<typename Fp::float_type>::remove_trailing_zeros(fp.significand);
  1542. }
  1543. template <class Fp>
  1544. static constexpr void no_trailing_zeros(Fp&) noexcept {}
  1545. };
  1546. struct report : base {
  1547. using trailing_zero_policy = report;
  1548. static constexpr bool report_trailing_zeros = true;
  1549. template <class Fp>
  1550. static constexpr void on_trailing_zeros(Fp& fp) noexcept {
  1551. fp.may_have_trailing_zeros = true;
  1552. }
  1553. template <class Fp>
  1554. static constexpr void no_trailing_zeros(Fp& fp) noexcept {
  1555. fp.may_have_trailing_zeros = false;
  1556. }
  1557. };
  1558. }
  1559. // Rounding mode policy
  1560. namespace rounding_mode {
  1561. struct base {};
  1562. enum class tag_t {
  1563. to_nearest,
  1564. left_closed_directed,
  1565. right_closed_directed
  1566. };
  1567. namespace interval_type {
  1568. struct symmetric_boundary {
  1569. static constexpr bool is_symmetric = true;
  1570. bool is_closed;
  1571. constexpr bool include_left_endpoint() const noexcept {
  1572. return is_closed;
  1573. }
  1574. constexpr bool include_right_endpoint() const noexcept {
  1575. return is_closed;
  1576. }
  1577. };
  1578. struct asymmetric_boundary {
  1579. static constexpr bool is_symmetric = false;
  1580. bool is_left_closed;
  1581. constexpr bool include_left_endpoint() const noexcept {
  1582. return is_left_closed;
  1583. }
  1584. constexpr bool include_right_endpoint() const noexcept {
  1585. return !is_left_closed;
  1586. }
  1587. };
  1588. struct closed {
  1589. static constexpr bool is_symmetric = true;
  1590. static constexpr bool include_left_endpoint() noexcept {
  1591. return true;
  1592. }
  1593. static constexpr bool include_right_endpoint() noexcept {
  1594. return true;
  1595. }
  1596. };
  1597. struct open {
  1598. static constexpr bool is_symmetric = true;
  1599. static constexpr bool include_left_endpoint() noexcept {
  1600. return false;
  1601. }
  1602. static constexpr bool include_right_endpoint() noexcept {
  1603. return false;
  1604. }
  1605. };
  1606. struct left_closed_right_open {
  1607. static constexpr bool is_symmetric = false;
  1608. static constexpr bool include_left_endpoint() noexcept {
  1609. return true;
  1610. }
  1611. static constexpr bool include_right_endpoint() noexcept {
  1612. return false;
  1613. }
  1614. };
  1615. struct right_closed_left_open {
  1616. static constexpr bool is_symmetric = false;
  1617. static constexpr bool include_left_endpoint() noexcept {
  1618. return false;
  1619. }
  1620. static constexpr bool include_right_endpoint() noexcept {
  1621. return true;
  1622. }
  1623. };
  1624. }
  1625. struct nearest_to_even : base {
  1626. using rounding_mode_policy = nearest_to_even;
  1627. static constexpr auto tag = tag_t::to_nearest;
  1628. template <class Float, class Func>
  1629. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1630. return f(nearest_to_even{});
  1631. }
  1632. template <class Float>
  1633. static constexpr interval_type::symmetric_boundary
  1634. interval_type_normal(ieee754_bits<Float> br) noexcept
  1635. {
  1636. return{ br.u % 2 == 0 };
  1637. }
  1638. template <class Float>
  1639. static constexpr interval_type::closed
  1640. interval_type_shorter(ieee754_bits<Float>) noexcept
  1641. {
  1642. return{};
  1643. }
  1644. };
  1645. struct nearest_to_odd : base {
  1646. using rounding_mode_policy = nearest_to_odd;
  1647. static constexpr auto tag = tag_t::to_nearest;
  1648. template <class Float, class Func>
  1649. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1650. return f(nearest_to_odd{});
  1651. }
  1652. template <class Float>
  1653. static constexpr interval_type::symmetric_boundary
  1654. interval_type_normal(ieee754_bits<Float> br) noexcept
  1655. {
  1656. return{ br.u % 2 != 0 };
  1657. }
  1658. template <class Float>
  1659. static constexpr interval_type::closed
  1660. interval_type_shorter(ieee754_bits<Float>) noexcept
  1661. {
  1662. return{};
  1663. }
  1664. };
  1665. struct nearest_toward_plus_infinity : base {
  1666. using rounding_mode_policy = nearest_toward_plus_infinity;
  1667. static constexpr auto tag = tag_t::to_nearest;
  1668. template <class Float, class Func>
  1669. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1670. return f(nearest_toward_plus_infinity{});
  1671. }
  1672. template <class Float>
  1673. static constexpr interval_type::asymmetric_boundary
  1674. interval_type_normal(ieee754_bits<Float> br) noexcept
  1675. {
  1676. return{ !br.is_negative() };
  1677. }
  1678. template <class Float>
  1679. static constexpr interval_type::asymmetric_boundary
  1680. interval_type_shorter(ieee754_bits<Float> br) noexcept
  1681. {
  1682. return{ !br.is_negative() };
  1683. }
  1684. };
  1685. struct nearest_toward_minus_infinity : base {
  1686. using rounding_mode_policy = nearest_toward_minus_infinity;
  1687. static constexpr auto tag = tag_t::to_nearest;
  1688. template <class Float, class Func>
  1689. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1690. return f(nearest_toward_minus_infinity{});
  1691. }
  1692. template <class Float>
  1693. static constexpr interval_type::asymmetric_boundary
  1694. interval_type_normal(ieee754_bits<Float> br) noexcept
  1695. {
  1696. return{ br.is_negative() };
  1697. }
  1698. template <class Float>
  1699. static constexpr interval_type::asymmetric_boundary
  1700. interval_type_shorter(ieee754_bits<Float> br) noexcept
  1701. {
  1702. return{ br.is_negative() };
  1703. }
  1704. };
  1705. struct nearest_toward_zero : base {
  1706. using rounding_mode_policy = nearest_toward_zero;
  1707. static constexpr auto tag = tag_t::to_nearest;
  1708. template <class Float, class Func>
  1709. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1710. return f(nearest_toward_zero{});
  1711. }
  1712. template <class Float>
  1713. static constexpr interval_type::right_closed_left_open
  1714. interval_type_normal(ieee754_bits<Float>) noexcept
  1715. {
  1716. return{};
  1717. }
  1718. template <class Float>
  1719. static constexpr interval_type::right_closed_left_open
  1720. interval_type_shorter(ieee754_bits<Float>) noexcept
  1721. {
  1722. return{};
  1723. }
  1724. };
  1725. struct nearest_away_from_zero : base {
  1726. using rounding_mode_policy = nearest_away_from_zero;
  1727. static constexpr auto tag = tag_t::to_nearest;
  1728. template <class Float, class Func>
  1729. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1730. return f(nearest_away_from_zero{});
  1731. }
  1732. template <class Float>
  1733. static constexpr interval_type::left_closed_right_open
  1734. interval_type_normal(ieee754_bits<Float>) noexcept
  1735. {
  1736. return{};
  1737. }
  1738. template <class Float>
  1739. static constexpr interval_type::left_closed_right_open
  1740. interval_type_shorter(ieee754_bits<Float>) noexcept
  1741. {
  1742. return{};
  1743. }
  1744. };
  1745. namespace detail {
  1746. struct nearest_always_closed {
  1747. static constexpr auto tag = tag_t::to_nearest;
  1748. template <class Float>
  1749. static constexpr interval_type::closed
  1750. interval_type_normal(ieee754_bits<Float>) noexcept
  1751. {
  1752. return{};
  1753. }
  1754. template <class Float>
  1755. static constexpr interval_type::closed
  1756. interval_type_shorter(ieee754_bits<Float>) noexcept
  1757. {
  1758. return{};
  1759. }
  1760. };
  1761. struct nearest_always_open {
  1762. static constexpr auto tag = tag_t::to_nearest;
  1763. template <class Float>
  1764. static constexpr interval_type::open
  1765. interval_type_normal(ieee754_bits<Float>) noexcept
  1766. {
  1767. return{};
  1768. }
  1769. template <class Float>
  1770. static constexpr interval_type::open
  1771. interval_type_shorter(ieee754_bits<Float>) noexcept
  1772. {
  1773. return{};
  1774. }
  1775. };
  1776. }
  1777. struct nearest_to_even_static_boundary : base {
  1778. using rounding_mode_policy = nearest_to_even_static_boundary;
  1779. template <class Float, class Func>
  1780. static auto delegate(ieee754_bits<Float> br, Func&& f) noexcept {
  1781. if (br.u % 2 == 0) {
  1782. return f(detail::nearest_always_closed{});
  1783. }
  1784. else {
  1785. return f(detail::nearest_always_open{});
  1786. }
  1787. }
  1788. };
  1789. struct nearest_to_odd_static_boundary : base {
  1790. using rounding_mode_policy = nearest_to_odd_static_boundary;
  1791. template <class Float, class Func>
  1792. static auto delegate(ieee754_bits<Float> br, Func&& f) noexcept {
  1793. if (br.u % 2 == 0) {
  1794. return f(detail::nearest_always_open{});
  1795. }
  1796. else {
  1797. return f(detail::nearest_always_closed{});
  1798. }
  1799. }
  1800. };
  1801. struct nearest_toward_plus_infinity_static_boundary : base {
  1802. using rounding_mode_policy = nearest_toward_plus_infinity_static_boundary;
  1803. template <class Float, class Func>
  1804. static auto delegate(ieee754_bits<Float> br, Func&& f) noexcept {
  1805. if (br.is_negative()) {
  1806. return f(nearest_toward_zero{});
  1807. }
  1808. else {
  1809. return f(nearest_away_from_zero{});
  1810. }
  1811. }
  1812. };
  1813. struct nearest_toward_minus_infinity_static_boundary : base {
  1814. using rounding_mode_policy = nearest_toward_minus_infinity_static_boundary;
  1815. template <class Float, class Func>
  1816. static auto delegate(ieee754_bits<Float> br, Func&& f) noexcept {
  1817. if (br.is_negative()) {
  1818. return f(nearest_away_from_zero{});
  1819. }
  1820. else {
  1821. return f(nearest_toward_zero{});
  1822. }
  1823. }
  1824. };
  1825. namespace detail {
  1826. struct left_closed_directed {
  1827. static constexpr auto tag = tag_t::left_closed_directed;
  1828. template <class Float>
  1829. static constexpr interval_type::left_closed_right_open
  1830. interval_type_normal(ieee754_bits<Float>) noexcept
  1831. {
  1832. return{};
  1833. }
  1834. };
  1835. struct right_closed_directed {
  1836. static constexpr auto tag = tag_t::right_closed_directed;
  1837. template <class Float>
  1838. static constexpr interval_type::right_closed_left_open
  1839. interval_type_normal(ieee754_bits<Float>) noexcept
  1840. {
  1841. return{};
  1842. }
  1843. };
  1844. }
  1845. struct toward_plus_infinity : base {
  1846. using rounding_mode_policy = toward_plus_infinity;
  1847. template <class Float, class Func>
  1848. static auto delegate(ieee754_bits<Float> br, Func&& f) noexcept {
  1849. if (br.is_negative()) {
  1850. return f(detail::left_closed_directed{});
  1851. }
  1852. else {
  1853. return f(detail::right_closed_directed{});
  1854. }
  1855. }
  1856. };
  1857. struct toward_minus_infinity : base {
  1858. using rounding_mode_policy = toward_minus_infinity;
  1859. template <class Float, class Func>
  1860. static auto delegate(ieee754_bits<Float> br, Func&& f) noexcept {
  1861. if (br.is_negative()) {
  1862. return f(detail::right_closed_directed{});
  1863. }
  1864. else {
  1865. return f(detail::left_closed_directed{});
  1866. }
  1867. }
  1868. };
  1869. struct toward_zero : base {
  1870. using rounding_mode_policy = toward_zero;
  1871. template <class Float, class Func>
  1872. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1873. return f(detail::left_closed_directed{});
  1874. }
  1875. };
  1876. struct away_from_zero : base {
  1877. using rounding_mode_policy = away_from_zero;
  1878. template <class Float, class Func>
  1879. static auto delegate(ieee754_bits<Float>, Func&& f) noexcept {
  1880. return f(detail::right_closed_directed{});
  1881. }
  1882. };
  1883. }
  1884. // Correct rounding policy
  1885. namespace correct_rounding {
  1886. struct base {};
  1887. enum class tag_t {
  1888. do_not_care,
  1889. to_even,
  1890. to_odd,
  1891. away_from_zero,
  1892. toward_zero
  1893. };
  1894. struct do_not_care : base {
  1895. using correct_rounding_policy = do_not_care;
  1896. static constexpr auto tag = tag_t::do_not_care;
  1897. template <class Fp>
  1898. static constexpr void break_rounding_tie(Fp&) noexcept {}
  1899. };
  1900. struct to_even : base {
  1901. using correct_rounding_policy = to_even;
  1902. static constexpr auto tag = tag_t::to_even;
  1903. template <class Fp>
  1904. static constexpr void break_rounding_tie(Fp& fp) noexcept
  1905. {
  1906. fp.significand = fp.significand % 2 == 0 ?
  1907. fp.significand : fp.significand - 1;
  1908. }
  1909. };
  1910. struct to_odd : base {
  1911. using correct_rounding_policy = to_odd;
  1912. static constexpr auto tag = tag_t::to_odd;
  1913. template <class Fp>
  1914. static constexpr void break_rounding_tie(Fp& fp) noexcept
  1915. {
  1916. fp.significand = fp.significand % 2 != 0 ?
  1917. fp.significand : fp.significand - 1;
  1918. }
  1919. };
  1920. struct away_from_zero : base {
  1921. using correct_rounding_policy = away_from_zero;
  1922. static constexpr auto tag = tag_t::away_from_zero;
  1923. template <class Fp>
  1924. static constexpr void break_rounding_tie(Fp& fp) noexcept {}
  1925. };
  1926. struct toward_zero : base {
  1927. using correct_rounding_policy = toward_zero;
  1928. static constexpr auto tag = tag_t::toward_zero;
  1929. template <class Fp>
  1930. static constexpr void break_rounding_tie(Fp& fp) noexcept
  1931. {
  1932. --fp.significand;
  1933. }
  1934. };
  1935. }
  1936. namespace cache {
  1937. struct base {};
  1938. struct normal : base {
  1939. using cache_policy = normal;
  1940. template <ieee754_format format>
  1941. static constexpr typename cache_holder<format>::cache_entry_type get_cache(int k) noexcept {
  1942. assert(k >= cache_holder<format>::min_k && k <= cache_holder<format>::max_k);
  1943. return cache_holder<format>::cache[std::size_t(k - cache_holder<format>::min_k)];
  1944. }
  1945. };
  1946. struct compressed : base {
  1947. using cache_policy = compressed;
  1948. template <ieee754_format format>
  1949. static constexpr typename cache_holder<format>::cache_entry_type get_cache(int k) noexcept {
  1950. assert(k >= cache_holder<format>::min_k && k <= cache_holder<format>::max_k);
  1951. if constexpr (format == ieee754_format::binary64)
  1952. {
  1953. // Compute base index
  1954. auto cache_index = (k - cache_holder<format>::min_k) /
  1955. compressed_cache_detail::compression_ratio;
  1956. auto kb = cache_index * compressed_cache_detail::compression_ratio
  1957. + cache_holder<format>::min_k;
  1958. auto offset = k - kb;
  1959. // Get base cache
  1960. auto base_cache = compressed_cache_detail::cache.table[cache_index];
  1961. if (offset == 0) {
  1962. return base_cache;
  1963. }
  1964. else {
  1965. // Compute the required amount of bit-shift
  1966. auto alpha = log::floor_log2_pow10(kb + offset) - log::floor_log2_pow10(kb) - offset;
  1967. assert(alpha > 0 && alpha < 64);
  1968. // Try to recover the real cache
  1969. auto pow5 = compressed_cache_detail::pow5.table[offset];
  1970. auto recovered_cache = wuint::umul128(base_cache.high(), pow5);
  1971. auto middle_low = wuint::umul128(base_cache.low() - (kb < 0 ? 1 : 0), pow5);
  1972. recovered_cache += middle_low.high();
  1973. auto high_to_middle = recovered_cache.high() << (64 - alpha);
  1974. auto middle_to_low = recovered_cache.low() << (64 - alpha);
  1975. recovered_cache = wuint::uint128{
  1976. (recovered_cache.low() >> alpha) | high_to_middle,
  1977. ((middle_low.low() >> alpha) | middle_to_low)
  1978. };
  1979. if (kb < 0) {
  1980. recovered_cache += 1;
  1981. }
  1982. // Get error
  1983. auto error_idx = (k - cache_holder<format>::min_k) / 16;
  1984. auto error = (compressed_cache_detail::errors[error_idx] >>
  1985. ((k - cache_holder<format>::min_k) % 16) * 2) & 0x3;
  1986. // Add the error back
  1987. assert(recovered_cache.low() + error >= recovered_cache.low());
  1988. recovered_cache = {
  1989. recovered_cache.high(),
  1990. recovered_cache.low() + error
  1991. };
  1992. return recovered_cache;
  1993. }
  1994. }
  1995. else
  1996. {
  1997. return cache_holder<format>::cache[std::size_t(k - cache_holder<format>::min_k)];
  1998. }
  1999. }
  2000. };
  2001. }
  2002. namespace input_validation {
  2003. struct base {};
  2004. struct assert_finite : base {
  2005. using input_validation_policy = assert_finite;
  2006. template <class Float>
  2007. static void validate_input([[maybe_unused]] ieee754_bits<Float> br) noexcept
  2008. {
  2009. assert(br.is_finite());
  2010. }
  2011. };
  2012. struct do_nothing : base {
  2013. using input_validation_policy = do_nothing;
  2014. template <class Float>
  2015. static void validate_input(ieee754_bits<Float>) noexcept {}
  2016. };
  2017. }
  2018. }
  2019. }
  2020. namespace policy {
  2021. namespace sign {
  2022. static constexpr auto ignore = detail::policy_impl::sign::ignore{};
  2023. static constexpr auto return_sign = detail::policy_impl::sign::return_sign{};
  2024. }
  2025. namespace trailing_zero {
  2026. static constexpr auto ignore = detail::policy_impl::trailing_zero::ignore{};
  2027. static constexpr auto remove = detail::policy_impl::trailing_zero::remove{};
  2028. static constexpr auto report = detail::policy_impl::trailing_zero::report{};
  2029. }
  2030. namespace rounding_mode {
  2031. static constexpr auto nearest_to_even =
  2032. detail::policy_impl::rounding_mode::nearest_to_even{};
  2033. static constexpr auto nearest_to_odd =
  2034. detail::policy_impl::rounding_mode::nearest_to_odd{};
  2035. static constexpr auto nearest_toward_plus_infinity =
  2036. detail::policy_impl::rounding_mode::nearest_toward_plus_infinity{};
  2037. static constexpr auto nearest_toward_minus_infinity =
  2038. detail::policy_impl::rounding_mode::nearest_toward_minus_infinity{};
  2039. static constexpr auto nearest_toward_zero =
  2040. detail::policy_impl::rounding_mode::nearest_toward_zero{};
  2041. static constexpr auto nearest_away_from_zero =
  2042. detail::policy_impl::rounding_mode::nearest_away_from_zero{};
  2043. static constexpr auto nearest_to_even_static_boundary =
  2044. detail::policy_impl::rounding_mode::nearest_to_even_static_boundary{};
  2045. static constexpr auto nearest_to_odd_static_boundary =
  2046. detail::policy_impl::rounding_mode::nearest_to_odd_static_boundary{};
  2047. static constexpr auto nearest_toward_plus_infinity_static_boundary =
  2048. detail::policy_impl::rounding_mode::nearest_toward_plus_infinity_static_boundary{};
  2049. static constexpr auto nearest_toward_minus_infinity_static_boundary =
  2050. detail::policy_impl::rounding_mode::nearest_toward_minus_infinity_static_boundary{};
  2051. static constexpr auto toward_plus_infinity =
  2052. detail::policy_impl::rounding_mode::toward_plus_infinity{};
  2053. static constexpr auto toward_minus_infinity =
  2054. detail::policy_impl::rounding_mode::toward_minus_infinity{};
  2055. static constexpr auto toward_zero =
  2056. detail::policy_impl::rounding_mode::toward_zero{};
  2057. static constexpr auto away_from_zero =
  2058. detail::policy_impl::rounding_mode::away_from_zero{};
  2059. }
  2060. namespace correct_rounding {
  2061. static constexpr auto do_not_care = detail::policy_impl::correct_rounding::do_not_care{};
  2062. static constexpr auto to_even = detail::policy_impl::correct_rounding::to_even{};
  2063. static constexpr auto to_odd = detail::policy_impl::correct_rounding::to_odd{};
  2064. static constexpr auto away_from_zero = detail::policy_impl::correct_rounding::away_from_zero{};
  2065. static constexpr auto toward_zero = detail::policy_impl::correct_rounding::toward_zero{};
  2066. }
  2067. namespace cache {
  2068. static constexpr auto normal = detail::policy_impl::cache::normal{};
  2069. static constexpr auto compressed = detail::policy_impl::cache::compressed{};
  2070. }
  2071. namespace input_validation {
  2072. static constexpr auto assert_finite = detail::policy_impl::input_validation::assert_finite{};
  2073. static constexpr auto do_nothing = detail::policy_impl::input_validation::do_nothing{};
  2074. }
  2075. }
  2076. namespace detail {
  2077. ////////////////////////////////////////////////////////////////////////////////////////
  2078. // The main algorithm
  2079. ////////////////////////////////////////////////////////////////////////////////////////
  2080. // Get sign/decimal significand/decimal exponent from
  2081. // the bit representation of a floating-point number
  2082. template <class Float>
  2083. struct impl : private ieee754_traits<Float>,
  2084. private ieee754_format_info<ieee754_traits<Float>::format>
  2085. {
  2086. using carrier_uint = typename ieee754_traits<Float>::carrier_uint;
  2087. using ieee754_traits<Float>::format;
  2088. using ieee754_traits<Float>::carrier_bits;
  2089. using ieee754_format_info<format>::significand_bits;
  2090. using ieee754_format_info<format>::min_exponent;
  2091. using ieee754_format_info<format>::max_exponent;
  2092. using ieee754_format_info<format>::exponent_bias;
  2093. using ieee754_format_info<format>::decimal_digits;
  2094. static constexpr int kappa = format == ieee754_format::binary32 ? 1 : 2;
  2095. static_assert(kappa >= 1);
  2096. static_assert(carrier_bits >= significand_bits + 2 + log::floor_log2_pow10(kappa + 1));
  2097. static constexpr int min_k = [] {
  2098. constexpr auto a = -log::floor_log10_pow2_minus_log10_4_over_3(
  2099. int(max_exponent - significand_bits));
  2100. constexpr auto b = -log::floor_log10_pow2(
  2101. int(max_exponent - significand_bits)) + kappa;
  2102. return a < b ? a : b;
  2103. }();
  2104. static_assert(min_k >= cache_holder<format>::min_k);
  2105. static constexpr int max_k = [] {
  2106. constexpr auto a = -log::floor_log10_pow2_minus_log10_4_over_3(
  2107. int(min_exponent - significand_bits + 1));
  2108. constexpr auto b = -log::floor_log10_pow2(
  2109. int(min_exponent - significand_bits)) + kappa;
  2110. return a > b ? a : b;
  2111. }();
  2112. static_assert(max_k <= cache_holder<format>::max_k);
  2113. using cache_entry_type =
  2114. typename cache_holder<format>::cache_entry_type;
  2115. static constexpr auto cache_bits =
  2116. cache_holder<format>::cache_bits;
  2117. static constexpr int max_power_of_factor_of_5 = log::floor_log5_pow2(int(significand_bits + 2));
  2118. static constexpr int divisibility_check_by_5_threshold =
  2119. log::floor_log2_pow10(max_power_of_factor_of_5 + kappa + 1);
  2120. static constexpr int case_fc_pm_half_lower_threshold = -kappa - log::floor_log5_pow2(kappa);
  2121. static constexpr int case_fc_pm_half_upper_threshold = log::floor_log2_pow10(kappa + 1);
  2122. static constexpr int case_fc_lower_threshold = -kappa - 1 - log::floor_log5_pow2(kappa + 1);
  2123. static constexpr int case_fc_upper_threshold = log::floor_log2_pow10(kappa + 1);
  2124. static constexpr int case_shorter_interval_left_endpoint_lower_threshold = 2;
  2125. static constexpr int case_shorter_interval_left_endpoint_upper_threshold = 2 +
  2126. log::floor_log2(compute_power<
  2127. count_factors<5>((carrier_uint(1) << (significand_bits + 2)) - 1) + 1
  2128. >(10) / 3);
  2129. static constexpr int case_shorter_interval_right_endpoint_lower_threshold = 0;
  2130. static constexpr int case_shorter_interval_right_endpoint_upper_threshold = 2 +
  2131. log::floor_log2(compute_power<
  2132. count_factors<5>((carrier_uint(1) << (significand_bits + 1)) + 1) + 1
  2133. >(10) / 3);
  2134. static constexpr int shorter_interval_tie_lower_threshold =
  2135. -log::floor_log5_pow2_minus_log5_3(significand_bits + 4) - 2 - significand_bits;
  2136. static constexpr int shorter_interval_tie_upper_threshold =
  2137. -log::floor_log5_pow2(significand_bits + 2) - 2 - significand_bits;
  2138. //// The main algorithm assumes the input is a normal/subnormal finite number
  2139. template <class ReturnType, class IntervalTypeProvider, class SignPolicy,
  2140. class TrailingZeroPolicy, class CorrectRoundingPolicy, class CachePolicy>
  2141. JKJ_SAFEBUFFERS static ReturnType compute_nearest(ieee754_bits<Float> const br) noexcept
  2142. {
  2143. //////////////////////////////////////////////////////////////////////
  2144. // Step 1: integer promotion & Schubfach multiplier calculation
  2145. //////////////////////////////////////////////////////////////////////
  2146. ReturnType ret_value;
  2147. SignPolicy::handle_sign(br, ret_value);
  2148. auto significand = br.extract_significand_bits();
  2149. auto exponent = int(br.extract_exponent_bits());
  2150. // Deal with normal/subnormal dichotomy
  2151. if (exponent != 0) {
  2152. exponent += exponent_bias - significand_bits;
  2153. // Shorter interval case; proceed like Schubfach
  2154. if (significand == 0) {
  2155. shorter_interval_case<TrailingZeroPolicy, CorrectRoundingPolicy, CachePolicy>(
  2156. ret_value, exponent,
  2157. IntervalTypeProvider::interval_type_shorter(br));
  2158. return ret_value;
  2159. }
  2160. significand |= (carrier_uint(1) << significand_bits);
  2161. }
  2162. // Subnormal case; interval is always regular
  2163. else {
  2164. exponent = min_exponent - significand_bits;
  2165. }
  2166. auto const interval_type = IntervalTypeProvider::interval_type_normal(br);
  2167. // Compute k and beta
  2168. int const minus_k = log::floor_log10_pow2(exponent) - kappa;
  2169. auto const cache = CachePolicy::template get_cache<format>(-minus_k);
  2170. int const beta_minus_1 = exponent + log::floor_log2_pow10(-minus_k);
  2171. // Compute zi and deltai
  2172. // 10^kappa <= deltai < 10^(kappa + 1)
  2173. auto const deltai = compute_delta(cache, beta_minus_1);
  2174. carrier_uint const two_fc = significand << 1;
  2175. carrier_uint const two_fr = two_fc | 1;
  2176. carrier_uint const zi = compute_mul(two_fr << beta_minus_1, cache);
  2177. //////////////////////////////////////////////////////////////////////
  2178. // Step 2: Try larger divisor; remove trailing zeros if necessary
  2179. //////////////////////////////////////////////////////////////////////
  2180. constexpr auto big_divisor = compute_power<kappa + 1>(std::uint32_t(10));
  2181. constexpr auto small_divisor = compute_power<kappa>(std::uint32_t(10));
  2182. // Using an upper bound on zi, we might be able to optimize the division
  2183. // better than the compiler; we are computing zi / big_divisor here
  2184. ret_value.significand = div::divide_by_pow10<kappa + 1,
  2185. significand_bits + kappa + 2, kappa + 1>(zi);
  2186. auto r = std::uint32_t(zi - big_divisor * ret_value.significand);
  2187. if (r > deltai) {
  2188. goto small_divisor_case_label;
  2189. }
  2190. else if (r < deltai) {
  2191. // Exclude the right endpoint if necessary
  2192. if (r == 0 && !interval_type.include_right_endpoint() &&
  2193. is_product_integer<integer_check_case_id::fc_pm_half>(two_fr, exponent, minus_k))
  2194. {
  2195. if constexpr (CorrectRoundingPolicy::tag ==
  2196. policy_impl::correct_rounding::tag_t::do_not_care)
  2197. {
  2198. ret_value.significand *= 10;
  2199. ret_value.exponent = minus_k + kappa;
  2200. --ret_value.significand;
  2201. return ret_value;
  2202. }
  2203. else {
  2204. --ret_value.significand;
  2205. r = big_divisor;
  2206. goto small_divisor_case_label;
  2207. }
  2208. }
  2209. }
  2210. else {
  2211. // r == deltai; compare fractional parts
  2212. // Check conditions in the order different from the paper
  2213. // to take advantage of short-circuiting
  2214. auto const two_fl = two_fc - 1;
  2215. if ((!interval_type.include_left_endpoint() ||
  2216. !is_product_integer<integer_check_case_id::fc_pm_half>(
  2217. two_fl, exponent, minus_k)) &&
  2218. !compute_mul_parity(two_fl, cache, beta_minus_1))
  2219. {
  2220. goto small_divisor_case_label;
  2221. }
  2222. }
  2223. ret_value.exponent = minus_k + kappa + 1;
  2224. // We may need to remove trailing zeros
  2225. TrailingZeroPolicy::on_trailing_zeros(ret_value);
  2226. return ret_value;
  2227. //////////////////////////////////////////////////////////////////////
  2228. // Step 3: Find the significand with the smaller divisor
  2229. //////////////////////////////////////////////////////////////////////
  2230. small_divisor_case_label:
  2231. TrailingZeroPolicy::no_trailing_zeros(ret_value);
  2232. ret_value.significand *= 10;
  2233. ret_value.exponent = minus_k + kappa;
  2234. constexpr auto mask = (std::uint32_t(1) << kappa) - 1;
  2235. if constexpr (CorrectRoundingPolicy::tag ==
  2236. policy_impl::correct_rounding::tag_t::do_not_care)
  2237. {
  2238. // Normally, we want to compute
  2239. // ret_value.significand += r / small_divisor
  2240. // and return, but we need to take care of the case that the resulting
  2241. // value is exactly the right endpoint, while that is not included in the interval
  2242. if (!interval_type.include_right_endpoint()) {
  2243. // Is r divisible by 2^kappa?
  2244. if ((r & mask) == 0) {
  2245. r >>= kappa;
  2246. // Is r divisible by 5^kappa?
  2247. if (div::check_divisibility_and_divide_by_pow5<kappa>(r) &&
  2248. is_product_integer<integer_check_case_id::fc_pm_half>(two_fr, exponent, minus_k))
  2249. {
  2250. // This should be in the interval
  2251. ret_value.significand += r - 1;
  2252. }
  2253. else {
  2254. ret_value.significand += r;
  2255. }
  2256. }
  2257. else {
  2258. ret_value.significand += div::small_division_by_pow10<kappa>(r);
  2259. }
  2260. }
  2261. else {
  2262. ret_value.significand += div::small_division_by_pow10<kappa>(r);
  2263. }
  2264. }
  2265. else
  2266. {
  2267. auto dist = r - (deltai / 2) + (small_divisor / 2);
  2268. // Is dist divisible by 2^kappa?
  2269. if ((dist & mask) == 0) {
  2270. bool const approx_y_parity = ((dist ^ (small_divisor / 2)) & 1) != 0;
  2271. dist >>= kappa;
  2272. // Is dist divisible by 5^kappa?
  2273. if (div::check_divisibility_and_divide_by_pow5<kappa>(dist)) {
  2274. ret_value.significand += dist;
  2275. // Check z^(f) >= epsilon^(f)
  2276. // We have either yi == zi - epsiloni or yi == (zi - epsiloni) - 1,
  2277. // where yi == zi - epsiloni if and only if z^(f) >= epsilon^(f)
  2278. // Since there are only 2 possibilities, we only need to care about the parity
  2279. // Also, zi and r should have the same parity since the divisor
  2280. // is an even number
  2281. if (compute_mul_parity(two_fc, cache, beta_minus_1) != approx_y_parity) {
  2282. --ret_value.significand;
  2283. }
  2284. else {
  2285. // If z^(f) >= epsilon^(f), we might have a tie
  2286. // when z^(f) == epsilon^(f), or equivalently, when y is an integer
  2287. // For tie-to-up case, we can just choose the upper one
  2288. if constexpr (CorrectRoundingPolicy::tag !=
  2289. policy_impl::correct_rounding::tag_t::away_from_zero)
  2290. {
  2291. if (is_product_integer<integer_check_case_id::fc>(
  2292. two_fc, exponent, minus_k))
  2293. {
  2294. CorrectRoundingPolicy::break_rounding_tie(ret_value);
  2295. }
  2296. }
  2297. }
  2298. }
  2299. // Is dist not divisible by 5^kappa?
  2300. else {
  2301. ret_value.significand += dist;
  2302. }
  2303. }
  2304. // Is dist not divisible by 2^kappa?
  2305. else {
  2306. // Since we know dist is small, we might be able to optimize the division
  2307. // better than the compiler; we are computing dist / small_divisor here
  2308. ret_value.significand += div::small_division_by_pow10<kappa>(dist);
  2309. }
  2310. }
  2311. return ret_value;
  2312. }
  2313. template <class TrailingZeroPolicy, class CorrectRoundingPolicy,
  2314. class CachePolicy, class ReturnType, class IntervalType>
  2315. JKJ_FORCEINLINE JKJ_SAFEBUFFERS static void shorter_interval_case(
  2316. ReturnType& ret_value, int const exponent, IntervalType const interval_type) noexcept
  2317. {
  2318. // Compute k and beta
  2319. int const minus_k = log::floor_log10_pow2_minus_log10_4_over_3(exponent);
  2320. int const beta_minus_1 = exponent + log::floor_log2_pow10(-minus_k);
  2321. // Compute xi and zi
  2322. auto const cache = CachePolicy::template get_cache<format>(-minus_k);
  2323. auto xi = compute_left_endpoint_for_shorter_interval_case(cache, beta_minus_1);
  2324. auto zi = compute_right_endpoint_for_shorter_interval_case(cache, beta_minus_1);
  2325. // If we don't accept the right endpoint and
  2326. // if the right endpoint is an integer, decrease it
  2327. if (!interval_type.include_right_endpoint() &&
  2328. is_right_endpoint_integer_shorter_interval(exponent))
  2329. {
  2330. --zi;
  2331. }
  2332. // If we don't accept the left endpoint or
  2333. // if the left endpoint is not an integer, increase it
  2334. if (!interval_type.include_left_endpoint() ||
  2335. !is_left_endpoint_integer_shorter_interval(exponent))
  2336. {
  2337. ++xi;
  2338. }
  2339. // Try bigger divisor
  2340. ret_value.significand = zi / 10;
  2341. // If succeed, remove trailing zeros if necessary and return
  2342. if (ret_value.significand * 10 >= xi) {
  2343. ret_value.exponent = minus_k + 1;
  2344. TrailingZeroPolicy::on_trailing_zeros(ret_value);
  2345. return;
  2346. }
  2347. // Otherwise, compute the round-up of y
  2348. TrailingZeroPolicy::no_trailing_zeros(ret_value);
  2349. ret_value.significand = compute_round_up_for_shorter_interval_case(cache, beta_minus_1);
  2350. ret_value.exponent = minus_k;
  2351. // When tie occurs, choose one of them according to the rule
  2352. if constexpr (CorrectRoundingPolicy::tag !=
  2353. policy_impl::correct_rounding::tag_t::do_not_care &&
  2354. CorrectRoundingPolicy::tag !=
  2355. policy_impl::correct_rounding::tag_t::away_from_zero)
  2356. {
  2357. if (exponent >= shorter_interval_tie_lower_threshold &&
  2358. exponent <= shorter_interval_tie_upper_threshold)
  2359. {
  2360. CorrectRoundingPolicy::break_rounding_tie(ret_value);
  2361. }
  2362. else if (ret_value.significand < xi) {
  2363. ++ret_value.significand;
  2364. }
  2365. }
  2366. else
  2367. {
  2368. if (ret_value.significand < xi) {
  2369. ++ret_value.significand;
  2370. }
  2371. }
  2372. }
  2373. template <class ReturnType, class SignPolicy, class TrailingZeroPolicy, class CachePolicy>
  2374. JKJ_SAFEBUFFERS static ReturnType
  2375. compute_left_closed_directed(ieee754_bits<Float> const br) noexcept
  2376. {
  2377. //////////////////////////////////////////////////////////////////////
  2378. // Step 1: integer promotion & Schubfach multiplier calculation
  2379. //////////////////////////////////////////////////////////////////////
  2380. ReturnType ret_value;
  2381. SignPolicy::handle_sign(br, ret_value);
  2382. auto significand = br.extract_significand_bits();
  2383. auto exponent = int(br.extract_exponent_bits());
  2384. // Deal with normal/subnormal dichotomy
  2385. if (exponent != 0) {
  2386. exponent += exponent_bias - significand_bits;
  2387. significand |= (carrier_uint(1) << significand_bits);
  2388. }
  2389. // Subnormal case; interval is always regular
  2390. else {
  2391. exponent = min_exponent - significand_bits;
  2392. }
  2393. // Compute k and beta
  2394. int const minus_k = log::floor_log10_pow2(exponent) - kappa;
  2395. auto const cache = CachePolicy::template get_cache<format>(-minus_k);
  2396. int const beta = exponent + log::floor_log2_pow10(-minus_k) + 1;
  2397. // Compute xi and deltai
  2398. // 10^kappa <= deltai < 10^(kappa + 1)
  2399. auto const deltai = compute_delta(cache, beta - 1);
  2400. carrier_uint xi = compute_mul(significand << beta, cache);
  2401. if (!is_product_integer<integer_check_case_id::fc>(significand, exponent + 1, minus_k)) {
  2402. ++xi;
  2403. }
  2404. //////////////////////////////////////////////////////////////////////
  2405. // Step 2: Try larger divisor; remove trailing zeros if necessary
  2406. //////////////////////////////////////////////////////////////////////
  2407. constexpr auto big_divisor = compute_power<kappa + 1>(std::uint32_t(10));
  2408. constexpr auto small_divisor = compute_power<kappa>(std::uint32_t(10));
  2409. // Using an upper bound on xi, we might be able to optimize the division
  2410. // better than the compiler; we are computing xi / big_divisor here
  2411. ret_value.significand = div::divide_by_pow10<kappa + 1,
  2412. significand_bits + kappa + 2, kappa + 1>(xi);
  2413. auto r = std::uint32_t(xi - big_divisor * ret_value.significand);
  2414. if (r != 0) {
  2415. ++ret_value.significand;
  2416. r = big_divisor - r;
  2417. }
  2418. if (r > deltai) {
  2419. goto small_divisor_case_label;
  2420. }
  2421. else if (r == deltai) {
  2422. // Compare the fractional parts
  2423. if (compute_mul_parity(significand + 1, cache, beta) ||
  2424. is_product_integer<integer_check_case_id::fc>(significand + 1, exponent + 1, minus_k))
  2425. {
  2426. goto small_divisor_case_label;
  2427. }
  2428. }
  2429. // The ceiling is inside, so we are done
  2430. ret_value.exponent = minus_k + kappa + 1;
  2431. TrailingZeroPolicy::on_trailing_zeros(ret_value);
  2432. return ret_value;
  2433. //////////////////////////////////////////////////////////////////////
  2434. // Step 3: Find the significand with the smaller divisor
  2435. //////////////////////////////////////////////////////////////////////
  2436. small_divisor_case_label:
  2437. ret_value.significand *= 10;
  2438. ret_value.significand -= div::small_division_by_pow10<kappa>(r);
  2439. ret_value.exponent = minus_k + kappa;
  2440. TrailingZeroPolicy::no_trailing_zeros(ret_value);
  2441. return ret_value;
  2442. }
  2443. template <class ReturnType, class SignPolicy, class TrailingZeroPolicy, class CachePolicy>
  2444. JKJ_SAFEBUFFERS static ReturnType
  2445. compute_right_closed_directed(ieee754_bits<Float> const br) noexcept
  2446. {
  2447. //////////////////////////////////////////////////////////////////////
  2448. // Step 1: integer promotion & Schubfach multiplier calculation
  2449. //////////////////////////////////////////////////////////////////////
  2450. ReturnType ret_value;
  2451. SignPolicy::handle_sign(br, ret_value);
  2452. auto significand = br.extract_significand_bits();
  2453. auto exponent = int(br.extract_exponent_bits());
  2454. // Deal with normal/subnormal dichotomy
  2455. bool closer_boundary = false;
  2456. if (exponent != 0) {
  2457. exponent += exponent_bias - significand_bits;
  2458. if (significand == 0) {
  2459. closer_boundary = true;
  2460. }
  2461. significand |= (carrier_uint(1) << significand_bits);
  2462. }
  2463. // Subnormal case; interval is always regular
  2464. else {
  2465. exponent = min_exponent - significand_bits;
  2466. }
  2467. // Compute k and beta
  2468. int const minus_k = log::floor_log10_pow2(exponent - (closer_boundary ? 1 : 0)) - kappa;
  2469. auto const cache = CachePolicy::template get_cache<format>(-minus_k);
  2470. int const beta = exponent + log::floor_log2_pow10(-minus_k) + 1;
  2471. // Compute zi and deltai
  2472. // 10^kappa <= deltai < 10^(kappa + 1)
  2473. auto const deltai = closer_boundary ?
  2474. compute_delta(cache, beta - 2) :
  2475. compute_delta(cache, beta - 1);
  2476. carrier_uint const zi = compute_mul(significand << beta, cache);
  2477. //////////////////////////////////////////////////////////////////////
  2478. // Step 2: Try larger divisor; remove trailing zeros if necessary
  2479. //////////////////////////////////////////////////////////////////////
  2480. constexpr auto big_divisor = compute_power<kappa + 1>(std::uint32_t(10));
  2481. constexpr auto small_divisor = compute_power<kappa>(std::uint32_t(10));
  2482. // Using an upper bound on zi, we might be able to optimize the division
  2483. // better than the compiler; we are computing zi / big_divisor here
  2484. ret_value.significand = div::divide_by_pow10<kappa + 1,
  2485. significand_bits + kappa + 2, kappa + 1>(zi);
  2486. auto const r = std::uint32_t(zi - big_divisor * ret_value.significand);
  2487. if (r > deltai) {
  2488. goto small_divisor_case_label;
  2489. }
  2490. else if (r == deltai) {
  2491. // Compare the fractional parts
  2492. if (closer_boundary) {
  2493. if (!compute_mul_parity((significand * 2) - 1, cache, beta - 1))
  2494. {
  2495. goto small_divisor_case_label;
  2496. }
  2497. }
  2498. else {
  2499. if (!compute_mul_parity(significand - 1, cache, beta))
  2500. {
  2501. goto small_divisor_case_label;
  2502. }
  2503. }
  2504. }
  2505. // The floor is inside, so we are done
  2506. ret_value.exponent = minus_k + kappa + 1;
  2507. TrailingZeroPolicy::on_trailing_zeros(ret_value);
  2508. return ret_value;
  2509. //////////////////////////////////////////////////////////////////////
  2510. // Step 3: Find the significand with the small divisor
  2511. //////////////////////////////////////////////////////////////////////
  2512. small_divisor_case_label:
  2513. ret_value.significand *= 10;
  2514. ret_value.significand += div::small_division_by_pow10<kappa>(r);
  2515. ret_value.exponent = minus_k + kappa;
  2516. TrailingZeroPolicy::no_trailing_zeros(ret_value);
  2517. return ret_value;
  2518. }
  2519. // Remove trailing zeros from n and return the number of zeros removed
  2520. JKJ_FORCEINLINE static int remove_trailing_zeros(carrier_uint& n) noexcept {
  2521. constexpr auto max_power = [] {
  2522. auto max_possible_significand =
  2523. std::numeric_limits<carrier_uint>::max() /
  2524. compute_power<kappa + 1>(std::uint32_t(10));
  2525. int k = 0;
  2526. carrier_uint p = 1;
  2527. while (p < max_possible_significand / 10) {
  2528. p *= 10;
  2529. ++k;
  2530. }
  2531. return k;
  2532. }();
  2533. auto t = bits::countr_zero(n);
  2534. if (t > max_power) {
  2535. t = max_power;
  2536. }
  2537. if constexpr (format == ieee754_format::binary32) {
  2538. constexpr auto const& divtable =
  2539. div::table_holder<carrier_uint, 5, decimal_digits>::table;
  2540. int s = 0;
  2541. for (; s < t - 1; s += 2) {
  2542. if (n * divtable.mod_inv[2] > divtable.max_quotients[2]) {
  2543. break;
  2544. }
  2545. n *= divtable.mod_inv[2];
  2546. }
  2547. if (s < t && n * divtable.mod_inv[1] <= divtable.max_quotients[1])
  2548. {
  2549. n *= divtable.mod_inv[1];
  2550. ++s;
  2551. }
  2552. n >>= s;
  2553. return s;
  2554. }
  2555. else {
  2556. static_assert(format == ieee754_format::binary64);
  2557. static_assert(kappa >= 2);
  2558. // Divide by 10^8 and reduce to 32-bits
  2559. // Since ret_value.significand <= (2^64 - 1) / 1000 < 10^17,
  2560. // both of the quotient and the r should fit in 32-bits
  2561. constexpr auto const& divtable =
  2562. div::table_holder<carrier_uint, 5, decimal_digits>::table;
  2563. // If the number is divisible by 1'0000'0000, work with the quotient
  2564. if (t >= 8) {
  2565. auto quotient_candidate = n * divtable.mod_inv[8];
  2566. if (quotient_candidate <= divtable.max_quotients[8]) {
  2567. auto quotient = std::uint32_t(quotient_candidate >> 8);
  2568. constexpr auto mod_inverse = std::uint32_t(divtable.mod_inv[1]);
  2569. constexpr auto max_quotient =
  2570. std::numeric_limits<std::uint32_t>::max() / 5;
  2571. int s = 8;
  2572. for (; s < t; ++s) {
  2573. if (quotient * mod_inverse > max_quotient) {
  2574. break;
  2575. }
  2576. quotient *= mod_inverse;
  2577. }
  2578. quotient >>= (s - 8);
  2579. n = quotient;
  2580. return s;
  2581. }
  2582. }
  2583. // Otherwise, work with the remainder
  2584. auto quotient = std::uint32_t(div::divide_by_pow10<8, 54, 0>(n));
  2585. auto remainder = std::uint32_t(n - 1'0000'0000 * quotient);
  2586. constexpr auto mod_inverse = std::uint32_t(divtable.mod_inv[1]);
  2587. constexpr auto max_quotient =
  2588. std::numeric_limits<std::uint32_t>::max() / 5;
  2589. if (t == 0 || remainder * mod_inverse > max_quotient) {
  2590. return 0;
  2591. }
  2592. remainder *= mod_inverse;
  2593. if (t == 1 || remainder * mod_inverse > max_quotient) {
  2594. n = (remainder >> 1)
  2595. + quotient * carrier_uint(1000'0000);
  2596. return 1;
  2597. }
  2598. remainder *= mod_inverse;
  2599. if (t == 2 || remainder * mod_inverse > max_quotient) {
  2600. n = (remainder >> 2)
  2601. + quotient * carrier_uint(100'0000);
  2602. return 2;
  2603. }
  2604. remainder *= mod_inverse;
  2605. if (t == 3 || remainder * mod_inverse > max_quotient) {
  2606. n = (remainder >> 3)
  2607. + quotient * carrier_uint(10'0000);
  2608. return 3;
  2609. }
  2610. remainder *= mod_inverse;
  2611. if (t == 4 || remainder * mod_inverse > max_quotient) {
  2612. n = (remainder >> 4)
  2613. + quotient * carrier_uint(1'0000);
  2614. return 4;
  2615. }
  2616. remainder *= mod_inverse;
  2617. if (t == 5 || remainder * mod_inverse > max_quotient) {
  2618. n = (remainder >> 5)
  2619. + quotient * carrier_uint(1000);
  2620. return 5;
  2621. }
  2622. remainder *= mod_inverse;
  2623. if (t == 6 || remainder * mod_inverse > max_quotient) {
  2624. n = (remainder >> 6)
  2625. + quotient * carrier_uint(100);
  2626. return 6;
  2627. }
  2628. remainder *= mod_inverse;
  2629. n = (remainder >> 7)
  2630. + quotient * carrier_uint(10);
  2631. return 7;
  2632. }
  2633. }
  2634. static carrier_uint compute_mul(carrier_uint u, cache_entry_type const& cache) noexcept
  2635. {
  2636. if constexpr (format == ieee754_format::binary32) {
  2637. return wuint::umul96_upper32(u, cache);
  2638. }
  2639. else {
  2640. return wuint::umul192_upper64(u, cache);
  2641. }
  2642. }
  2643. static std::uint32_t compute_delta(cache_entry_type const& cache, int beta_minus_1) noexcept
  2644. {
  2645. if constexpr (format == ieee754_format::binary32) {
  2646. return std::uint32_t(cache >> (cache_bits - 1 - beta_minus_1));
  2647. }
  2648. else {
  2649. return std::uint32_t(cache.high() >> (carrier_bits - 1 - beta_minus_1));
  2650. }
  2651. }
  2652. static bool compute_mul_parity(carrier_uint two_f, cache_entry_type const& cache, int beta_minus_1) noexcept
  2653. {
  2654. assert(beta_minus_1 >= 1);
  2655. assert(beta_minus_1 < 64);
  2656. if constexpr (format == ieee754_format::binary32) {
  2657. return ((wuint::umul96_lower64(two_f, cache) >>
  2658. (64 - beta_minus_1)) & 1) != 0;
  2659. }
  2660. else {
  2661. return ((wuint::umul192_middle64(two_f, cache) >>
  2662. (64 - beta_minus_1)) & 1) != 0;
  2663. }
  2664. }
  2665. static carrier_uint compute_left_endpoint_for_shorter_interval_case(
  2666. cache_entry_type const& cache, int beta_minus_1) noexcept
  2667. {
  2668. if constexpr (format == ieee754_format::binary32) {
  2669. return carrier_uint(
  2670. (cache - (cache >> (significand_bits + 2))) >>
  2671. (cache_bits - significand_bits - 1 - beta_minus_1));
  2672. }
  2673. else {
  2674. return (cache.high() - (cache.high() >> (significand_bits + 2))) >>
  2675. (carrier_bits - significand_bits - 1 - beta_minus_1);
  2676. }
  2677. }
  2678. static carrier_uint compute_right_endpoint_for_shorter_interval_case(
  2679. cache_entry_type const& cache, int beta_minus_1) noexcept
  2680. {
  2681. if constexpr (format == ieee754_format::binary32) {
  2682. return carrier_uint(
  2683. (cache + (cache >> (significand_bits + 1))) >>
  2684. (cache_bits - significand_bits - 1 - beta_minus_1));
  2685. }
  2686. else {
  2687. return (cache.high() + (cache.high() >> (significand_bits + 1))) >>
  2688. (carrier_bits - significand_bits - 1 - beta_minus_1);
  2689. }
  2690. }
  2691. static carrier_uint compute_round_up_for_shorter_interval_case(
  2692. cache_entry_type const& cache, int beta_minus_1) noexcept
  2693. {
  2694. if constexpr (format == ieee754_format::binary32) {
  2695. return (carrier_uint(cache >> (cache_bits - significand_bits - 2 - beta_minus_1)) + 1) / 2;
  2696. }
  2697. else {
  2698. return ((cache.high() >> (carrier_bits - significand_bits - 2 - beta_minus_1)) + 1) / 2;
  2699. }
  2700. }
  2701. static bool is_right_endpoint_integer_shorter_interval(int exponent) noexcept {
  2702. return exponent >= case_shorter_interval_right_endpoint_lower_threshold &&
  2703. exponent <= case_shorter_interval_right_endpoint_upper_threshold;
  2704. }
  2705. static bool is_left_endpoint_integer_shorter_interval(int exponent) noexcept {
  2706. return exponent >= case_shorter_interval_left_endpoint_lower_threshold &&
  2707. exponent <= case_shorter_interval_left_endpoint_upper_threshold;
  2708. }
  2709. enum class integer_check_case_id {
  2710. fc_pm_half,
  2711. fc
  2712. };
  2713. template <integer_check_case_id case_id>
  2714. static bool is_product_integer(carrier_uint two_f, int exponent, int minus_k) noexcept
  2715. {
  2716. // Case I: f = fc +- 1/2
  2717. if constexpr (case_id == integer_check_case_id::fc_pm_half)
  2718. {
  2719. if (exponent < case_fc_pm_half_lower_threshold) {
  2720. return false;
  2721. }
  2722. // For k >= 0
  2723. else if (exponent <= case_fc_pm_half_upper_threshold) {
  2724. return true;
  2725. }
  2726. // For k < 0
  2727. else if (exponent > divisibility_check_by_5_threshold) {
  2728. return false;
  2729. }
  2730. else {
  2731. return div::divisible_by_power_of_5<max_power_of_factor_of_5 + 1>(two_f, minus_k);
  2732. }
  2733. }
  2734. // Case II: f = fc + 1
  2735. // Case III: f = fc
  2736. else
  2737. {
  2738. // Exponent for 5 is negative
  2739. if (exponent > divisibility_check_by_5_threshold) {
  2740. return false;
  2741. }
  2742. else if (exponent > case_fc_upper_threshold) {
  2743. return div::divisible_by_power_of_5<max_power_of_factor_of_5 + 1>(two_f, minus_k);
  2744. }
  2745. // Both exponents are nonnegative
  2746. else if (exponent >= case_fc_lower_threshold) {
  2747. return true;
  2748. }
  2749. // Exponent for 2 is negative
  2750. else {
  2751. return div::divisible_by_power_of_2(two_f, minus_k - exponent + 1);
  2752. }
  2753. }
  2754. }
  2755. };
  2756. ////////////////////////////////////////////////////////////////////////////////////////
  2757. // Policy holder
  2758. ////////////////////////////////////////////////////////////////////////////////////////
  2759. namespace policy_impl {
  2760. // The library will specify a list of accepted kinds of policies and their defaults,
  2761. // and the user will pass a list of policies. The aim of helper classes/functions here
  2762. // is to do the following:
  2763. // 1. Check if the policy parameters given by the user are all valid; that means,
  2764. // each of them should be of the kinds specified by the library.
  2765. // If that's not the case, then the compilation fails.
  2766. // 2. Check if multiple policy parameters for the same kind is specified by the user.
  2767. // If that's the case, then the compilation fails.
  2768. // 3. Build a class deriving from all policies the user have given, and also from
  2769. // the default policies if the user did not specify one for some kinds.
  2770. // A policy belongs to a certain kind if it is deriving from a base class.
  2771. // For a given kind, find a policy belonging to that kind.
  2772. // Check if there are more than one such policies.
  2773. enum class policy_found_info {
  2774. not_found, unique, repeated
  2775. };
  2776. template <class Policy, policy_found_info info>
  2777. struct found_policy_pair {
  2778. using policy = Policy;
  2779. static constexpr auto found_info = info;
  2780. };
  2781. template <class Base, class DefaultPolicy>
  2782. struct base_default_pair {
  2783. using base = Base;
  2784. template <class FoundPolicyInfo>
  2785. static constexpr FoundPolicyInfo get_policy_impl(FoundPolicyInfo) {
  2786. return{};
  2787. }
  2788. template <class FoundPolicyInfo, class FirstPolicy, class... RemainingPolicies>
  2789. static constexpr auto get_policy_impl(FoundPolicyInfo, FirstPolicy, RemainingPolicies... remainings) {
  2790. if constexpr (std::is_base_of_v<Base, FirstPolicy>) {
  2791. if constexpr (FoundPolicyInfo::found_info == policy_found_info::not_found) {
  2792. return get_policy_impl(
  2793. found_policy_pair<FirstPolicy, policy_found_info::unique>{},
  2794. remainings...);
  2795. }
  2796. else {
  2797. return get_policy_impl(
  2798. found_policy_pair<FirstPolicy, policy_found_info::repeated>{},
  2799. remainings...);
  2800. }
  2801. }
  2802. else {
  2803. return get_policy_impl(FoundPolicyInfo{},
  2804. remainings...);
  2805. }
  2806. }
  2807. template <class... Policies>
  2808. static constexpr auto get_policy(Policies... policies) {
  2809. return get_policy_impl(
  2810. found_policy_pair<DefaultPolicy, policy_found_info::not_found>{},
  2811. policies...);
  2812. }
  2813. };
  2814. template <class... BaseDefaultPairs>
  2815. struct base_default_pair_list {};
  2816. // Check if a given policy belongs to one of the kinds specified by the library
  2817. template <class Policy>
  2818. constexpr bool check_policy_validity(Policy, base_default_pair_list<>)
  2819. {
  2820. return false;
  2821. }
  2822. template <class Policy, class FirstBaseDefaultPair, class... RemainingBaseDefaultPairs>
  2823. constexpr bool check_policy_validity(Policy,
  2824. base_default_pair_list<FirstBaseDefaultPair, RemainingBaseDefaultPairs...>)
  2825. {
  2826. return std::is_base_of_v<typename FirstBaseDefaultPair::base, Policy> ||
  2827. check_policy_validity(Policy{}, base_default_pair_list< RemainingBaseDefaultPairs...>{});
  2828. }
  2829. template <class BaseDefaultPairList>
  2830. constexpr bool check_policy_list_validity(BaseDefaultPairList) {
  2831. return true;
  2832. }
  2833. template <class BaseDefaultPairList, class FirstPolicy, class... RemainingPolicies>
  2834. constexpr bool check_policy_list_validity(BaseDefaultPairList,
  2835. FirstPolicy, RemainingPolicies... remaining_policies)
  2836. {
  2837. return check_policy_validity(FirstPolicy{}, BaseDefaultPairList{}) &&
  2838. check_policy_list_validity(BaseDefaultPairList{}, remaining_policies...);
  2839. }
  2840. // Build policy_holder
  2841. template <bool repeated_, class... FoundPolicyPairs>
  2842. struct found_policy_pair_list {
  2843. static constexpr bool repeated = repeated_;
  2844. };
  2845. template <class... Policies>
  2846. struct policy_holder : Policies... {};
  2847. template <bool repeated, class... FoundPolicyPairs, class... Policies>
  2848. constexpr auto make_policy_holder_impl(
  2849. base_default_pair_list<>,
  2850. found_policy_pair_list<repeated, FoundPolicyPairs...>,
  2851. Policies...)
  2852. {
  2853. return found_policy_pair_list<repeated, FoundPolicyPairs...>{};
  2854. }
  2855. template <class FirstBaseDefaultPair, class... RemainingBaseDefaultPairs,
  2856. bool repeated, class... FoundPolicyPairs, class... Policies>
  2857. constexpr auto make_policy_holder_impl(
  2858. base_default_pair_list<FirstBaseDefaultPair, RemainingBaseDefaultPairs...>,
  2859. found_policy_pair_list<repeated, FoundPolicyPairs...>,
  2860. Policies... policies)
  2861. {
  2862. using new_found_policy_pair = decltype(FirstBaseDefaultPair::get_policy(policies...));
  2863. return make_policy_holder_impl(
  2864. base_default_pair_list<RemainingBaseDefaultPairs...>{},
  2865. found_policy_pair_list<
  2866. repeated || new_found_policy_pair::found_info == policy_found_info::repeated,
  2867. new_found_policy_pair, FoundPolicyPairs...
  2868. >{}, policies...);
  2869. }
  2870. template <bool repeated, class... RawPolicies>
  2871. constexpr auto convert_to_policy_holder(found_policy_pair_list<repeated>, RawPolicies...) {
  2872. return policy_holder<RawPolicies...>{};
  2873. }
  2874. template <bool repeated, class FirstFoundPolicyPair, class... RemainingFoundPolicyPairs, class... RawPolicies>
  2875. constexpr auto convert_to_policy_holder(
  2876. found_policy_pair_list<repeated, FirstFoundPolicyPair, RemainingFoundPolicyPairs...>, RawPolicies... policies)
  2877. {
  2878. return convert_to_policy_holder(found_policy_pair_list<repeated, RemainingFoundPolicyPairs...>{},
  2879. typename FirstFoundPolicyPair::policy{}, policies...);
  2880. }
  2881. template <class BaseDefaultPairList, class... Policies>
  2882. constexpr auto make_policy_holder(BaseDefaultPairList, Policies... policies) {
  2883. static_assert(check_policy_list_validity(BaseDefaultPairList{}, Policies{}...),
  2884. "jkj::dragonbox: an invalid policy is specified");
  2885. using policy_pair_list = decltype(make_policy_holder_impl(BaseDefaultPairList{},
  2886. found_policy_pair_list<false>{}, policies...));
  2887. static_assert(!policy_pair_list::repeated,
  2888. "jkj::dragonbox: each policy should be specified at most once");
  2889. return convert_to_policy_holder(policy_pair_list{});
  2890. }
  2891. }
  2892. }
  2893. ////////////////////////////////////////////////////////////////////////////////////////
  2894. // The interface function
  2895. ////////////////////////////////////////////////////////////////////////////////////////
  2896. template <class Float, class... Policies>
  2897. JKJ_SAFEBUFFERS JKJ_FORCEINLINE auto to_decimal(Float x, Policies... policies)
  2898. {
  2899. // Build policy holder type
  2900. using namespace detail::policy_impl;
  2901. using policy_holder = decltype(make_policy_holder(
  2902. base_default_pair_list<
  2903. base_default_pair<sign::base, sign::return_sign>,
  2904. base_default_pair<trailing_zero::base, trailing_zero::remove>,
  2905. base_default_pair<rounding_mode::base, rounding_mode::nearest_to_even>,
  2906. base_default_pair<correct_rounding::base, correct_rounding::to_even>,
  2907. base_default_pair<cache::base, cache::normal>,
  2908. base_default_pair<input_validation::base, input_validation::assert_finite>
  2909. >{}, policies...));
  2910. using return_type = fp_t<Float,
  2911. policy_holder::return_has_sign,
  2912. policy_holder::report_trailing_zeros>;
  2913. auto br = ieee754_bits(x);
  2914. policy_holder::validate_input(br);
  2915. return policy_holder::delegate(br,
  2916. [br](auto interval_type_provider) {
  2917. constexpr auto tag = decltype(interval_type_provider)::tag;
  2918. if constexpr (tag == rounding_mode::tag_t::to_nearest) {
  2919. return detail::impl<Float>::template
  2920. compute_nearest<return_type, decltype(interval_type_provider),
  2921. typename policy_holder::sign_policy,
  2922. typename policy_holder::trailing_zero_policy,
  2923. typename policy_holder::correct_rounding_policy,
  2924. typename policy_holder::cache_policy
  2925. >(br);
  2926. }
  2927. else if constexpr (tag == rounding_mode::tag_t::left_closed_directed) {
  2928. return detail::impl<Float>::template
  2929. compute_left_closed_directed<return_type,
  2930. typename policy_holder::sign_policy,
  2931. typename policy_holder::trailing_zero_policy,
  2932. typename policy_holder::cache_policy
  2933. >(br);
  2934. }
  2935. else {
  2936. return detail::impl<Float>::template
  2937. compute_right_closed_directed<return_type,
  2938. typename policy_holder::sign_policy,
  2939. typename policy_holder::trailing_zero_policy,
  2940. typename policy_holder::cache_policy
  2941. >(br);
  2942. }
  2943. });
  2944. }
  2945. }
  2946. #undef JKJ_HAS_COUNTR_ZERO_INTRINSIC
  2947. #undef JKJ_FORCEINLINE
  2948. #undef JKJ_SAFEBUFFERS
  2949. #endif