cpp_dec_float.hpp 132 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // Copyright Christopher Kormanyos 2002 - 2013.
  3. // Copyright 2011 -2013 John Maddock. Distributed under the Boost
  4. // Software License, Version 1.0. (See accompanying file
  5. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. //
  7. // This work is based on an earlier work:
  8. // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
  9. // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
  10. //
  11. // Note that there are no "noexcept" specifications on the functions in this file: there are too many
  12. // calls to lexical_cast (and similar) to easily analyse the code for correctness. So until compilers
  13. // can detect noexcept misuse at compile time, the only realistic option is to simply not use it here.
  14. //
  15. #ifndef BOOST_MP_CPP_DEC_FLOAT_BACKEND_HPP
  16. #define BOOST_MP_CPP_DEC_FLOAT_BACKEND_HPP
  17. #include <boost/config.hpp>
  18. #include <boost/cstdint.hpp>
  19. #include <limits>
  20. #ifndef BOOST_NO_CXX11_HDR_ARRAY
  21. #include <array>
  22. #else
  23. #include <boost/array.hpp>
  24. #endif
  25. #include <boost/cstdint.hpp>
  26. #include <boost/functional/hash_fwd.hpp>
  27. #include <boost/multiprecision/number.hpp>
  28. #include <boost/multiprecision/detail/big_lanczos.hpp>
  29. #include <boost/multiprecision/detail/dynamic_array.hpp>
  30. //
  31. // Headers required for Boost.Math integration:
  32. //
  33. #include <boost/math/policies/policy.hpp>
  34. //
  35. // Some includes we need from Boost.Math, since we rely on that library to provide these functions:
  36. //
  37. #include <boost/math/special_functions/asinh.hpp>
  38. #include <boost/math/special_functions/acosh.hpp>
  39. #include <boost/math/special_functions/atanh.hpp>
  40. #include <boost/math/special_functions/cbrt.hpp>
  41. #include <boost/math/special_functions/expm1.hpp>
  42. #include <boost/math/special_functions/gamma.hpp>
  43. #ifdef BOOST_MSVC
  44. #pragma warning(push)
  45. #pragma warning(disable:6326) // comparison of two constants
  46. #endif
  47. namespace boost{
  48. namespace multiprecision{
  49. namespace backends{
  50. template <unsigned Digits10, class ExponentType = boost::int32_t, class Allocator = void>
  51. class cpp_dec_float;
  52. } // namespace
  53. template <unsigned Digits10, class ExponentType, class Allocator>
  54. struct number_category<backends::cpp_dec_float<Digits10, ExponentType, Allocator> > : public mpl::int_<number_kind_floating_point>{};
  55. namespace backends{
  56. template <unsigned Digits10, class ExponentType, class Allocator>
  57. class cpp_dec_float
  58. {
  59. private:
  60. static const boost::int32_t cpp_dec_float_digits10_setting = Digits10;
  61. // We need at least 16-bits in the exponent type to do anything sensible:
  62. BOOST_STATIC_ASSERT_MSG(boost::is_signed<ExponentType>::value, "ExponentType must be a signed built in integer type.");
  63. BOOST_STATIC_ASSERT_MSG(sizeof(ExponentType) > 1, "ExponentType is too small.");
  64. public:
  65. typedef mpl::list<boost::long_long_type> signed_types;
  66. typedef mpl::list<boost::ulong_long_type> unsigned_types;
  67. typedef mpl::list<long double> float_types;
  68. typedef ExponentType exponent_type;
  69. static const boost::int32_t cpp_dec_float_radix = 10L;
  70. static const boost::int32_t cpp_dec_float_digits10_limit_lo = 9L;
  71. static const boost::int32_t cpp_dec_float_digits10_limit_hi = boost::integer_traits<boost::int32_t>::const_max - 100;
  72. static const boost::int32_t cpp_dec_float_digits10 = ((cpp_dec_float_digits10_setting < cpp_dec_float_digits10_limit_lo) ? cpp_dec_float_digits10_limit_lo : ((cpp_dec_float_digits10_setting > cpp_dec_float_digits10_limit_hi) ? cpp_dec_float_digits10_limit_hi : cpp_dec_float_digits10_setting));
  73. static const ExponentType cpp_dec_float_max_exp10 = (static_cast<ExponentType>(1) << (std::numeric_limits<ExponentType>::digits - 5));
  74. static const ExponentType cpp_dec_float_min_exp10 = -cpp_dec_float_max_exp10;
  75. static const ExponentType cpp_dec_float_max_exp = cpp_dec_float_max_exp10;
  76. static const ExponentType cpp_dec_float_min_exp = cpp_dec_float_min_exp10;
  77. BOOST_STATIC_ASSERT((cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp10 == -cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp10));
  78. private:
  79. static const boost::int32_t cpp_dec_float_elem_digits10 = 8L;
  80. static const boost::int32_t cpp_dec_float_elem_mask = 100000000L;
  81. BOOST_STATIC_ASSERT(0 == cpp_dec_float_max_exp10 % cpp_dec_float_elem_digits10);
  82. // There are three guard limbs.
  83. // 1) The first limb has 'play' from 1...8 decimal digits.
  84. // 2) The last limb also has 'play' from 1...8 decimal digits.
  85. // 3) One limb can get lost when justifying after multiply,
  86. // as only half of the triangle is multiplied and a carry
  87. // from below is missing.
  88. static const boost::int32_t cpp_dec_float_elem_number_request = static_cast<boost::int32_t>((cpp_dec_float_digits10 / cpp_dec_float_elem_digits10) + (((cpp_dec_float_digits10 % cpp_dec_float_elem_digits10) != 0) ? 1 : 0));
  89. // The number of elements needed (with a minimum of two) plus three added guard limbs.
  90. static const boost::int32_t cpp_dec_float_elem_number = static_cast<boost::int32_t>(((cpp_dec_float_elem_number_request < 2L) ? 2L : cpp_dec_float_elem_number_request) + 3L);
  91. public:
  92. static const boost::int32_t cpp_dec_float_total_digits10 = static_cast<boost::int32_t>(cpp_dec_float_elem_number * cpp_dec_float_elem_digits10);
  93. private:
  94. typedef enum enum_fpclass_type
  95. {
  96. cpp_dec_float_finite,
  97. cpp_dec_float_inf,
  98. cpp_dec_float_NaN
  99. }
  100. fpclass_type;
  101. #ifndef BOOST_NO_CXX11_HDR_ARRAY
  102. typedef typename mpl::if_<is_void<Allocator>,
  103. std::array<boost::uint32_t, cpp_dec_float_elem_number>,
  104. detail::dynamic_array<boost::uint32_t, cpp_dec_float_elem_number, Allocator>
  105. >::type array_type;
  106. #else
  107. typedef typename mpl::if_<is_void<Allocator>,
  108. boost::array<boost::uint32_t, cpp_dec_float_elem_number>,
  109. detail::dynamic_array<boost::uint32_t, cpp_dec_float_elem_number, Allocator>
  110. >::type array_type;
  111. #endif
  112. array_type data;
  113. ExponentType exp;
  114. bool neg;
  115. fpclass_type fpclass;
  116. boost::int32_t prec_elem;
  117. //
  118. // Special values constructor:
  119. //
  120. cpp_dec_float(fpclass_type c) :
  121. data(),
  122. exp (static_cast<ExponentType>(0)),
  123. neg (false),
  124. fpclass (c),
  125. prec_elem(cpp_dec_float_elem_number) { }
  126. //
  127. // Static data initializer:
  128. //
  129. struct initializer
  130. {
  131. initializer()
  132. {
  133. cpp_dec_float<Digits10, ExponentType, Allocator>::nan();
  134. cpp_dec_float<Digits10, ExponentType, Allocator>::inf();
  135. (cpp_dec_float<Digits10, ExponentType, Allocator>::min)();
  136. (cpp_dec_float<Digits10, ExponentType, Allocator>::max)();
  137. cpp_dec_float<Digits10, ExponentType, Allocator>::zero();
  138. cpp_dec_float<Digits10, ExponentType, Allocator>::one();
  139. cpp_dec_float<Digits10, ExponentType, Allocator>::two();
  140. cpp_dec_float<Digits10, ExponentType, Allocator>::half();
  141. cpp_dec_float<Digits10, ExponentType, Allocator>::double_min();
  142. cpp_dec_float<Digits10, ExponentType, Allocator>::double_max();
  143. cpp_dec_float<Digits10, ExponentType, Allocator>::long_double_max();
  144. cpp_dec_float<Digits10, ExponentType, Allocator>::long_double_min();
  145. cpp_dec_float<Digits10, ExponentType, Allocator>::long_long_max();
  146. cpp_dec_float<Digits10, ExponentType, Allocator>::long_long_min();
  147. cpp_dec_float<Digits10, ExponentType, Allocator>::ulong_long_max();
  148. cpp_dec_float<Digits10, ExponentType, Allocator>::eps();
  149. cpp_dec_float<Digits10, ExponentType, Allocator>::pow2(0);
  150. }
  151. void do_nothing(){}
  152. };
  153. static initializer init;
  154. public:
  155. // Constructors
  156. cpp_dec_float() BOOST_MP_NOEXCEPT_IF(noexcept(array_type())) :
  157. data(),
  158. exp (static_cast<ExponentType>(0)),
  159. neg (false),
  160. fpclass (cpp_dec_float_finite),
  161. prec_elem(cpp_dec_float_elem_number) { }
  162. cpp_dec_float(const char* s) :
  163. data(),
  164. exp (static_cast<ExponentType>(0)),
  165. neg (false),
  166. fpclass (cpp_dec_float_finite),
  167. prec_elem(cpp_dec_float_elem_number)
  168. {
  169. *this = s;
  170. }
  171. template<class I>
  172. cpp_dec_float(I i, typename enable_if<is_unsigned<I> >::type* = 0) :
  173. data(),
  174. exp (static_cast<ExponentType>(0)),
  175. neg (false),
  176. fpclass (cpp_dec_float_finite),
  177. prec_elem(cpp_dec_float_elem_number)
  178. {
  179. from_unsigned_long_long(i);
  180. }
  181. template <class I>
  182. cpp_dec_float(I i, typename enable_if<is_signed<I> >::type* = 0) :
  183. data(),
  184. exp (static_cast<ExponentType>(0)),
  185. neg (false),
  186. fpclass (cpp_dec_float_finite),
  187. prec_elem(cpp_dec_float_elem_number)
  188. {
  189. if(i < 0)
  190. {
  191. from_unsigned_long_long(boost::multiprecision::detail::unsigned_abs(i));
  192. negate();
  193. }
  194. else
  195. from_unsigned_long_long(i);
  196. }
  197. cpp_dec_float(const cpp_dec_float& f) BOOST_MP_NOEXCEPT_IF(noexcept(array_type(std::declval<const array_type&>()))) :
  198. data (f.data),
  199. exp (f.exp),
  200. neg (f.neg),
  201. fpclass (f.fpclass),
  202. prec_elem(f.prec_elem) { }
  203. template <unsigned D, class ET, class A>
  204. cpp_dec_float(const cpp_dec_float<D, ET, A>& f, typename enable_if_c<D <= Digits10>::type* = 0) :
  205. data(),
  206. exp (f.exp),
  207. neg (f.neg),
  208. fpclass (static_cast<fpclass_type>(static_cast<int>(f.fpclass))),
  209. prec_elem(cpp_dec_float_elem_number)
  210. {
  211. std::copy(f.data.begin(), f.data.begin() + f.prec_elem, data.begin());
  212. }
  213. template <unsigned D, class ET, class A>
  214. explicit cpp_dec_float(const cpp_dec_float<D, ET, A>& f, typename disable_if_c<D <= Digits10>::type* = 0) :
  215. data(),
  216. exp (f.exp),
  217. neg (f.neg),
  218. fpclass (static_cast<fpclass_type>(static_cast<int>(f.fpclass))),
  219. prec_elem(cpp_dec_float_elem_number)
  220. {
  221. // TODO: this doesn't round!
  222. std::copy(f.data.begin(), f.data.begin() + prec_elem, data.begin());
  223. }
  224. template <class F>
  225. cpp_dec_float(const F val, typename enable_if_c<is_floating_point<F>::value
  226. #ifdef BOOST_HAS_FLOAT128
  227. && !boost::is_same<F, __float128>::value
  228. #endif
  229. >::type* = 0) :
  230. data(),
  231. exp (static_cast<ExponentType>(0)),
  232. neg (false),
  233. fpclass (cpp_dec_float_finite),
  234. prec_elem(cpp_dec_float_elem_number)
  235. {
  236. *this = val;
  237. }
  238. cpp_dec_float(const double mantissa, const ExponentType exponent);
  239. std::size_t hash()const
  240. {
  241. std::size_t result = 0;
  242. for(int i = 0; i < prec_elem; ++i)
  243. boost::hash_combine(result, data[i]);
  244. boost::hash_combine(result, exp);
  245. boost::hash_combine(result, neg);
  246. boost::hash_combine(result, fpclass);
  247. return result;
  248. }
  249. // Specific special values.
  250. static const cpp_dec_float& nan()
  251. {
  252. static const cpp_dec_float val(cpp_dec_float_NaN);
  253. init.do_nothing();
  254. return val;
  255. }
  256. static const cpp_dec_float& inf()
  257. {
  258. static const cpp_dec_float val(cpp_dec_float_inf);
  259. init.do_nothing();
  260. return val;
  261. }
  262. static const cpp_dec_float& (max)()
  263. {
  264. init.do_nothing();
  265. static cpp_dec_float val_max = std::string("1.0e" + boost::lexical_cast<std::string>(cpp_dec_float_max_exp10)).c_str();
  266. return val_max;
  267. }
  268. static const cpp_dec_float& (min)()
  269. {
  270. init.do_nothing();
  271. static cpp_dec_float val_min = std::string("1.0e" + boost::lexical_cast<std::string>(cpp_dec_float_min_exp10)).c_str();
  272. return val_min;
  273. }
  274. static const cpp_dec_float& zero()
  275. {
  276. init.do_nothing();
  277. static cpp_dec_float val(static_cast<boost::ulong_long_type>(0u));
  278. return val;
  279. }
  280. static const cpp_dec_float& one()
  281. {
  282. init.do_nothing();
  283. static cpp_dec_float val(static_cast<boost::ulong_long_type>(1u));
  284. return val;
  285. }
  286. static const cpp_dec_float& two()
  287. {
  288. init.do_nothing();
  289. static cpp_dec_float val(static_cast<boost::ulong_long_type>(2u));
  290. return val;
  291. }
  292. static const cpp_dec_float& half()
  293. {
  294. init.do_nothing();
  295. static cpp_dec_float val(0.5L);
  296. return val;
  297. }
  298. static const cpp_dec_float& double_min()
  299. {
  300. init.do_nothing();
  301. static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::min)()));
  302. return val;
  303. }
  304. static const cpp_dec_float& double_max()
  305. {
  306. init.do_nothing();
  307. static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::max)()));
  308. return val;
  309. }
  310. static const cpp_dec_float& long_double_min()
  311. {
  312. init.do_nothing();
  313. #ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  314. static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::min)()));
  315. #else
  316. static cpp_dec_float val((std::numeric_limits<long double>::min)());
  317. #endif
  318. return val;
  319. }
  320. static const cpp_dec_float& long_double_max()
  321. {
  322. init.do_nothing();
  323. #ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  324. static cpp_dec_float val(static_cast<long double>((std::numeric_limits<double>::max)()));
  325. #else
  326. static cpp_dec_float val((std::numeric_limits<long double>::max)());
  327. #endif
  328. return val;
  329. }
  330. static const cpp_dec_float& long_long_max()
  331. {
  332. init.do_nothing();
  333. static cpp_dec_float val((std::numeric_limits<boost::long_long_type>::max)());
  334. return val;
  335. }
  336. static const cpp_dec_float& long_long_min()
  337. {
  338. init.do_nothing();
  339. static cpp_dec_float val((std::numeric_limits<boost::long_long_type>::min)());
  340. return val;
  341. }
  342. static const cpp_dec_float& ulong_long_max()
  343. {
  344. init.do_nothing();
  345. static cpp_dec_float val((std::numeric_limits<boost::ulong_long_type>::max)());
  346. return val;
  347. }
  348. static const cpp_dec_float& eps()
  349. {
  350. init.do_nothing();
  351. static cpp_dec_float val(1.0, 1 - static_cast<int>(cpp_dec_float_digits10));
  352. return val;
  353. }
  354. // Basic operations.
  355. cpp_dec_float& operator=(const cpp_dec_float& v) BOOST_MP_NOEXCEPT_IF(noexcept(std::declval<array_type&>() = std::declval<const array_type&>()))
  356. {
  357. data = v.data;
  358. exp = v.exp;
  359. neg = v.neg;
  360. fpclass = v.fpclass;
  361. prec_elem = v.prec_elem;
  362. return *this;
  363. }
  364. template <unsigned D>
  365. cpp_dec_float& operator=(const cpp_dec_float<D>& f)
  366. {
  367. exp = f.exp;
  368. neg = f.neg;
  369. fpclass = static_cast<enum_fpclass_type>(static_cast<int>(f.fpclass));
  370. unsigned elems = (std::min)(f.prec_elem, cpp_dec_float_elem_number);
  371. std::copy(f.data.begin(), f.data.begin() + elems, data.begin());
  372. std::fill(data.begin() + elems, data.end(), 0);
  373. prec_elem = cpp_dec_float_elem_number;
  374. return *this;
  375. }
  376. cpp_dec_float& operator=(boost::long_long_type v)
  377. {
  378. if(v < 0)
  379. {
  380. from_unsigned_long_long(1u - boost::ulong_long_type(v + 1)); // Avoid undefined behaviour in negation of minimum value for long long
  381. negate();
  382. }
  383. else
  384. from_unsigned_long_long(v);
  385. return *this;
  386. }
  387. cpp_dec_float& operator=(boost::ulong_long_type v)
  388. {
  389. from_unsigned_long_long(v);
  390. return *this;
  391. }
  392. cpp_dec_float& operator=(long double v);
  393. cpp_dec_float& operator=(const char* v)
  394. {
  395. rd_string(v);
  396. return *this;
  397. }
  398. cpp_dec_float& operator+=(const cpp_dec_float& v);
  399. cpp_dec_float& operator-=(const cpp_dec_float& v);
  400. cpp_dec_float& operator*=(const cpp_dec_float& v);
  401. cpp_dec_float& operator/=(const cpp_dec_float& v);
  402. cpp_dec_float& add_unsigned_long_long(const boost::ulong_long_type n)
  403. {
  404. cpp_dec_float t;
  405. t.from_unsigned_long_long(n);
  406. return *this += t;
  407. }
  408. cpp_dec_float& sub_unsigned_long_long(const boost::ulong_long_type n)
  409. {
  410. cpp_dec_float t;
  411. t.from_unsigned_long_long(n);
  412. return *this -= t;
  413. }
  414. cpp_dec_float& mul_unsigned_long_long(const boost::ulong_long_type n);
  415. cpp_dec_float& div_unsigned_long_long(const boost::ulong_long_type n);
  416. // Elementary primitives.
  417. cpp_dec_float& calculate_inv ();
  418. cpp_dec_float& calculate_sqrt();
  419. void negate()
  420. {
  421. if(!iszero())
  422. neg = !neg;
  423. }
  424. // Comparison functions
  425. bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION() const { return (fpclass == cpp_dec_float_NaN); }
  426. bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION() const { return (fpclass == cpp_dec_float_inf); }
  427. bool isfinite BOOST_PREVENT_MACRO_SUBSTITUTION() const { return (fpclass == cpp_dec_float_finite); }
  428. bool iszero () const
  429. {
  430. return ((fpclass == cpp_dec_float_finite) && (data[0u] == 0u));
  431. }
  432. bool isone () const;
  433. bool isint () const;
  434. bool isneg () const { return neg; }
  435. // Operators pre-increment and pre-decrement
  436. cpp_dec_float& operator++()
  437. {
  438. return *this += one();
  439. }
  440. cpp_dec_float& operator--()
  441. {
  442. return *this -= one();
  443. }
  444. std::string str(boost::intmax_t digits, std::ios_base::fmtflags f)const;
  445. int compare(const cpp_dec_float& v)const;
  446. template <class V>
  447. int compare(const V& v)const
  448. {
  449. cpp_dec_float<Digits10, ExponentType, Allocator> t;
  450. t = v;
  451. return compare(t);
  452. }
  453. void swap(cpp_dec_float& v)
  454. {
  455. data.swap(v.data);
  456. std::swap(exp, v.exp);
  457. std::swap(neg, v.neg);
  458. std::swap(fpclass, v.fpclass);
  459. std::swap(prec_elem, v.prec_elem);
  460. }
  461. double extract_double() const;
  462. long double extract_long_double() const;
  463. boost::long_long_type extract_signed_long_long() const;
  464. boost::ulong_long_type extract_unsigned_long_long() const;
  465. void extract_parts(double& mantissa, ExponentType& exponent) const;
  466. cpp_dec_float extract_integer_part() const;
  467. void precision(const boost::int32_t prec_digits)
  468. {
  469. if(prec_digits >= cpp_dec_float_total_digits10)
  470. {
  471. prec_elem = cpp_dec_float_elem_number;
  472. }
  473. else
  474. {
  475. const boost::int32_t elems = static_cast<boost::int32_t>( static_cast<boost::int32_t>( (prec_digits + (cpp_dec_float_elem_digits10 / 2)) / cpp_dec_float_elem_digits10)
  476. + static_cast<boost::int32_t>(((prec_digits % cpp_dec_float_elem_digits10) != 0) ? 1 : 0));
  477. prec_elem = (std::min)(cpp_dec_float_elem_number, (std::max)(elems, static_cast<boost::int32_t>(2)));
  478. }
  479. }
  480. static cpp_dec_float pow2(boost::long_long_type i);
  481. ExponentType order()const
  482. {
  483. const bool bo_order_is_zero = ((!(isfinite)()) || (data[0] == static_cast<boost::uint32_t>(0u)));
  484. //
  485. // Binary search to find the order of the leading term:
  486. //
  487. ExponentType prefix = 0;
  488. if(data[0] >= 100000UL)
  489. {
  490. if(data[0] >= 10000000UL)
  491. {
  492. if(data[0] >= 100000000UL)
  493. {
  494. if(data[0] >= 1000000000UL)
  495. prefix = 9;
  496. else
  497. prefix = 8;
  498. }
  499. else
  500. prefix = 7;
  501. }
  502. else
  503. {
  504. if(data[0] >= 1000000UL)
  505. prefix = 6;
  506. else
  507. prefix = 5;
  508. }
  509. }
  510. else
  511. {
  512. if(data[0] >= 1000UL)
  513. {
  514. if(data[0] >= 10000UL)
  515. prefix = 4;
  516. else
  517. prefix = 3;
  518. }
  519. else
  520. {
  521. if(data[0] >= 100)
  522. prefix = 2;
  523. else if(data[0] >= 10)
  524. prefix = 1;
  525. }
  526. }
  527. return (bo_order_is_zero ? static_cast<ExponentType>(0) : static_cast<ExponentType>(exp + prefix));
  528. }
  529. template<class Archive>
  530. void serialize(Archive & ar, const unsigned int /*version*/)
  531. {
  532. for(unsigned i = 0; i < data.size(); ++i)
  533. ar & boost::serialization::make_nvp("digit", data[i]);
  534. ar & boost::serialization::make_nvp("exponent", exp);
  535. ar & boost::serialization::make_nvp("sign", neg);
  536. ar & boost::serialization::make_nvp("class-type", fpclass);
  537. ar & boost::serialization::make_nvp("precision", prec_elem);
  538. }
  539. private:
  540. static bool data_elem_is_non_zero_predicate(const boost::uint32_t& d) { return (d != static_cast<boost::uint32_t>(0u)); }
  541. static bool data_elem_is_non_nine_predicate(const boost::uint32_t& d) { return (d != static_cast<boost::uint32_t>(cpp_dec_float::cpp_dec_float_elem_mask - 1)); }
  542. static bool char_is_nonzero_predicate(const char& c) { return (c != static_cast<char>('0')); }
  543. void from_unsigned_long_long(const boost::ulong_long_type u);
  544. int cmp_data(const array_type& vd) const;
  545. static boost::uint32_t mul_loop_uv(boost::uint32_t* const u, const boost::uint32_t* const v, const boost::int32_t p);
  546. static boost::uint32_t mul_loop_n (boost::uint32_t* const u, boost::uint32_t n, const boost::int32_t p);
  547. static boost::uint32_t div_loop_n (boost::uint32_t* const u, boost::uint32_t n, const boost::int32_t p);
  548. bool rd_string(const char* const s);
  549. template <unsigned D, class ET, class A>
  550. friend class cpp_dec_float;
  551. };
  552. template <unsigned Digits10, class ExponentType, class Allocator>
  553. typename cpp_dec_float<Digits10, ExponentType, Allocator>::initializer cpp_dec_float<Digits10, ExponentType, Allocator>::init;
  554. template <unsigned Digits10, class ExponentType, class Allocator>
  555. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_radix;
  556. template <unsigned Digits10, class ExponentType, class Allocator>
  557. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10_setting;
  558. template <unsigned Digits10, class ExponentType, class Allocator>
  559. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10_limit_lo;
  560. template <unsigned Digits10, class ExponentType, class Allocator>
  561. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10_limit_hi;
  562. template <unsigned Digits10, class ExponentType, class Allocator>
  563. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10;
  564. template <unsigned Digits10, class ExponentType, class Allocator>
  565. const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp;
  566. template <unsigned Digits10, class ExponentType, class Allocator>
  567. const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp;
  568. template <unsigned Digits10, class ExponentType, class Allocator>
  569. const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp10;
  570. template <unsigned Digits10, class ExponentType, class Allocator>
  571. const ExponentType cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp10;
  572. template <unsigned Digits10, class ExponentType, class Allocator>
  573. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_digits10;
  574. template <unsigned Digits10, class ExponentType, class Allocator>
  575. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_number_request;
  576. template <unsigned Digits10, class ExponentType, class Allocator>
  577. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_number;
  578. template <unsigned Digits10, class ExponentType, class Allocator>
  579. const boost::int32_t cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_elem_mask;
  580. template <unsigned Digits10, class ExponentType, class Allocator>
  581. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator+=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)
  582. {
  583. if((isnan)())
  584. {
  585. return *this;
  586. }
  587. if((isinf)())
  588. {
  589. if((v.isinf)() && (isneg() != v.isneg()))
  590. {
  591. *this = nan();
  592. }
  593. return *this;
  594. }
  595. if(iszero())
  596. {
  597. return operator=(v);
  598. }
  599. if((v.isnan)() || (v.isinf)())
  600. {
  601. *this = v;
  602. return *this;
  603. }
  604. // Get the offset for the add/sub operation.
  605. static const ExponentType max_delta_exp = static_cast<ExponentType>((cpp_dec_float_elem_number - 1) * cpp_dec_float_elem_digits10);
  606. const ExponentType ofs_exp = static_cast<ExponentType>(exp - v.exp);
  607. // Check if the operation is out of range, requiring special handling.
  608. if(v.iszero() || (ofs_exp > max_delta_exp))
  609. {
  610. // Result is *this unchanged since v is negligible compared to *this.
  611. return *this;
  612. }
  613. else if(ofs_exp < -max_delta_exp)
  614. {
  615. // Result is *this = v since *this is negligible compared to v.
  616. return operator=(v);
  617. }
  618. // Do the add/sub operation.
  619. typename array_type::iterator p_u = data.begin();
  620. typename array_type::const_iterator p_v = v.data.begin();
  621. bool b_copy = false;
  622. const boost::int32_t ofs = static_cast<boost::int32_t>(static_cast<boost::int32_t>(ofs_exp) / cpp_dec_float_elem_digits10);
  623. array_type n_data;
  624. if(neg == v.neg)
  625. {
  626. // Add v to *this, where the data array of either *this or v
  627. // might have to be treated with a positive, negative or zero offset.
  628. // The result is stored in *this. The data are added one element
  629. // at a time, each element with carry.
  630. if(ofs >= static_cast<boost::int32_t>(0))
  631. {
  632. std::copy(v.data.begin(), v.data.end() - static_cast<size_t>(ofs), n_data.begin() + static_cast<size_t>(ofs));
  633. std::fill(n_data.begin(), n_data.begin() + static_cast<size_t>(ofs), static_cast<boost::uint32_t>(0u));
  634. p_v = n_data.begin();
  635. }
  636. else
  637. {
  638. std::copy(data.begin(), data.end() - static_cast<size_t>(-ofs), n_data.begin() + static_cast<size_t>(-ofs));
  639. std::fill(n_data.begin(), n_data.begin() + static_cast<size_t>(-ofs), static_cast<boost::uint32_t>(0u));
  640. p_u = n_data.begin();
  641. b_copy = true;
  642. }
  643. // Addition algorithm
  644. boost::uint32_t carry = static_cast<boost::uint32_t>(0u);
  645. for(boost::int32_t j = static_cast<boost::int32_t>(cpp_dec_float_elem_number - static_cast<boost::int32_t>(1)); j >= static_cast<boost::int32_t>(0); j--)
  646. {
  647. boost::uint32_t t = static_cast<boost::uint32_t>(static_cast<boost::uint32_t>(p_u[j] + p_v[j]) + carry);
  648. carry = t / static_cast<boost::uint32_t>(cpp_dec_float_elem_mask);
  649. p_u[j] = static_cast<boost::uint32_t>(t - static_cast<boost::uint32_t>(carry * static_cast<boost::uint32_t>(cpp_dec_float_elem_mask)));
  650. }
  651. if(b_copy)
  652. {
  653. data = n_data;
  654. exp = v.exp;
  655. }
  656. // There needs to be a carry into the element -1 of the array data
  657. if(carry != static_cast<boost::uint32_t>(0u))
  658. {
  659. std::copy_backward(data.begin(), data.end() - static_cast<std::size_t>(1u), data.end());
  660. data[0] = carry;
  661. exp += static_cast<ExponentType>(cpp_dec_float_elem_digits10);
  662. }
  663. }
  664. else
  665. {
  666. // Subtract v from *this, where the data array of either *this or v
  667. // might have to be treated with a positive, negative or zero offset.
  668. if((ofs > static_cast<boost::int32_t>(0))
  669. || ( (ofs == static_cast<boost::int32_t>(0))
  670. && (cmp_data(v.data) > static_cast<boost::int32_t>(0)))
  671. )
  672. {
  673. // In this case, |u| > |v| and ofs is positive.
  674. // Copy the data of v, shifted down to a lower value
  675. // into the data array m_n. Set the operand pointer p_v
  676. // to point to the copied, shifted data m_n.
  677. std::copy(v.data.begin(), v.data.end() - static_cast<size_t>(ofs), n_data.begin() + static_cast<size_t>(ofs));
  678. std::fill(n_data.begin(), n_data.begin() + static_cast<size_t>(ofs), static_cast<boost::uint32_t>(0u));
  679. p_v = n_data.begin();
  680. }
  681. else
  682. {
  683. if(ofs != static_cast<boost::int32_t>(0))
  684. {
  685. // In this case, |u| < |v| and ofs is negative.
  686. // Shift the data of u down to a lower value.
  687. std::copy_backward(data.begin(), data.end() - static_cast<size_t>(-ofs), data.end());
  688. std::fill(data.begin(), data.begin() + static_cast<size_t>(-ofs), static_cast<boost::uint32_t>(0u));
  689. }
  690. // Copy the data of v into the data array n_data.
  691. // Set the u-pointer p_u to point to m_n and the
  692. // operand pointer p_v to point to the shifted
  693. // data m_data.
  694. n_data = v.data;
  695. p_u = n_data.begin();
  696. p_v = data.begin();
  697. b_copy = true;
  698. }
  699. boost::int32_t j;
  700. // Subtraction algorithm
  701. boost::int32_t borrow = static_cast<boost::int32_t>(0);
  702. for(j = static_cast<boost::int32_t>(cpp_dec_float_elem_number - static_cast<boost::int32_t>(1)); j >= static_cast<boost::int32_t>(0); j--)
  703. {
  704. boost::int32_t t = static_cast<boost::int32_t>(static_cast<boost::int32_t>( static_cast<boost::int32_t>(p_u[j])
  705. - static_cast<boost::int32_t>(p_v[j])) - borrow);
  706. // Underflow? Borrow?
  707. if(t < static_cast<boost::int32_t>(0))
  708. {
  709. // Yes, underflow and borrow
  710. t += static_cast<boost::int32_t>(cpp_dec_float_elem_mask);
  711. borrow = static_cast<boost::int32_t>(1);
  712. }
  713. else
  714. {
  715. borrow = static_cast<boost::int32_t>(0);
  716. }
  717. p_u[j] = static_cast<boost::uint32_t>(static_cast<boost::uint32_t>(t) % static_cast<boost::uint32_t>(cpp_dec_float_elem_mask));
  718. }
  719. if(b_copy)
  720. {
  721. data = n_data;
  722. exp = v.exp;
  723. neg = v.neg;
  724. }
  725. // Is it necessary to justify the data?
  726. const typename array_type::const_iterator first_nonzero_elem = std::find_if(data.begin(), data.end(), data_elem_is_non_zero_predicate);
  727. if(first_nonzero_elem != data.begin())
  728. {
  729. if(first_nonzero_elem == data.end())
  730. {
  731. // This result of the subtraction is exactly zero.
  732. // Reset the sign and the exponent.
  733. neg = false;
  734. exp = static_cast<ExponentType>(0);
  735. }
  736. else
  737. {
  738. // Justify the data
  739. const std::size_t sj = static_cast<std::size_t>(std::distance<typename array_type::const_iterator>(data.begin(), first_nonzero_elem));
  740. std::copy(data.begin() + static_cast<std::size_t>(sj), data.end(), data.begin());
  741. std::fill(data.end() - sj, data.end(), static_cast<boost::uint32_t>(0u));
  742. exp -= static_cast<ExponentType>(sj * static_cast<std::size_t>(cpp_dec_float_elem_digits10));
  743. }
  744. }
  745. }
  746. // Handle underflow.
  747. if(iszero())
  748. return (*this = zero());
  749. // Check for potential overflow.
  750. const bool b_result_might_overflow = (exp >= static_cast<ExponentType>(cpp_dec_float_max_exp10));
  751. // Handle overflow.
  752. if(b_result_might_overflow)
  753. {
  754. const bool b_result_is_neg = neg;
  755. neg = false;
  756. if(compare((cpp_dec_float::max)()) > 0)
  757. *this = inf();
  758. neg = b_result_is_neg;
  759. }
  760. return *this;
  761. }
  762. template <unsigned Digits10, class ExponentType, class Allocator>
  763. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator-=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)
  764. {
  765. // Use *this - v = -(-*this + v).
  766. negate();
  767. *this += v;
  768. negate();
  769. return *this;
  770. }
  771. template <unsigned Digits10, class ExponentType, class Allocator>
  772. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator*=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)
  773. {
  774. // Evaluate the sign of the result.
  775. const bool b_result_is_neg = (neg != v.neg);
  776. // Artificially set the sign of the result to be positive.
  777. neg = false;
  778. // Handle special cases like zero, inf and NaN.
  779. const bool b_u_is_inf = (isinf)();
  780. const bool b_v_is_inf = (v.isinf)();
  781. const bool b_u_is_zero = iszero();
  782. const bool b_v_is_zero = v.iszero();
  783. if( ((isnan)() || (v.isnan)())
  784. || (b_u_is_inf && b_v_is_zero)
  785. || (b_v_is_inf && b_u_is_zero)
  786. )
  787. {
  788. *this = nan();
  789. return *this;
  790. }
  791. if(b_u_is_inf || b_v_is_inf)
  792. {
  793. *this = inf();
  794. if(b_result_is_neg)
  795. negate();
  796. return *this;
  797. }
  798. if(b_u_is_zero || b_v_is_zero)
  799. {
  800. return *this = zero();
  801. }
  802. // Check for potential overflow or underflow.
  803. const bool b_result_might_overflow = ((exp + v.exp) >= static_cast<ExponentType>(cpp_dec_float_max_exp10));
  804. const bool b_result_might_underflow = ((exp + v.exp) <= static_cast<ExponentType>(cpp_dec_float_min_exp10));
  805. // Set the exponent of the result.
  806. exp += v.exp;
  807. const boost::int32_t prec_mul = (std::min)(prec_elem, v.prec_elem);
  808. const boost::uint32_t carry = mul_loop_uv(data.data(), v.data.data(), prec_mul);
  809. // Handle a potential carry.
  810. if(carry != static_cast<boost::uint32_t>(0u))
  811. {
  812. exp += cpp_dec_float_elem_digits10;
  813. // Shift the result of the multiplication one element to the right...
  814. std::copy_backward(data.begin(),
  815. data.begin() + static_cast<std::size_t>(prec_elem - static_cast<boost::int32_t>(1)),
  816. data.begin() + static_cast<std::size_t>(prec_elem));
  817. // ... And insert the carry.
  818. data.front() = carry;
  819. }
  820. // Handle overflow.
  821. if(b_result_might_overflow && (compare((cpp_dec_float::max)()) > 0))
  822. {
  823. *this = inf();
  824. }
  825. // Handle underflow.
  826. if(b_result_might_underflow && (compare((cpp_dec_float::min)()) < 0))
  827. {
  828. *this = zero();
  829. return *this;
  830. }
  831. // Set the sign of the result.
  832. neg = b_result_is_neg;
  833. return *this;
  834. }
  835. template <unsigned Digits10, class ExponentType, class Allocator>
  836. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator/=(const cpp_dec_float<Digits10, ExponentType, Allocator>& v)
  837. {
  838. if(iszero())
  839. {
  840. if((v.isnan)())
  841. {
  842. return *this = v;
  843. }
  844. else if(v.iszero())
  845. {
  846. return *this = nan();
  847. }
  848. }
  849. const bool u_and_v_are_finite_and_identical = ( (isfinite)()
  850. && (fpclass == v.fpclass)
  851. && (exp == v.exp)
  852. && (cmp_data(v.data) == static_cast<boost::int32_t>(0)));
  853. if(u_and_v_are_finite_and_identical)
  854. {
  855. if(neg != v.neg)
  856. {
  857. *this = one();
  858. negate();
  859. }
  860. else
  861. *this = one();
  862. return *this;
  863. }
  864. else
  865. {
  866. cpp_dec_float t(v);
  867. t.calculate_inv();
  868. return operator*=(t);
  869. }
  870. }
  871. template <unsigned Digits10, class ExponentType, class Allocator>
  872. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::mul_unsigned_long_long(const boost::ulong_long_type n)
  873. {
  874. // Multiply *this with a constant boost::ulong_long_type.
  875. // Evaluate the sign of the result.
  876. const bool b_neg = neg;
  877. // Artificially set the sign of the result to be positive.
  878. neg = false;
  879. // Handle special cases like zero, inf and NaN.
  880. const bool b_u_is_inf = (isinf)();
  881. const bool b_n_is_zero = (n == static_cast<boost::int32_t>(0));
  882. if((isnan)() || (b_u_is_inf && b_n_is_zero))
  883. {
  884. return (*this = nan());
  885. }
  886. if(b_u_is_inf)
  887. {
  888. *this = inf();
  889. if(b_neg)
  890. negate();
  891. return *this;
  892. }
  893. if(iszero() || b_n_is_zero)
  894. {
  895. // Multiplication by zero.
  896. return *this = zero();
  897. }
  898. if(n >= static_cast<boost::ulong_long_type>(cpp_dec_float_elem_mask))
  899. {
  900. neg = b_neg;
  901. cpp_dec_float t;
  902. t = n;
  903. return operator*=(t);
  904. }
  905. if(n == static_cast<boost::ulong_long_type>(1u))
  906. {
  907. neg = b_neg;
  908. return *this;
  909. }
  910. // Set up the multiplication loop.
  911. const boost::uint32_t nn = static_cast<boost::uint32_t>(n);
  912. const boost::uint32_t carry = mul_loop_n(data.data(), nn, prec_elem);
  913. // Handle the carry and adjust the exponent.
  914. if(carry != static_cast<boost::uint32_t>(0u))
  915. {
  916. exp += static_cast<ExponentType>(cpp_dec_float_elem_digits10);
  917. // Shift the result of the multiplication one element to the right.
  918. std::copy_backward(data.begin(),
  919. data.begin() + static_cast<std::size_t>(prec_elem - static_cast<boost::int32_t>(1)),
  920. data.begin() + static_cast<std::size_t>(prec_elem));
  921. data.front() = static_cast<boost::uint32_t>(carry);
  922. }
  923. // Check for potential overflow.
  924. const bool b_result_might_overflow = (exp >= cpp_dec_float_max_exp10);
  925. // Handle overflow.
  926. if(b_result_might_overflow && (compare((cpp_dec_float::max)()) > 0))
  927. {
  928. *this = inf();
  929. }
  930. // Set the sign.
  931. neg = b_neg;
  932. return *this;
  933. }
  934. template <unsigned Digits10, class ExponentType, class Allocator>
  935. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::div_unsigned_long_long(const boost::ulong_long_type n)
  936. {
  937. // Divide *this by a constant boost::ulong_long_type.
  938. // Evaluate the sign of the result.
  939. const bool b_neg = neg;
  940. // Artificially set the sign of the result to be positive.
  941. neg = false;
  942. // Handle special cases like zero, inf and NaN.
  943. if((isnan)())
  944. {
  945. return *this;
  946. }
  947. if((isinf)())
  948. {
  949. *this = inf();
  950. if(b_neg)
  951. negate();
  952. return *this;
  953. }
  954. if(n == static_cast<boost::ulong_long_type>(0u))
  955. {
  956. // Divide by 0.
  957. if(iszero())
  958. {
  959. *this = nan();
  960. return *this;
  961. }
  962. else
  963. {
  964. *this = inf();
  965. if(isneg())
  966. negate();
  967. return *this;
  968. }
  969. }
  970. if(iszero())
  971. {
  972. return *this;
  973. }
  974. if(n >= static_cast<boost::ulong_long_type>(cpp_dec_float_elem_mask))
  975. {
  976. neg = b_neg;
  977. cpp_dec_float t;
  978. t = n;
  979. return operator/=(t);
  980. }
  981. const boost::uint32_t nn = static_cast<boost::uint32_t>(n);
  982. if(nn > static_cast<boost::uint32_t>(1u))
  983. {
  984. // Do the division loop.
  985. const boost::uint32_t prev = div_loop_n(data.data(), nn, prec_elem);
  986. // Determine if one leading zero is in the result data.
  987. if(data[0] == static_cast<boost::uint32_t>(0u))
  988. {
  989. // Adjust the exponent
  990. exp -= static_cast<ExponentType>(cpp_dec_float_elem_digits10);
  991. // Shift result of the division one element to the left.
  992. std::copy(data.begin() + static_cast<std::size_t>(1u),
  993. data.begin() + static_cast<std::size_t>(prec_elem - static_cast<boost::int32_t>(1)),
  994. data.begin());
  995. data[prec_elem - static_cast<boost::int32_t>(1)] = static_cast<boost::uint32_t>(static_cast<boost::uint64_t>(prev * static_cast<boost::uint64_t>(cpp_dec_float_elem_mask)) / nn);
  996. }
  997. }
  998. // Check for potential underflow.
  999. const bool b_result_might_underflow = (exp <= cpp_dec_float_min_exp10);
  1000. // Handle underflow.
  1001. if(b_result_might_underflow && (compare((cpp_dec_float::min)()) < 0))
  1002. return (*this = zero());
  1003. // Set the sign of the result.
  1004. neg = b_neg;
  1005. return *this;
  1006. }
  1007. template <unsigned Digits10, class ExponentType, class Allocator>
  1008. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::calculate_inv()
  1009. {
  1010. // Compute the inverse of *this.
  1011. const bool b_neg = neg;
  1012. neg = false;
  1013. // Handle special cases like zero, inf and NaN.
  1014. if(iszero())
  1015. {
  1016. *this = inf();
  1017. if(b_neg)
  1018. negate();
  1019. return *this;
  1020. }
  1021. if((isnan)())
  1022. {
  1023. return *this;
  1024. }
  1025. if((isinf)())
  1026. {
  1027. return *this = zero();
  1028. }
  1029. if(isone())
  1030. {
  1031. if(b_neg)
  1032. negate();
  1033. return *this;
  1034. }
  1035. // Save the original *this.
  1036. cpp_dec_float<Digits10, ExponentType, Allocator> x(*this);
  1037. // Generate the initial estimate using division.
  1038. // Extract the mantissa and exponent for a "manual"
  1039. // computation of the estimate.
  1040. double dd;
  1041. ExponentType ne;
  1042. x.extract_parts(dd, ne);
  1043. // Do the inverse estimate using double precision estimates of mantissa and exponent.
  1044. operator=(cpp_dec_float<Digits10, ExponentType, Allocator>(1.0 / dd, -ne));
  1045. // Compute the inverse of *this. Quadratically convergent Newton-Raphson iteration
  1046. // is used. During the iterative steps, the precision of the calculation is limited
  1047. // to the minimum required in order to minimize the run-time.
  1048. static const boost::int32_t double_digits10_minus_a_few = std::numeric_limits<double>::digits10 - 3;
  1049. for(boost::int32_t digits = double_digits10_minus_a_few; digits <= cpp_dec_float_total_digits10; digits *= static_cast<boost::int32_t>(2))
  1050. {
  1051. // Adjust precision of the terms.
  1052. precision(static_cast<boost::int32_t>((digits + 10) * static_cast<boost::int32_t>(2)));
  1053. x.precision(static_cast<boost::int32_t>((digits + 10) * static_cast<boost::int32_t>(2)));
  1054. // Next iteration.
  1055. cpp_dec_float t(*this);
  1056. t *= x;
  1057. t -= two();
  1058. t.negate();
  1059. *this *= t;
  1060. }
  1061. neg = b_neg;
  1062. prec_elem = cpp_dec_float_elem_number;
  1063. return *this;
  1064. }
  1065. template <unsigned Digits10, class ExponentType, class Allocator>
  1066. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::calculate_sqrt()
  1067. {
  1068. // Compute the square root of *this.
  1069. if((isinf)() && !isneg())
  1070. {
  1071. return *this;
  1072. }
  1073. if(isneg() || (!(isfinite)()))
  1074. {
  1075. *this = nan();
  1076. errno = EDOM;
  1077. return *this;
  1078. }
  1079. if(iszero() || isone())
  1080. {
  1081. return *this;
  1082. }
  1083. // Save the original *this.
  1084. cpp_dec_float<Digits10, ExponentType, Allocator> x(*this);
  1085. // Generate the initial estimate using division.
  1086. // Extract the mantissa and exponent for a "manual"
  1087. // computation of the estimate.
  1088. double dd;
  1089. ExponentType ne;
  1090. extract_parts(dd, ne);
  1091. // Force the exponent to be an even multiple of two.
  1092. if((ne % static_cast<ExponentType>(2)) != static_cast<ExponentType>(0))
  1093. {
  1094. ++ne;
  1095. dd /= 10.0;
  1096. }
  1097. // Setup the iteration.
  1098. // Estimate the square root using simple manipulations.
  1099. const double sqd = std::sqrt(dd);
  1100. *this = cpp_dec_float<Digits10, ExponentType, Allocator>(sqd, static_cast<ExponentType>(ne / static_cast<ExponentType>(2)));
  1101. // Estimate 1.0 / (2.0 * x0) using simple manipulations.
  1102. cpp_dec_float<Digits10, ExponentType, Allocator> vi(0.5 / sqd, static_cast<ExponentType>(-ne / static_cast<ExponentType>(2)));
  1103. // Compute the square root of x. Coupled Newton iteration
  1104. // as described in "Pi Unleashed" is used. During the
  1105. // iterative steps, the precision of the calculation is
  1106. // limited to the minimum required in order to minimize
  1107. // the run-time.
  1108. //
  1109. // Book references:
  1110. // http://www.jjj.de/pibook/pibook.html
  1111. // http://www.amazon.com/exec/obidos/tg/detail/-/3540665722/qid=1035535482/sr=8-7/ref=sr_8_7/104-3357872-6059916?v=glance&n=507846
  1112. static const boost::uint32_t double_digits10_minus_a_few = std::numeric_limits<double>::digits10 - 3;
  1113. for(boost::int32_t digits = double_digits10_minus_a_few; digits <= cpp_dec_float_total_digits10; digits *= 2u)
  1114. {
  1115. // Adjust precision of the terms.
  1116. precision((digits + 10) * 2);
  1117. vi.precision((digits + 10) * 2);
  1118. // Next iteration of vi
  1119. cpp_dec_float t(*this);
  1120. t *= vi;
  1121. t.negate();
  1122. t.mul_unsigned_long_long(2u);
  1123. t += one();
  1124. t *= vi;
  1125. vi += t;
  1126. // Next iteration of *this
  1127. t = *this;
  1128. t *= *this;
  1129. t.negate();
  1130. t += x;
  1131. t *= vi;
  1132. *this += t;
  1133. }
  1134. prec_elem = cpp_dec_float_elem_number;
  1135. return *this;
  1136. }
  1137. template <unsigned Digits10, class ExponentType, class Allocator>
  1138. int cpp_dec_float<Digits10, ExponentType, Allocator>::cmp_data(const array_type& vd) const
  1139. {
  1140. // Compare the data of *this with those of v.
  1141. // Return +1 for *this > v
  1142. // 0 for *this = v
  1143. // -1 for *this < v
  1144. const std::pair<typename array_type::const_iterator, typename array_type::const_iterator> mismatch_pair = std::mismatch(data.begin(), data.end(), vd.begin());
  1145. const bool is_equal = ((mismatch_pair.first == data.end()) && (mismatch_pair.second == vd.end()));
  1146. if(is_equal)
  1147. {
  1148. return 0;
  1149. }
  1150. else
  1151. {
  1152. return ((*mismatch_pair.first > *mismatch_pair.second) ? 1 : -1);
  1153. }
  1154. }
  1155. template <unsigned Digits10, class ExponentType, class Allocator>
  1156. int cpp_dec_float<Digits10, ExponentType, Allocator>::compare(const cpp_dec_float& v) const
  1157. {
  1158. // Compare v with *this.
  1159. // Return +1 for *this > v
  1160. // 0 for *this = v
  1161. // -1 for *this < v
  1162. // Handle all non-finite cases.
  1163. if((!(isfinite)()) || (!(v.isfinite)()))
  1164. {
  1165. // NaN can never equal NaN. Return an implementation-dependent
  1166. // signed result. Also note that comparison of NaN with NaN
  1167. // using operators greater-than or less-than is undefined.
  1168. if((isnan)() || (v.isnan)()) { return ((isnan)() ? 1 : -1); }
  1169. if((isinf)() && (v.isinf)())
  1170. {
  1171. // Both *this and v are infinite. They are equal if they have the same sign.
  1172. // Otherwise, *this is less than v if and only if *this is negative.
  1173. return ((neg == v.neg) ? 0 : (neg ? -1 : 1));
  1174. }
  1175. if((isinf)())
  1176. {
  1177. // *this is infinite, but v is finite.
  1178. // So negative infinite *this is less than any finite v.
  1179. // Whereas positive infinite *this is greater than any finite v.
  1180. return (isneg() ? -1 : 1);
  1181. }
  1182. else
  1183. {
  1184. // *this is finite, and v is infinite.
  1185. // So any finite *this is greater than negative infinite v.
  1186. // Whereas any finite *this is less than positive infinite v.
  1187. return (v.neg ? 1 : -1);
  1188. }
  1189. }
  1190. // And now handle all *finite* cases.
  1191. if(iszero())
  1192. {
  1193. // The value of *this is zero and v is either zero or non-zero.
  1194. return (v.iszero() ? 0
  1195. : (v.neg ? 1 : -1));
  1196. }
  1197. else if(v.iszero())
  1198. {
  1199. // The value of v is zero and *this is non-zero.
  1200. return (neg ? -1 : 1);
  1201. }
  1202. else
  1203. {
  1204. // Both *this and v are non-zero.
  1205. if(neg != v.neg)
  1206. {
  1207. // The signs are different.
  1208. return (neg ? -1 : 1);
  1209. }
  1210. else if(exp != v.exp)
  1211. {
  1212. // The signs are the same and the exponents are different.
  1213. const int val_cexpression = ((exp < v.exp) ? 1 : -1);
  1214. return (neg ? val_cexpression : -val_cexpression);
  1215. }
  1216. else
  1217. {
  1218. // The signs are the same and the exponents are the same.
  1219. // Compare the data.
  1220. const int val_cmp_data = cmp_data(v.data);
  1221. return ((!neg) ? val_cmp_data : -val_cmp_data);
  1222. }
  1223. }
  1224. }
  1225. template <unsigned Digits10, class ExponentType, class Allocator>
  1226. bool cpp_dec_float<Digits10, ExponentType, Allocator>::isone() const
  1227. {
  1228. // Check if the value of *this is identically 1 or very close to 1.
  1229. const bool not_negative_and_is_finite = ((!neg) && (isfinite)());
  1230. if(not_negative_and_is_finite)
  1231. {
  1232. if((data[0u] == static_cast<boost::uint32_t>(1u)) && (exp == static_cast<ExponentType>(0)))
  1233. {
  1234. const typename array_type::const_iterator it_non_zero = std::find_if(data.begin(), data.end(), data_elem_is_non_zero_predicate);
  1235. return (it_non_zero == data.end());
  1236. }
  1237. else if((data[0u] == static_cast<boost::uint32_t>(cpp_dec_float_elem_mask - 1)) && (exp == static_cast<ExponentType>(-cpp_dec_float_elem_digits10)))
  1238. {
  1239. const typename array_type::const_iterator it_non_nine = std::find_if(data.begin(), data.end(), data_elem_is_non_nine_predicate);
  1240. return (it_non_nine == data.end());
  1241. }
  1242. }
  1243. return false;
  1244. }
  1245. template <unsigned Digits10, class ExponentType, class Allocator>
  1246. bool cpp_dec_float<Digits10, ExponentType, Allocator>::isint() const
  1247. {
  1248. if(fpclass != cpp_dec_float_finite) { return false; }
  1249. if(iszero()) { return true; }
  1250. if(exp < static_cast<ExponentType>(0)) { return false; } // |*this| < 1.
  1251. const typename array_type::size_type offset_decimal_part = static_cast<typename array_type::size_type>(exp / cpp_dec_float_elem_digits10) + 1u;
  1252. if(offset_decimal_part >= static_cast<typename array_type::size_type>(cpp_dec_float_elem_number))
  1253. {
  1254. // The number is too large to resolve the integer part.
  1255. // It considered to be a pure integer.
  1256. return true;
  1257. }
  1258. typename array_type::const_iterator it_non_zero = std::find_if(data.begin() + offset_decimal_part, data.end(), data_elem_is_non_zero_predicate);
  1259. return (it_non_zero == data.end());
  1260. }
  1261. template <unsigned Digits10, class ExponentType, class Allocator>
  1262. void cpp_dec_float<Digits10, ExponentType, Allocator>::extract_parts(double& mantissa, ExponentType& exponent) const
  1263. {
  1264. // Extract the approximate parts mantissa and base-10 exponent from the input cpp_dec_float<Digits10, ExponentType, Allocator> value x.
  1265. // Extracts the mantissa and exponent.
  1266. exponent = exp;
  1267. boost::uint32_t p10 = static_cast<boost::uint32_t>(1u);
  1268. boost::uint32_t test = data[0u];
  1269. for(;;)
  1270. {
  1271. test /= static_cast<boost::uint32_t>(10u);
  1272. if(test == static_cast<boost::uint32_t>(0u))
  1273. {
  1274. break;
  1275. }
  1276. p10 *= static_cast<boost::uint32_t>(10u);
  1277. ++exponent;
  1278. }
  1279. // Establish the upper bound of limbs for extracting the double.
  1280. const int max_elem_in_double_count = static_cast<int>(static_cast<boost::int32_t>(std::numeric_limits<double>::digits10) / cpp_dec_float_elem_digits10)
  1281. + (static_cast<int>(static_cast<boost::int32_t>(std::numeric_limits<double>::digits10) % cpp_dec_float_elem_digits10) != 0 ? 1 : 0)
  1282. + 1;
  1283. // And make sure this upper bound stays within bounds of the elems.
  1284. const std::size_t max_elem_extract_count = static_cast<std::size_t>((std::min)(static_cast<boost::int32_t>(max_elem_in_double_count), cpp_dec_float_elem_number));
  1285. // Extract into the mantissa the first limb, extracted as a double.
  1286. mantissa = static_cast<double>(data[0]);
  1287. double scale = 1.0;
  1288. // Extract the rest of the mantissa piecewise from the limbs.
  1289. for(std::size_t i = 1u; i < max_elem_extract_count; i++)
  1290. {
  1291. scale /= static_cast<double>(cpp_dec_float_elem_mask);
  1292. mantissa += (static_cast<double>(data[i]) * scale);
  1293. }
  1294. mantissa /= static_cast<double>(p10);
  1295. if(neg) { mantissa = -mantissa; }
  1296. }
  1297. template <unsigned Digits10, class ExponentType, class Allocator>
  1298. double cpp_dec_float<Digits10, ExponentType, Allocator>::extract_double() const
  1299. {
  1300. // Returns the double conversion of a cpp_dec_float<Digits10, ExponentType, Allocator>.
  1301. // Check for non-normal cpp_dec_float<Digits10, ExponentType, Allocator>.
  1302. if(!(isfinite)())
  1303. {
  1304. if((isnan)())
  1305. {
  1306. return std::numeric_limits<double>::quiet_NaN();
  1307. }
  1308. else
  1309. {
  1310. return ((!neg) ? std::numeric_limits<double>::infinity()
  1311. : -std::numeric_limits<double>::infinity());
  1312. }
  1313. }
  1314. cpp_dec_float<Digits10, ExponentType, Allocator> xx(*this);
  1315. if(xx.isneg())
  1316. xx.negate();
  1317. // Check if *this cpp_dec_float<Digits10, ExponentType, Allocator> is zero.
  1318. if(iszero() || (xx.compare(double_min()) < 0))
  1319. {
  1320. return 0.0;
  1321. }
  1322. // Check if *this cpp_dec_float<Digits10, ExponentType, Allocator> exceeds the maximum of double.
  1323. if(xx.compare(double_max()) > 0)
  1324. {
  1325. return ((!neg) ? std::numeric_limits<double>::infinity()
  1326. : -std::numeric_limits<double>::infinity());
  1327. }
  1328. std::stringstream ss;
  1329. ss << str(std::numeric_limits<double>::digits10 + (2 + 1), std::ios_base::scientific);
  1330. double d;
  1331. ss >> d;
  1332. return d;
  1333. }
  1334. template <unsigned Digits10, class ExponentType, class Allocator>
  1335. long double cpp_dec_float<Digits10, ExponentType, Allocator>::extract_long_double() const
  1336. {
  1337. // Returns the long double conversion of a cpp_dec_float<Digits10, ExponentType, Allocator>.
  1338. // Check if *this cpp_dec_float<Digits10, ExponentType, Allocator> is subnormal.
  1339. if(!(isfinite)())
  1340. {
  1341. if((isnan)())
  1342. {
  1343. return std::numeric_limits<long double>::quiet_NaN();
  1344. }
  1345. else
  1346. {
  1347. return ((!neg) ? std::numeric_limits<long double>::infinity()
  1348. : -std::numeric_limits<long double>::infinity());
  1349. }
  1350. }
  1351. cpp_dec_float<Digits10, ExponentType, Allocator> xx(*this);
  1352. if(xx.isneg())
  1353. xx.negate();
  1354. // Check if *this cpp_dec_float<Digits10, ExponentType, Allocator> is zero.
  1355. if(iszero() || (xx.compare(long_double_min()) < 0))
  1356. {
  1357. return static_cast<long double>(0.0);
  1358. }
  1359. // Check if *this cpp_dec_float<Digits10, ExponentType, Allocator> exceeds the maximum of double.
  1360. if(xx.compare(long_double_max()) > 0)
  1361. {
  1362. return ((!neg) ? std::numeric_limits<long double>::infinity()
  1363. : -std::numeric_limits<long double>::infinity());
  1364. }
  1365. std::stringstream ss;
  1366. ss << str(std::numeric_limits<long double>::digits10 + (2 + 1), std::ios_base::scientific);
  1367. long double ld;
  1368. ss >> ld;
  1369. return ld;
  1370. }
  1371. template <unsigned Digits10, class ExponentType, class Allocator>
  1372. boost::long_long_type cpp_dec_float<Digits10, ExponentType, Allocator>::extract_signed_long_long() const
  1373. {
  1374. // Extracts a signed long long from *this.
  1375. // If (x > maximum of long long) or (x < minimum of long long),
  1376. // then the maximum or minimum of long long is returned accordingly.
  1377. if(exp < static_cast<ExponentType>(0))
  1378. {
  1379. return static_cast<boost::long_long_type>(0);
  1380. }
  1381. const bool b_neg = isneg();
  1382. boost::ulong_long_type val;
  1383. if((!b_neg) && (compare(long_long_max()) > 0))
  1384. {
  1385. return (std::numeric_limits<boost::long_long_type>::max)();
  1386. }
  1387. else if(b_neg && (compare(long_long_min()) < 0))
  1388. {
  1389. return (std::numeric_limits<boost::long_long_type>::min)();
  1390. }
  1391. else
  1392. {
  1393. // Extract the data into an boost::ulong_long_type value.
  1394. cpp_dec_float<Digits10, ExponentType, Allocator> xn(extract_integer_part());
  1395. if(xn.isneg())
  1396. xn.negate();
  1397. val = static_cast<boost::ulong_long_type>(xn.data[0]);
  1398. const boost::int32_t imax = (std::min)(static_cast<boost::int32_t>(static_cast<boost::int32_t>(xn.exp) / cpp_dec_float_elem_digits10), static_cast<boost::int32_t>(cpp_dec_float_elem_number - static_cast<boost::int32_t>(1)));
  1399. for(boost::int32_t i = static_cast<boost::int32_t>(1); i <= imax; i++)
  1400. {
  1401. val *= static_cast<boost::ulong_long_type>(cpp_dec_float_elem_mask);
  1402. val += static_cast<boost::ulong_long_type>(xn.data[i]);
  1403. }
  1404. }
  1405. if (!b_neg)
  1406. {
  1407. return static_cast<boost::long_long_type>(val);
  1408. }
  1409. else
  1410. {
  1411. // This strange expression avoids a hardware trap in the corner case
  1412. // that val is the most negative value permitted in boost::long_long_type.
  1413. // See https://svn.boost.org/trac/boost/ticket/9740.
  1414. //
  1415. boost::long_long_type sval = static_cast<boost::long_long_type>(val - 1);
  1416. sval = -sval;
  1417. --sval;
  1418. return sval;
  1419. }
  1420. }
  1421. template <unsigned Digits10, class ExponentType, class Allocator>
  1422. boost::ulong_long_type cpp_dec_float<Digits10, ExponentType, Allocator>::extract_unsigned_long_long() const
  1423. {
  1424. // Extracts an boost::ulong_long_type from *this.
  1425. // If x exceeds the maximum of boost::ulong_long_type,
  1426. // then the maximum of boost::ulong_long_type is returned.
  1427. // If x is negative, then the boost::ulong_long_type cast of
  1428. // the long long extracted value is returned.
  1429. if(isneg())
  1430. {
  1431. return static_cast<boost::ulong_long_type>(extract_signed_long_long());
  1432. }
  1433. if(exp < static_cast<ExponentType>(0))
  1434. {
  1435. return static_cast<boost::ulong_long_type>(0u);
  1436. }
  1437. const cpp_dec_float<Digits10, ExponentType, Allocator> xn(extract_integer_part());
  1438. boost::ulong_long_type val;
  1439. if(xn.compare(ulong_long_max()) > 0)
  1440. {
  1441. return (std::numeric_limits<boost::ulong_long_type>::max)();
  1442. }
  1443. else
  1444. {
  1445. // Extract the data into an boost::ulong_long_type value.
  1446. val = static_cast<boost::ulong_long_type>(xn.data[0]);
  1447. const boost::int32_t imax = (std::min)(static_cast<boost::int32_t>(static_cast<boost::int32_t>(xn.exp) / cpp_dec_float_elem_digits10), static_cast<boost::int32_t>(cpp_dec_float_elem_number - static_cast<boost::int32_t>(1)));
  1448. for(boost::int32_t i = static_cast<boost::int32_t>(1); i <= imax; i++)
  1449. {
  1450. val *= static_cast<boost::ulong_long_type>(cpp_dec_float_elem_mask);
  1451. val += static_cast<boost::ulong_long_type>(xn.data[i]);
  1452. }
  1453. }
  1454. return val;
  1455. }
  1456. template <unsigned Digits10, class ExponentType, class Allocator>
  1457. cpp_dec_float<Digits10, ExponentType, Allocator> cpp_dec_float<Digits10, ExponentType, Allocator>::extract_integer_part() const
  1458. {
  1459. // Compute the signed integer part of x.
  1460. if(!(isfinite)())
  1461. {
  1462. return *this;
  1463. }
  1464. if(exp < static_cast<ExponentType>(0))
  1465. {
  1466. // The absolute value of the number is smaller than 1.
  1467. // Thus the integer part is zero.
  1468. return zero();
  1469. }
  1470. // Truncate the digits from the decimal part, including guard digits
  1471. // that do not belong to the integer part.
  1472. // Make a local copy.
  1473. cpp_dec_float<Digits10, ExponentType, Allocator> x = *this;
  1474. // Clear out the decimal portion
  1475. const size_t first_clear = (static_cast<size_t>(x.exp) / static_cast<size_t>(cpp_dec_float_elem_digits10)) + 1u;
  1476. const size_t last_clear = static_cast<size_t>(cpp_dec_float_elem_number);
  1477. if(first_clear < last_clear)
  1478. std::fill(x.data.begin() + first_clear, x.data.begin() + last_clear, static_cast<boost::uint32_t>(0u));
  1479. return x;
  1480. }
  1481. template <unsigned Digits10, class ExponentType, class Allocator>
  1482. std::string cpp_dec_float<Digits10, ExponentType, Allocator>::str(boost::intmax_t number_of_digits, std::ios_base::fmtflags f) const
  1483. {
  1484. if((this->isinf)())
  1485. {
  1486. if(this->isneg())
  1487. return "-inf";
  1488. else if(f & std::ios_base::showpos)
  1489. return "+inf";
  1490. else
  1491. return "inf";
  1492. }
  1493. else if((this->isnan)())
  1494. {
  1495. return "nan";
  1496. }
  1497. std::string str;
  1498. boost::intmax_t org_digits(number_of_digits);
  1499. ExponentType my_exp = order();
  1500. if(number_of_digits == 0)
  1501. number_of_digits = cpp_dec_float_total_digits10;
  1502. if(f & std::ios_base::fixed)
  1503. {
  1504. number_of_digits += my_exp + 1;
  1505. }
  1506. else if(f & std::ios_base::scientific)
  1507. ++number_of_digits;
  1508. // Determine the number of elements needed to provide the requested digits from cpp_dec_float<Digits10, ExponentType, Allocator>.
  1509. const std::size_t number_of_elements = (std::min)(static_cast<std::size_t>((number_of_digits / static_cast<std::size_t>(cpp_dec_float_elem_digits10)) + 2u),
  1510. static_cast<std::size_t>(cpp_dec_float_elem_number));
  1511. // Extract the remaining digits from cpp_dec_float<Digits10, ExponentType, Allocator> after the decimal point.
  1512. str = boost::lexical_cast<std::string>(data[0]);
  1513. // Extract all of the digits from cpp_dec_float<Digits10, ExponentType, Allocator>, beginning with the first data element.
  1514. for(std::size_t i = static_cast<std::size_t>(1u); i < number_of_elements; i++)
  1515. {
  1516. std::stringstream ss;
  1517. ss << std::setw(static_cast<std::streamsize>(cpp_dec_float_elem_digits10))
  1518. << std::setfill(static_cast<char>('0'))
  1519. << data[i];
  1520. str += ss.str();
  1521. }
  1522. bool have_leading_zeros = false;
  1523. if(number_of_digits == 0)
  1524. {
  1525. // We only get here if the output format is "fixed" and we just need to
  1526. // round the first non-zero digit.
  1527. number_of_digits -= my_exp + 1; // reset to original value
  1528. str.insert(static_cast<std::string::size_type>(0), std::string::size_type(number_of_digits), '0');
  1529. have_leading_zeros = true;
  1530. }
  1531. if(number_of_digits < 0)
  1532. {
  1533. str = "0";
  1534. if(isneg())
  1535. str.insert(static_cast<std::string::size_type>(0), 1, '-');
  1536. boost::multiprecision::detail::format_float_string(str, 0, number_of_digits - my_exp - 1, f, this->iszero());
  1537. return str;
  1538. }
  1539. else
  1540. {
  1541. // Cut the output to the size of the precision.
  1542. if(str.length() > static_cast<std::string::size_type>(number_of_digits))
  1543. {
  1544. // Get the digit after the last needed digit for rounding
  1545. const boost::uint32_t round = static_cast<boost::uint32_t>(static_cast<boost::uint32_t>(str[static_cast<std::string::size_type>(number_of_digits)]) - static_cast<boost::uint32_t>('0'));
  1546. bool need_round_up = round >= 5u;
  1547. if(round == 5u)
  1548. {
  1549. const boost::uint32_t ix = static_cast<boost::uint32_t>(static_cast<boost::uint32_t>(str[static_cast<std::string::size_type>(number_of_digits - 1)]) - static_cast<boost::uint32_t>('0'));
  1550. if((ix & 1u) == 0)
  1551. {
  1552. // We have an even digit followed by a 5, so we might not actually need to round up
  1553. // if all the remaining digits are zero:
  1554. if(str.find_first_not_of('0', static_cast<std::string::size_type>(number_of_digits + 1)) == std::string::npos)
  1555. {
  1556. bool all_zeros = true;
  1557. // No none-zero trailing digits in the string, now check whatever parts we didn't convert to the string:
  1558. for(std::size_t i = number_of_elements; i < data.size(); i++)
  1559. {
  1560. if(data[i])
  1561. {
  1562. all_zeros = false;
  1563. break;
  1564. }
  1565. }
  1566. if(all_zeros)
  1567. need_round_up = false; // tie break - round to even.
  1568. }
  1569. }
  1570. }
  1571. // Truncate the string
  1572. str.erase(static_cast<std::string::size_type>(number_of_digits));
  1573. if(need_round_up)
  1574. {
  1575. std::size_t ix = static_cast<std::size_t>(str.length() - 1u);
  1576. // Every trailing 9 must be rounded up
  1577. while(ix && (static_cast<boost::int32_t>(str.at(ix)) - static_cast<boost::int32_t>('0') == static_cast<boost::int32_t>(9)))
  1578. {
  1579. str.at(ix) = static_cast<char>('0');
  1580. --ix;
  1581. }
  1582. if(!ix)
  1583. {
  1584. // There were nothing but trailing nines.
  1585. if(static_cast<boost::int32_t>(static_cast<boost::int32_t>(str.at(ix)) - static_cast<boost::int32_t>(0x30)) == static_cast<boost::int32_t>(9))
  1586. {
  1587. // Increment up to the next order and adjust exponent.
  1588. str.at(ix) = static_cast<char>('1');
  1589. ++my_exp;
  1590. }
  1591. else
  1592. {
  1593. // Round up this digit.
  1594. ++str.at(ix);
  1595. }
  1596. }
  1597. else
  1598. {
  1599. // Round up the last digit.
  1600. ++str[ix];
  1601. }
  1602. }
  1603. }
  1604. }
  1605. if(have_leading_zeros)
  1606. {
  1607. // We need to take the zeros back out again, and correct the exponent
  1608. // if we rounded up:
  1609. if(str[std::string::size_type(number_of_digits - 1)] != '0')
  1610. {
  1611. ++my_exp;
  1612. str.erase(0, std::string::size_type(number_of_digits - 1));
  1613. }
  1614. else
  1615. str.erase(0, std::string::size_type(number_of_digits));
  1616. }
  1617. if(isneg())
  1618. str.insert(static_cast<std::string::size_type>(0), 1, '-');
  1619. boost::multiprecision::detail::format_float_string(str, my_exp, org_digits, f, this->iszero());
  1620. return str;
  1621. }
  1622. template <unsigned Digits10, class ExponentType, class Allocator>
  1623. bool cpp_dec_float<Digits10, ExponentType, Allocator>::rd_string(const char* const s)
  1624. {
  1625. #ifndef BOOST_NO_EXCEPTIONS
  1626. try{
  1627. #endif
  1628. std::string str(s);
  1629. // TBD: Using several regular expressions may significantly reduce
  1630. // the code complexity (and perhaps the run-time) of rd_string().
  1631. // Get a possible exponent and remove it.
  1632. exp = static_cast<ExponentType>(0);
  1633. std::size_t pos;
  1634. if( ((pos = str.find('e')) != std::string::npos)
  1635. || ((pos = str.find('E')) != std::string::npos)
  1636. )
  1637. {
  1638. // Remove the exponent part from the string.
  1639. exp = boost::lexical_cast<ExponentType>(static_cast<const char*>(str.c_str() + (pos + 1u)));
  1640. str = str.substr(static_cast<std::size_t>(0u), pos);
  1641. }
  1642. // Get a possible +/- sign and remove it.
  1643. neg = false;
  1644. if(str.size())
  1645. {
  1646. if(str[0] == '-')
  1647. {
  1648. neg = true;
  1649. str.erase(0, 1);
  1650. }
  1651. else if(str[0] == '+')
  1652. {
  1653. str.erase(0, 1);
  1654. }
  1655. }
  1656. //
  1657. // Special cases for infinities and NaN's:
  1658. //
  1659. if((str == "inf") || (str == "INF") || (str == "infinity") || (str == "INFINITY"))
  1660. {
  1661. if(neg)
  1662. {
  1663. *this = this->inf();
  1664. this->negate();
  1665. }
  1666. else
  1667. *this = this->inf();
  1668. return true;
  1669. }
  1670. if((str.size() >= 3) && ((str.substr(0, 3) == "nan") || (str.substr(0, 3) == "NAN") || (str.substr(0, 3) == "NaN")))
  1671. {
  1672. *this = this->nan();
  1673. return true;
  1674. }
  1675. // Remove the leading zeros for all input types.
  1676. const std::string::iterator fwd_it_leading_zero = std::find_if(str.begin(), str.end(), char_is_nonzero_predicate);
  1677. if(fwd_it_leading_zero != str.begin())
  1678. {
  1679. if(fwd_it_leading_zero == str.end())
  1680. {
  1681. // The string contains nothing but leading zeros.
  1682. // This string represents zero.
  1683. operator=(zero());
  1684. return true;
  1685. }
  1686. else
  1687. {
  1688. str.erase(str.begin(), fwd_it_leading_zero);
  1689. }
  1690. }
  1691. // Put the input string into the standard cpp_dec_float<Digits10, ExponentType, Allocator> input form
  1692. // aaa.bbbbE+/-n, where aaa has 1...cpp_dec_float_elem_digits10, bbbb has an
  1693. // even multiple of cpp_dec_float_elem_digits10 which are possibly zero padded
  1694. // on the right-end, and n is a signed 64-bit integer which is an
  1695. // even multiple of cpp_dec_float_elem_digits10.
  1696. // Find a possible decimal point.
  1697. pos = str.find(static_cast<char>('.'));
  1698. if(pos != std::string::npos)
  1699. {
  1700. // Remove all trailing insignificant zeros.
  1701. const std::string::const_reverse_iterator rit_non_zero = std::find_if(str.rbegin(), str.rend(), char_is_nonzero_predicate);
  1702. if(rit_non_zero != static_cast<std::string::const_reverse_iterator>(str.rbegin()))
  1703. {
  1704. const std::string::size_type ofs = str.length() - std::distance<std::string::const_reverse_iterator>(str.rbegin(), rit_non_zero);
  1705. str.erase(str.begin() + ofs, str.end());
  1706. }
  1707. // Check if the input is identically zero.
  1708. if(str == std::string("."))
  1709. {
  1710. operator=(zero());
  1711. return true;
  1712. }
  1713. // Remove leading significant zeros just after the decimal point
  1714. // and adjust the exponent accordingly.
  1715. // Note that the while-loop operates only on strings of the form ".000abcd..."
  1716. // and peels away the zeros just after the decimal point.
  1717. if(str.at(static_cast<std::size_t>(0u)) == static_cast<char>('.'))
  1718. {
  1719. const std::string::iterator it_non_zero = std::find_if(str.begin() + 1u, str.end(), char_is_nonzero_predicate);
  1720. std::size_t delta_exp = static_cast<std::size_t>(0u);
  1721. if(str.at(static_cast<std::size_t>(1u)) == static_cast<char>('0'))
  1722. {
  1723. delta_exp = std::distance<std::string::const_iterator>(str.begin() + 1u, it_non_zero);
  1724. }
  1725. // Bring one single digit into the mantissa and adjust the exponent accordingly.
  1726. str.erase(str.begin(), it_non_zero);
  1727. str.insert(static_cast<std::string::size_type>(1u), ".");
  1728. exp -= static_cast<ExponentType>(delta_exp + 1u);
  1729. }
  1730. }
  1731. else
  1732. {
  1733. // Input string has no decimal point: Append decimal point.
  1734. str.append(".");
  1735. }
  1736. // Shift the decimal point such that the exponent is an even multiple of cpp_dec_float_elem_digits10.
  1737. std::size_t n_shift = static_cast<std::size_t>(0u);
  1738. const std::size_t n_exp_rem = static_cast<std::size_t>(exp % static_cast<ExponentType>(cpp_dec_float_elem_digits10));
  1739. if((exp % static_cast<ExponentType>(cpp_dec_float_elem_digits10)) != static_cast<ExponentType>(0))
  1740. {
  1741. n_shift = ((exp < static_cast<ExponentType>(0))
  1742. ? static_cast<std::size_t>(n_exp_rem + static_cast<std::size_t>(cpp_dec_float_elem_digits10))
  1743. : static_cast<std::size_t>(n_exp_rem));
  1744. }
  1745. // Make sure that there are enough digits for the decimal point shift.
  1746. pos = str.find(static_cast<char>('.'));
  1747. std::size_t pos_plus_one = static_cast<std::size_t>(pos + 1u);
  1748. if((str.length() - pos_plus_one) < n_shift)
  1749. {
  1750. const std::size_t sz = static_cast<std::size_t>(n_shift - (str.length() - pos_plus_one));
  1751. str.append(std::string(sz, static_cast<char>('0')));
  1752. }
  1753. // Do the decimal point shift.
  1754. if(n_shift != static_cast<std::size_t>(0u))
  1755. {
  1756. str.insert(static_cast<std::string::size_type>(pos_plus_one + n_shift), ".");
  1757. str.erase(pos, static_cast<std::string::size_type>(1u));
  1758. exp -= static_cast<ExponentType>(n_shift);
  1759. }
  1760. // Cut the size of the mantissa to <= cpp_dec_float_elem_digits10.
  1761. pos = str.find(static_cast<char>('.'));
  1762. pos_plus_one = static_cast<std::size_t>(pos + 1u);
  1763. if(pos > static_cast<std::size_t>(cpp_dec_float_elem_digits10))
  1764. {
  1765. const boost::int32_t n_pos = static_cast<boost::int32_t>(pos);
  1766. const boost::int32_t n_rem_is_zero = ((static_cast<boost::int32_t>(n_pos % cpp_dec_float_elem_digits10) == static_cast<boost::int32_t>(0)) ? static_cast<boost::int32_t>(1) : static_cast<boost::int32_t>(0));
  1767. const boost::int32_t n = static_cast<boost::int32_t>(static_cast<boost::int32_t>(n_pos / cpp_dec_float_elem_digits10) - n_rem_is_zero);
  1768. str.insert(static_cast<std::size_t>(static_cast<boost::int32_t>(n_pos - static_cast<boost::int32_t>(n * cpp_dec_float_elem_digits10))), ".");
  1769. str.erase(pos_plus_one, static_cast<std::size_t>(1u));
  1770. exp += static_cast<ExponentType>(static_cast<ExponentType>(n) * static_cast<ExponentType>(cpp_dec_float_elem_digits10));
  1771. }
  1772. // Pad the decimal part such that its value is an even
  1773. // multiple of cpp_dec_float_elem_digits10.
  1774. pos = str.find(static_cast<char>('.'));
  1775. pos_plus_one = static_cast<std::size_t>(pos + 1u);
  1776. const boost::int32_t n_dec = static_cast<boost::int32_t>(static_cast<boost::int32_t>(str.length() - 1u) - static_cast<boost::int32_t>(pos));
  1777. const boost::int32_t n_rem = static_cast<boost::int32_t>(n_dec % cpp_dec_float_elem_digits10);
  1778. boost::int32_t n_cnt = ((n_rem != static_cast<boost::int32_t>(0))
  1779. ? static_cast<boost::int32_t>(cpp_dec_float_elem_digits10 - n_rem)
  1780. : static_cast<boost::int32_t>(0));
  1781. if(n_cnt != static_cast<boost::int32_t>(0))
  1782. {
  1783. str.append(static_cast<std::size_t>(n_cnt), static_cast<char>('0'));
  1784. }
  1785. // Truncate decimal part if it is too long.
  1786. const std::size_t max_dec = static_cast<std::size_t>((cpp_dec_float_elem_number - 1) * cpp_dec_float_elem_digits10);
  1787. if(static_cast<std::size_t>(str.length() - pos) > max_dec)
  1788. {
  1789. str = str.substr(static_cast<std::size_t>(0u),
  1790. static_cast<std::size_t>(pos_plus_one + max_dec));
  1791. }
  1792. // Now the input string has the standard cpp_dec_float<Digits10, ExponentType, Allocator> input form.
  1793. // (See the comment above.)
  1794. // Set all the data elements to 0.
  1795. std::fill(data.begin(), data.end(), static_cast<boost::uint32_t>(0u));
  1796. // Extract the data.
  1797. // First get the digits to the left of the decimal point...
  1798. data[0u] = boost::lexical_cast<boost::uint32_t>(str.substr(static_cast<std::size_t>(0u), pos));
  1799. // ...then get the remaining digits to the right of the decimal point.
  1800. const std::string::size_type i_end = ((str.length() - pos_plus_one) / static_cast<std::string::size_type>(cpp_dec_float_elem_digits10));
  1801. for(std::string::size_type i = static_cast<std::string::size_type>(0u); i < i_end; i++)
  1802. {
  1803. const std::string::const_iterator it = str.begin()
  1804. + pos_plus_one
  1805. + (i * static_cast<std::string::size_type>(cpp_dec_float_elem_digits10));
  1806. data[i + 1u] = boost::lexical_cast<boost::uint32_t>(std::string(it, it + static_cast<std::string::size_type>(cpp_dec_float_elem_digits10)));
  1807. }
  1808. // Check for overflow...
  1809. if(exp > cpp_dec_float_max_exp10)
  1810. {
  1811. const bool b_result_is_neg = neg;
  1812. *this = inf();
  1813. if(b_result_is_neg)
  1814. negate();
  1815. }
  1816. // ...and check for underflow.
  1817. if(exp <= cpp_dec_float_min_exp10)
  1818. {
  1819. if(exp == cpp_dec_float_min_exp10)
  1820. {
  1821. // Check for identity with the minimum value.
  1822. cpp_dec_float<Digits10, ExponentType, Allocator> test = *this;
  1823. test.exp = static_cast<ExponentType>(0);
  1824. if(test.isone())
  1825. {
  1826. *this = zero();
  1827. }
  1828. }
  1829. else
  1830. {
  1831. *this = zero();
  1832. }
  1833. }
  1834. #ifndef BOOST_NO_EXCEPTIONS
  1835. }
  1836. catch(const bad_lexical_cast&)
  1837. {
  1838. // Rethrow with better error message:
  1839. std::string msg = "Unable to parse the string \"";
  1840. msg += s;
  1841. msg += "\" as a floating point value.";
  1842. throw std::runtime_error(msg);
  1843. }
  1844. #endif
  1845. return true;
  1846. }
  1847. template <unsigned Digits10, class ExponentType, class Allocator>
  1848. cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float(const double mantissa, const ExponentType exponent)
  1849. : data (),
  1850. exp (static_cast<ExponentType>(0)),
  1851. neg (false),
  1852. fpclass (cpp_dec_float_finite),
  1853. prec_elem(cpp_dec_float_elem_number)
  1854. {
  1855. // Create *this cpp_dec_float<Digits10, ExponentType, Allocator> from a given mantissa and exponent.
  1856. // Note: This constructor does not maintain the full precision of double.
  1857. const bool mantissa_is_iszero = (::fabs(mantissa) < ((std::numeric_limits<double>::min)() * (1.0 + std::numeric_limits<double>::epsilon())));
  1858. if(mantissa_is_iszero)
  1859. {
  1860. std::fill(data.begin(), data.end(), static_cast<boost::uint32_t>(0u));
  1861. return;
  1862. }
  1863. const bool b_neg = (mantissa < 0.0);
  1864. double d = ((!b_neg) ? mantissa : -mantissa);
  1865. ExponentType e = exponent;
  1866. while(d > 10.0) { d /= 10.0; ++e; }
  1867. while(d < 1.0) { d *= 10.0; --e; }
  1868. boost::int32_t shift = static_cast<boost::int32_t>(e % static_cast<boost::int32_t>(cpp_dec_float_elem_digits10));
  1869. while(static_cast<boost::int32_t>(shift-- % cpp_dec_float_elem_digits10) != static_cast<boost::int32_t>(0))
  1870. {
  1871. d *= 10.0;
  1872. --e;
  1873. }
  1874. exp = e;
  1875. neg = b_neg;
  1876. std::fill(data.begin(), data.end(), static_cast<boost::uint32_t>(0u));
  1877. static const boost::int32_t digit_ratio = static_cast<boost::int32_t>(static_cast<boost::int32_t>(std::numeric_limits<double>::digits10) / static_cast<boost::int32_t>(cpp_dec_float_elem_digits10));
  1878. static const boost::int32_t digit_loops = static_cast<boost::int32_t>(digit_ratio + static_cast<boost::int32_t>(2));
  1879. for(boost::int32_t i = static_cast<boost::int32_t>(0); i < digit_loops; i++)
  1880. {
  1881. boost::uint32_t n = static_cast<boost::uint32_t>(static_cast<boost::uint64_t>(d));
  1882. data[i] = static_cast<boost::uint32_t>(n);
  1883. d -= static_cast<double>(n);
  1884. d *= static_cast<double>(cpp_dec_float_elem_mask);
  1885. }
  1886. }
  1887. template <unsigned Digits10, class ExponentType, class Allocator>
  1888. cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::operator= (long double a)
  1889. {
  1890. // Christopher Kormanyos's original code used a cast to boost::long_long_type here, but that fails
  1891. // when long double has more digits than a boost::long_long_type.
  1892. using std::frexp;
  1893. using std::ldexp;
  1894. using std::floor;
  1895. if(a == 0)
  1896. return *this = zero();
  1897. if(a == 1)
  1898. return *this = one();
  1899. if((boost::math::isinf)(a))
  1900. {
  1901. *this = inf();
  1902. if(a < 0)
  1903. this->negate();
  1904. return *this;
  1905. }
  1906. if((boost::math::isnan)(a))
  1907. return *this = nan();
  1908. int e;
  1909. long double f, term;
  1910. *this = zero();
  1911. f = frexp(a, &e);
  1912. // See https://svn.boost.org/trac/boost/ticket/10924 for an example of why this may go wrong:
  1913. BOOST_ASSERT((boost::math::isfinite)(f));
  1914. static const int shift = std::numeric_limits<int>::digits - 1;
  1915. while(f)
  1916. {
  1917. // extract int sized bits from f:
  1918. f = ldexp(f, shift);
  1919. BOOST_ASSERT((boost::math::isfinite)(f));
  1920. term = floor(f);
  1921. e -= shift;
  1922. *this *= pow2(shift);
  1923. if(term > 0)
  1924. add_unsigned_long_long(static_cast<unsigned>(term));
  1925. else
  1926. sub_unsigned_long_long(static_cast<unsigned>(-term));
  1927. f -= term;
  1928. }
  1929. if(e != 0)
  1930. *this *= pow2(e);
  1931. return *this;
  1932. }
  1933. template <unsigned Digits10, class ExponentType, class Allocator>
  1934. void cpp_dec_float<Digits10, ExponentType, Allocator>::from_unsigned_long_long(const boost::ulong_long_type u)
  1935. {
  1936. std::fill(data.begin(), data.end(), static_cast<boost::uint32_t>(0u));
  1937. exp = static_cast<ExponentType>(0);
  1938. neg = false;
  1939. fpclass = cpp_dec_float_finite;
  1940. prec_elem = cpp_dec_float_elem_number;
  1941. if(u == 0)
  1942. {
  1943. return;
  1944. }
  1945. std::size_t i =static_cast<std::size_t>(0u);
  1946. boost::ulong_long_type uu = u;
  1947. boost::uint32_t temp[(std::numeric_limits<boost::ulong_long_type>::digits10 / static_cast<int>(cpp_dec_float_elem_digits10)) + 3] = { static_cast<boost::uint32_t>(0u) };
  1948. while(uu != static_cast<boost::ulong_long_type>(0u))
  1949. {
  1950. temp[i] = static_cast<boost::uint32_t>(uu % static_cast<boost::ulong_long_type>(cpp_dec_float_elem_mask));
  1951. uu = static_cast<boost::ulong_long_type>(uu / static_cast<boost::ulong_long_type>(cpp_dec_float_elem_mask));
  1952. ++i;
  1953. }
  1954. if(i > static_cast<std::size_t>(1u))
  1955. {
  1956. exp += static_cast<ExponentType>((i - 1u) * static_cast<std::size_t>(cpp_dec_float_elem_digits10));
  1957. }
  1958. std::reverse(temp, temp + i);
  1959. std::copy(temp, temp + (std::min)(i, static_cast<std::size_t>(cpp_dec_float_elem_number)), data.begin());
  1960. }
  1961. template <unsigned Digits10, class ExponentType, class Allocator>
  1962. boost::uint32_t cpp_dec_float<Digits10, ExponentType, Allocator>::mul_loop_uv(boost::uint32_t* const u, const boost::uint32_t* const v, const boost::int32_t p)
  1963. {
  1964. //
  1965. // There is a limit on how many limbs this algorithm can handle without dropping digits
  1966. // due to overflow in the carry, it is:
  1967. //
  1968. // FLOOR( (2^64 - 1) / (10^8 * 10^8) ) == 1844
  1969. //
  1970. BOOST_STATIC_ASSERT_MSG(cpp_dec_float_elem_number < 1800, "Too many limbs in the data type for the multiplication algorithm - unsupported precision in cpp_dec_float.");
  1971. boost::uint64_t carry = static_cast<boost::uint64_t>(0u);
  1972. for(boost::int32_t j = static_cast<boost::int32_t>(p - 1u); j >= static_cast<boost::int32_t>(0); j--)
  1973. {
  1974. boost::uint64_t sum = carry;
  1975. for(boost::int32_t i = j; i >= static_cast<boost::int32_t>(0); i--)
  1976. {
  1977. sum += static_cast<boost::uint64_t>(u[j - i] * static_cast<boost::uint64_t>(v[i]));
  1978. }
  1979. u[j] = static_cast<boost::uint32_t>(sum % static_cast<boost::uint32_t>(cpp_dec_float_elem_mask));
  1980. carry = static_cast<boost::uint64_t>(sum / static_cast<boost::uint32_t>(cpp_dec_float_elem_mask));
  1981. }
  1982. return static_cast<boost::uint32_t>(carry);
  1983. }
  1984. template <unsigned Digits10, class ExponentType, class Allocator>
  1985. boost::uint32_t cpp_dec_float<Digits10, ExponentType, Allocator>::mul_loop_n(boost::uint32_t* const u, boost::uint32_t n, const boost::int32_t p)
  1986. {
  1987. boost::uint64_t carry = static_cast<boost::uint64_t>(0u);
  1988. // Multiplication loop.
  1989. for(boost::int32_t j = p - 1; j >= static_cast<boost::int32_t>(0); j--)
  1990. {
  1991. const boost::uint64_t t = static_cast<boost::uint64_t>(carry + static_cast<boost::uint64_t>(u[j] * static_cast<boost::uint64_t>(n)));
  1992. carry = static_cast<boost::uint64_t>(t / static_cast<boost::uint32_t>(cpp_dec_float_elem_mask));
  1993. u[j] = static_cast<boost::uint32_t>(t - static_cast<boost::uint64_t>(static_cast<boost::uint32_t>(cpp_dec_float_elem_mask) * static_cast<boost::uint64_t>(carry)));
  1994. }
  1995. return static_cast<boost::uint32_t>(carry);
  1996. }
  1997. template <unsigned Digits10, class ExponentType, class Allocator>
  1998. boost::uint32_t cpp_dec_float<Digits10, ExponentType, Allocator>::div_loop_n(boost::uint32_t* const u, boost::uint32_t n, const boost::int32_t p)
  1999. {
  2000. boost::uint64_t prev = static_cast<boost::uint64_t>(0u);
  2001. for(boost::int32_t j = static_cast<boost::int32_t>(0); j < p; j++)
  2002. {
  2003. const boost::uint64_t t = static_cast<boost::uint64_t>(u[j] + static_cast<boost::uint64_t>(prev * static_cast<boost::uint32_t>(cpp_dec_float_elem_mask)));
  2004. u[j] = static_cast<boost::uint32_t>(t / n);
  2005. prev = static_cast<boost::uint64_t>(t - static_cast<boost::uint64_t>(n * static_cast<boost::uint64_t>(u[j])));
  2006. }
  2007. return static_cast<boost::uint32_t>(prev);
  2008. }
  2009. template <unsigned Digits10, class ExponentType, class Allocator>
  2010. cpp_dec_float<Digits10, ExponentType, Allocator> cpp_dec_float<Digits10, ExponentType, Allocator>::pow2(const boost::long_long_type p)
  2011. {
  2012. // Create a static const table of p^2 for -128 < p < +128.
  2013. // Note: The size of this table must be odd-numbered and
  2014. // symmetric about 0.
  2015. init.do_nothing();
  2016. static const boost::array<cpp_dec_float<Digits10, ExponentType, Allocator>, 255u> p2_data =
  2017. {{
  2018. cpp_dec_float("5.877471754111437539843682686111228389093327783860437607543758531392086297273635864257812500000000000e-39"),
  2019. cpp_dec_float("1.175494350822287507968736537222245677818665556772087521508751706278417259454727172851562500000000000e-38"),
  2020. cpp_dec_float("2.350988701644575015937473074444491355637331113544175043017503412556834518909454345703125000000000000e-38"),
  2021. cpp_dec_float("4.701977403289150031874946148888982711274662227088350086035006825113669037818908691406250000000000000e-38"),
  2022. cpp_dec_float("9.403954806578300063749892297777965422549324454176700172070013650227338075637817382812500000000000000e-38"),
  2023. cpp_dec_float("1.880790961315660012749978459555593084509864890835340034414002730045467615127563476562500000000000000e-37"),
  2024. cpp_dec_float("3.761581922631320025499956919111186169019729781670680068828005460090935230255126953125000000000000000e-37"),
  2025. cpp_dec_float("7.523163845262640050999913838222372338039459563341360137656010920181870460510253906250000000000000000e-37"),
  2026. cpp_dec_float("1.504632769052528010199982767644474467607891912668272027531202184036374092102050781250000000000000000e-36"),
  2027. cpp_dec_float("3.009265538105056020399965535288948935215783825336544055062404368072748184204101562500000000000000000e-36"),
  2028. cpp_dec_float("6.018531076210112040799931070577897870431567650673088110124808736145496368408203125000000000000000000e-36"),
  2029. cpp_dec_float("1.203706215242022408159986214115579574086313530134617622024961747229099273681640625000000000000000000e-35"),
  2030. cpp_dec_float("2.407412430484044816319972428231159148172627060269235244049923494458198547363281250000000000000000000e-35"),
  2031. cpp_dec_float("4.814824860968089632639944856462318296345254120538470488099846988916397094726562500000000000000000000e-35"),
  2032. cpp_dec_float("9.629649721936179265279889712924636592690508241076940976199693977832794189453125000000000000000000000e-35"),
  2033. cpp_dec_float("1.925929944387235853055977942584927318538101648215388195239938795566558837890625000000000000000000000e-34"),
  2034. cpp_dec_float("3.851859888774471706111955885169854637076203296430776390479877591133117675781250000000000000000000000e-34"),
  2035. cpp_dec_float("7.703719777548943412223911770339709274152406592861552780959755182266235351562500000000000000000000000e-34"),
  2036. cpp_dec_float("1.540743955509788682444782354067941854830481318572310556191951036453247070312500000000000000000000000e-33"),
  2037. cpp_dec_float("3.081487911019577364889564708135883709660962637144621112383902072906494140625000000000000000000000000e-33"),
  2038. cpp_dec_float("6.162975822039154729779129416271767419321925274289242224767804145812988281250000000000000000000000000e-33"),
  2039. cpp_dec_float("1.232595164407830945955825883254353483864385054857848444953560829162597656250000000000000000000000000e-32"),
  2040. cpp_dec_float("2.465190328815661891911651766508706967728770109715696889907121658325195312500000000000000000000000000e-32"),
  2041. cpp_dec_float("4.930380657631323783823303533017413935457540219431393779814243316650390625000000000000000000000000000e-32"),
  2042. cpp_dec_float("9.860761315262647567646607066034827870915080438862787559628486633300781250000000000000000000000000000e-32"),
  2043. cpp_dec_float("1.972152263052529513529321413206965574183016087772557511925697326660156250000000000000000000000000000e-31"),
  2044. cpp_dec_float("3.944304526105059027058642826413931148366032175545115023851394653320312500000000000000000000000000000e-31"),
  2045. cpp_dec_float("7.888609052210118054117285652827862296732064351090230047702789306640625000000000000000000000000000000e-31"),
  2046. cpp_dec_float("1.577721810442023610823457130565572459346412870218046009540557861328125000000000000000000000000000000e-30"),
  2047. cpp_dec_float("3.155443620884047221646914261131144918692825740436092019081115722656250000000000000000000000000000000e-30"),
  2048. cpp_dec_float("6.310887241768094443293828522262289837385651480872184038162231445312500000000000000000000000000000000e-30"),
  2049. cpp_dec_float("1.262177448353618888658765704452457967477130296174436807632446289062500000000000000000000000000000000e-29"),
  2050. cpp_dec_float("2.524354896707237777317531408904915934954260592348873615264892578125000000000000000000000000000000000e-29"),
  2051. cpp_dec_float("5.048709793414475554635062817809831869908521184697747230529785156250000000000000000000000000000000000e-29"),
  2052. cpp_dec_float("1.009741958682895110927012563561966373981704236939549446105957031250000000000000000000000000000000000e-28"),
  2053. cpp_dec_float("2.019483917365790221854025127123932747963408473879098892211914062500000000000000000000000000000000000e-28"),
  2054. cpp_dec_float("4.038967834731580443708050254247865495926816947758197784423828125000000000000000000000000000000000000e-28"),
  2055. cpp_dec_float("8.077935669463160887416100508495730991853633895516395568847656250000000000000000000000000000000000000e-28"),
  2056. cpp_dec_float("1.615587133892632177483220101699146198370726779103279113769531250000000000000000000000000000000000000e-27"),
  2057. cpp_dec_float("3.231174267785264354966440203398292396741453558206558227539062500000000000000000000000000000000000000e-27"),
  2058. cpp_dec_float("6.462348535570528709932880406796584793482907116413116455078125000000000000000000000000000000000000000e-27"),
  2059. cpp_dec_float("1.292469707114105741986576081359316958696581423282623291015625000000000000000000000000000000000000000e-26"),
  2060. cpp_dec_float("2.584939414228211483973152162718633917393162846565246582031250000000000000000000000000000000000000000e-26"),
  2061. cpp_dec_float("5.169878828456422967946304325437267834786325693130493164062500000000000000000000000000000000000000000e-26"),
  2062. cpp_dec_float("1.033975765691284593589260865087453566957265138626098632812500000000000000000000000000000000000000000e-25"),
  2063. cpp_dec_float("2.067951531382569187178521730174907133914530277252197265625000000000000000000000000000000000000000000e-25"),
  2064. cpp_dec_float("4.135903062765138374357043460349814267829060554504394531250000000000000000000000000000000000000000000e-25"),
  2065. cpp_dec_float("8.271806125530276748714086920699628535658121109008789062500000000000000000000000000000000000000000000e-25"),
  2066. cpp_dec_float("1.654361225106055349742817384139925707131624221801757812500000000000000000000000000000000000000000000e-24"),
  2067. cpp_dec_float("3.308722450212110699485634768279851414263248443603515625000000000000000000000000000000000000000000000e-24"),
  2068. cpp_dec_float("6.617444900424221398971269536559702828526496887207031250000000000000000000000000000000000000000000000e-24"),
  2069. cpp_dec_float("1.323488980084844279794253907311940565705299377441406250000000000000000000000000000000000000000000000e-23"),
  2070. cpp_dec_float("2.646977960169688559588507814623881131410598754882812500000000000000000000000000000000000000000000000e-23"),
  2071. cpp_dec_float("5.293955920339377119177015629247762262821197509765625000000000000000000000000000000000000000000000000e-23"),
  2072. cpp_dec_float("1.058791184067875423835403125849552452564239501953125000000000000000000000000000000000000000000000000e-22"),
  2073. cpp_dec_float("2.117582368135750847670806251699104905128479003906250000000000000000000000000000000000000000000000000e-22"),
  2074. cpp_dec_float("4.235164736271501695341612503398209810256958007812500000000000000000000000000000000000000000000000000e-22"),
  2075. cpp_dec_float("8.470329472543003390683225006796419620513916015625000000000000000000000000000000000000000000000000000e-22"),
  2076. cpp_dec_float("1.694065894508600678136645001359283924102783203125000000000000000000000000000000000000000000000000000e-21"),
  2077. cpp_dec_float("3.388131789017201356273290002718567848205566406250000000000000000000000000000000000000000000000000000e-21"),
  2078. cpp_dec_float("6.776263578034402712546580005437135696411132812500000000000000000000000000000000000000000000000000000e-21"),
  2079. cpp_dec_float("1.355252715606880542509316001087427139282226562500000000000000000000000000000000000000000000000000000e-20"),
  2080. cpp_dec_float("2.710505431213761085018632002174854278564453125000000000000000000000000000000000000000000000000000000e-20"),
  2081. cpp_dec_float("5.421010862427522170037264004349708557128906250000000000000000000000000000000000000000000000000000000e-20"),
  2082. cpp_dec_float("1.084202172485504434007452800869941711425781250000000000000000000000000000000000000000000000000000000e-19"),
  2083. cpp_dec_float("2.168404344971008868014905601739883422851562500000000000000000000000000000000000000000000000000000000e-19"),
  2084. cpp_dec_float("4.336808689942017736029811203479766845703125000000000000000000000000000000000000000000000000000000000e-19"),
  2085. cpp_dec_float("8.673617379884035472059622406959533691406250000000000000000000000000000000000000000000000000000000000e-19"),
  2086. cpp_dec_float("1.734723475976807094411924481391906738281250000000000000000000000000000000000000000000000000000000000e-18"),
  2087. cpp_dec_float("3.469446951953614188823848962783813476562500000000000000000000000000000000000000000000000000000000000e-18"),
  2088. cpp_dec_float("6.938893903907228377647697925567626953125000000000000000000000000000000000000000000000000000000000000e-18"),
  2089. cpp_dec_float("1.387778780781445675529539585113525390625000000000000000000000000000000000000000000000000000000000000e-17"),
  2090. cpp_dec_float("2.775557561562891351059079170227050781250000000000000000000000000000000000000000000000000000000000000e-17"),
  2091. cpp_dec_float("5.551115123125782702118158340454101562500000000000000000000000000000000000000000000000000000000000000e-17"),
  2092. cpp_dec_float("1.110223024625156540423631668090820312500000000000000000000000000000000000000000000000000000000000000e-16"),
  2093. cpp_dec_float("2.220446049250313080847263336181640625000000000000000000000000000000000000000000000000000000000000000e-16"),
  2094. cpp_dec_float("4.440892098500626161694526672363281250000000000000000000000000000000000000000000000000000000000000000e-16"),
  2095. cpp_dec_float("8.881784197001252323389053344726562500000000000000000000000000000000000000000000000000000000000000000e-16"),
  2096. cpp_dec_float("1.776356839400250464677810668945312500000000000000000000000000000000000000000000000000000000000000000e-15"),
  2097. cpp_dec_float("3.552713678800500929355621337890625000000000000000000000000000000000000000000000000000000000000000000e-15"),
  2098. cpp_dec_float("7.105427357601001858711242675781250000000000000000000000000000000000000000000000000000000000000000000e-15"),
  2099. cpp_dec_float("1.421085471520200371742248535156250000000000000000000000000000000000000000000000000000000000000000000e-14"),
  2100. cpp_dec_float("2.842170943040400743484497070312500000000000000000000000000000000000000000000000000000000000000000000e-14"),
  2101. cpp_dec_float("5.684341886080801486968994140625000000000000000000000000000000000000000000000000000000000000000000000e-14"),
  2102. cpp_dec_float("1.136868377216160297393798828125000000000000000000000000000000000000000000000000000000000000000000000e-13"),
  2103. cpp_dec_float("2.273736754432320594787597656250000000000000000000000000000000000000000000000000000000000000000000000e-13"),
  2104. cpp_dec_float("4.547473508864641189575195312500000000000000000000000000000000000000000000000000000000000000000000000e-13"),
  2105. cpp_dec_float("9.094947017729282379150390625000000000000000000000000000000000000000000000000000000000000000000000000e-13"),
  2106. cpp_dec_float("1.818989403545856475830078125000000000000000000000000000000000000000000000000000000000000000000000000e-12"),
  2107. cpp_dec_float("3.637978807091712951660156250000000000000000000000000000000000000000000000000000000000000000000000000e-12"),
  2108. cpp_dec_float("7.275957614183425903320312500000000000000000000000000000000000000000000000000000000000000000000000000e-12"),
  2109. cpp_dec_float("1.455191522836685180664062500000000000000000000000000000000000000000000000000000000000000000000000000e-11"),
  2110. cpp_dec_float("2.910383045673370361328125000000000000000000000000000000000000000000000000000000000000000000000000000e-11"),
  2111. cpp_dec_float("5.820766091346740722656250000000000000000000000000000000000000000000000000000000000000000000000000000e-11"),
  2112. cpp_dec_float("1.164153218269348144531250000000000000000000000000000000000000000000000000000000000000000000000000000e-10"),
  2113. cpp_dec_float("2.328306436538696289062500000000000000000000000000000000000000000000000000000000000000000000000000000e-10"),
  2114. cpp_dec_float("4.656612873077392578125000000000000000000000000000000000000000000000000000000000000000000000000000000e-10"),
  2115. cpp_dec_float("9.313225746154785156250000000000000000000000000000000000000000000000000000000000000000000000000000000e-10"),
  2116. cpp_dec_float("1.862645149230957031250000000000000000000000000000000000000000000000000000000000000000000000000000000e-9"),
  2117. cpp_dec_float("3.725290298461914062500000000000000000000000000000000000000000000000000000000000000000000000000000000e-9"),
  2118. cpp_dec_float("7.450580596923828125000000000000000000000000000000000000000000000000000000000000000000000000000000000e-9"),
  2119. cpp_dec_float("1.490116119384765625000000000000000000000000000000000000000000000000000000000000000000000000000000000e-8"),
  2120. cpp_dec_float("2.980232238769531250000000000000000000000000000000000000000000000000000000000000000000000000000000000e-8"),
  2121. cpp_dec_float("5.960464477539062500000000000000000000000000000000000000000000000000000000000000000000000000000000000e-8"),
  2122. cpp_dec_float("1.192092895507812500000000000000000000000000000000000000000000000000000000000000000000000000000000000e-7"),
  2123. cpp_dec_float("2.384185791015625000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-7"),
  2124. cpp_dec_float("4.768371582031250000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-7"),
  2125. cpp_dec_float("9.536743164062500000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-7"),
  2126. cpp_dec_float("1.907348632812500000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-6"),
  2127. cpp_dec_float("3.814697265625000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-6"),
  2128. cpp_dec_float("7.629394531250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-6"),
  2129. cpp_dec_float("0.000015258789062500000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2130. cpp_dec_float("0.000030517578125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2131. cpp_dec_float("0.000061035156250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2132. cpp_dec_float("0.000122070312500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2133. cpp_dec_float("0.000244140625000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2134. cpp_dec_float("0.000488281250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2135. cpp_dec_float("0.000976562500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2136. cpp_dec_float("0.001953125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2137. cpp_dec_float("0.003906250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2138. cpp_dec_float("0.007812500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2139. cpp_dec_float("0.01562500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2140. cpp_dec_float("0.03125000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2141. cpp_dec_float("0.06250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"),
  2142. cpp_dec_float("0.125"),
  2143. cpp_dec_float("0.25"),
  2144. cpp_dec_float("0.5"),
  2145. one(),
  2146. two(),
  2147. cpp_dec_float(static_cast<boost::ulong_long_type>(4)),
  2148. cpp_dec_float(static_cast<boost::ulong_long_type>(8)),
  2149. cpp_dec_float(static_cast<boost::ulong_long_type>(16)),
  2150. cpp_dec_float(static_cast<boost::ulong_long_type>(32)),
  2151. cpp_dec_float(static_cast<boost::ulong_long_type>(64)),
  2152. cpp_dec_float(static_cast<boost::ulong_long_type>(128)),
  2153. cpp_dec_float(static_cast<boost::ulong_long_type>(256)),
  2154. cpp_dec_float(static_cast<boost::ulong_long_type>(512)),
  2155. cpp_dec_float(static_cast<boost::ulong_long_type>(1024)),
  2156. cpp_dec_float(static_cast<boost::ulong_long_type>(2048)),
  2157. cpp_dec_float(static_cast<boost::ulong_long_type>(4096)),
  2158. cpp_dec_float(static_cast<boost::ulong_long_type>(8192)),
  2159. cpp_dec_float(static_cast<boost::ulong_long_type>(16384)),
  2160. cpp_dec_float(static_cast<boost::ulong_long_type>(32768)),
  2161. cpp_dec_float(static_cast<boost::ulong_long_type>(65536)),
  2162. cpp_dec_float(static_cast<boost::ulong_long_type>(131072)),
  2163. cpp_dec_float(static_cast<boost::ulong_long_type>(262144)),
  2164. cpp_dec_float(static_cast<boost::ulong_long_type>(524288)),
  2165. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 20u)),
  2166. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 21u)),
  2167. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 22u)),
  2168. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 23u)),
  2169. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 24u)),
  2170. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 25u)),
  2171. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 26u)),
  2172. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 27u)),
  2173. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 28u)),
  2174. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 29u)),
  2175. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 30u)),
  2176. cpp_dec_float(static_cast<boost::uint64_t>(1uL << 31u)),
  2177. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 32u)),
  2178. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 33u)),
  2179. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 34u)),
  2180. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 35u)),
  2181. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 36u)),
  2182. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 37u)),
  2183. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 38u)),
  2184. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 39u)),
  2185. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 40u)),
  2186. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 41u)),
  2187. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 42u)),
  2188. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 43u)),
  2189. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 44u)),
  2190. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 45u)),
  2191. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 46u)),
  2192. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 47u)),
  2193. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 48u)),
  2194. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 49u)),
  2195. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 50u)),
  2196. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 51u)),
  2197. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 52u)),
  2198. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 53u)),
  2199. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 54u)),
  2200. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 55u)),
  2201. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 56u)),
  2202. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 57u)),
  2203. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 58u)),
  2204. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 59u)),
  2205. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 60u)),
  2206. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 61u)),
  2207. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 62u)),
  2208. cpp_dec_float(static_cast<boost::uint64_t>(1uLL << 63u)),
  2209. cpp_dec_float("1.844674407370955161600000000000000000000000000000000000000000000000000000000000000000000000000000000e19"),
  2210. cpp_dec_float("3.689348814741910323200000000000000000000000000000000000000000000000000000000000000000000000000000000e19"),
  2211. cpp_dec_float("7.378697629483820646400000000000000000000000000000000000000000000000000000000000000000000000000000000e19"),
  2212. cpp_dec_float("1.475739525896764129280000000000000000000000000000000000000000000000000000000000000000000000000000000e20"),
  2213. cpp_dec_float("2.951479051793528258560000000000000000000000000000000000000000000000000000000000000000000000000000000e20"),
  2214. cpp_dec_float("5.902958103587056517120000000000000000000000000000000000000000000000000000000000000000000000000000000e20"),
  2215. cpp_dec_float("1.180591620717411303424000000000000000000000000000000000000000000000000000000000000000000000000000000e21"),
  2216. cpp_dec_float("2.361183241434822606848000000000000000000000000000000000000000000000000000000000000000000000000000000e21"),
  2217. cpp_dec_float("4.722366482869645213696000000000000000000000000000000000000000000000000000000000000000000000000000000e21"),
  2218. cpp_dec_float("9.444732965739290427392000000000000000000000000000000000000000000000000000000000000000000000000000000e21"),
  2219. cpp_dec_float("1.888946593147858085478400000000000000000000000000000000000000000000000000000000000000000000000000000e22"),
  2220. cpp_dec_float("3.777893186295716170956800000000000000000000000000000000000000000000000000000000000000000000000000000e22"),
  2221. cpp_dec_float("7.555786372591432341913600000000000000000000000000000000000000000000000000000000000000000000000000000e22"),
  2222. cpp_dec_float("1.511157274518286468382720000000000000000000000000000000000000000000000000000000000000000000000000000e23"),
  2223. cpp_dec_float("3.022314549036572936765440000000000000000000000000000000000000000000000000000000000000000000000000000e23"),
  2224. cpp_dec_float("6.044629098073145873530880000000000000000000000000000000000000000000000000000000000000000000000000000e23"),
  2225. cpp_dec_float("1.208925819614629174706176000000000000000000000000000000000000000000000000000000000000000000000000000e24"),
  2226. cpp_dec_float("2.417851639229258349412352000000000000000000000000000000000000000000000000000000000000000000000000000e24"),
  2227. cpp_dec_float("4.835703278458516698824704000000000000000000000000000000000000000000000000000000000000000000000000000e24"),
  2228. cpp_dec_float("9.671406556917033397649408000000000000000000000000000000000000000000000000000000000000000000000000000e24"),
  2229. cpp_dec_float("1.934281311383406679529881600000000000000000000000000000000000000000000000000000000000000000000000000e25"),
  2230. cpp_dec_float("3.868562622766813359059763200000000000000000000000000000000000000000000000000000000000000000000000000e25"),
  2231. cpp_dec_float("7.737125245533626718119526400000000000000000000000000000000000000000000000000000000000000000000000000e25"),
  2232. cpp_dec_float("1.547425049106725343623905280000000000000000000000000000000000000000000000000000000000000000000000000e26"),
  2233. cpp_dec_float("3.094850098213450687247810560000000000000000000000000000000000000000000000000000000000000000000000000e26"),
  2234. cpp_dec_float("6.189700196426901374495621120000000000000000000000000000000000000000000000000000000000000000000000000e26"),
  2235. cpp_dec_float("1.237940039285380274899124224000000000000000000000000000000000000000000000000000000000000000000000000e27"),
  2236. cpp_dec_float("2.475880078570760549798248448000000000000000000000000000000000000000000000000000000000000000000000000e27"),
  2237. cpp_dec_float("4.951760157141521099596496896000000000000000000000000000000000000000000000000000000000000000000000000e27"),
  2238. cpp_dec_float("9.903520314283042199192993792000000000000000000000000000000000000000000000000000000000000000000000000e27"),
  2239. cpp_dec_float("1.980704062856608439838598758400000000000000000000000000000000000000000000000000000000000000000000000e28"),
  2240. cpp_dec_float("3.961408125713216879677197516800000000000000000000000000000000000000000000000000000000000000000000000e28"),
  2241. cpp_dec_float("7.922816251426433759354395033600000000000000000000000000000000000000000000000000000000000000000000000e28"),
  2242. cpp_dec_float("1.584563250285286751870879006720000000000000000000000000000000000000000000000000000000000000000000000e29"),
  2243. cpp_dec_float("3.169126500570573503741758013440000000000000000000000000000000000000000000000000000000000000000000000e29"),
  2244. cpp_dec_float("6.338253001141147007483516026880000000000000000000000000000000000000000000000000000000000000000000000e29"),
  2245. cpp_dec_float("1.267650600228229401496703205376000000000000000000000000000000000000000000000000000000000000000000000e30"),
  2246. cpp_dec_float("2.535301200456458802993406410752000000000000000000000000000000000000000000000000000000000000000000000e30"),
  2247. cpp_dec_float("5.070602400912917605986812821504000000000000000000000000000000000000000000000000000000000000000000000e30"),
  2248. cpp_dec_float("1.014120480182583521197362564300800000000000000000000000000000000000000000000000000000000000000000000e31"),
  2249. cpp_dec_float("2.028240960365167042394725128601600000000000000000000000000000000000000000000000000000000000000000000e31"),
  2250. cpp_dec_float("4.056481920730334084789450257203200000000000000000000000000000000000000000000000000000000000000000000e31"),
  2251. cpp_dec_float("8.112963841460668169578900514406400000000000000000000000000000000000000000000000000000000000000000000e31"),
  2252. cpp_dec_float("1.622592768292133633915780102881280000000000000000000000000000000000000000000000000000000000000000000e32"),
  2253. cpp_dec_float("3.245185536584267267831560205762560000000000000000000000000000000000000000000000000000000000000000000e32"),
  2254. cpp_dec_float("6.490371073168534535663120411525120000000000000000000000000000000000000000000000000000000000000000000e32"),
  2255. cpp_dec_float("1.298074214633706907132624082305024000000000000000000000000000000000000000000000000000000000000000000e33"),
  2256. cpp_dec_float("2.596148429267413814265248164610048000000000000000000000000000000000000000000000000000000000000000000e33"),
  2257. cpp_dec_float("5.192296858534827628530496329220096000000000000000000000000000000000000000000000000000000000000000000e33"),
  2258. cpp_dec_float("1.038459371706965525706099265844019200000000000000000000000000000000000000000000000000000000000000000e34"),
  2259. cpp_dec_float("2.076918743413931051412198531688038400000000000000000000000000000000000000000000000000000000000000000e34"),
  2260. cpp_dec_float("4.153837486827862102824397063376076800000000000000000000000000000000000000000000000000000000000000000e34"),
  2261. cpp_dec_float("8.307674973655724205648794126752153600000000000000000000000000000000000000000000000000000000000000000e34"),
  2262. cpp_dec_float("1.661534994731144841129758825350430720000000000000000000000000000000000000000000000000000000000000000e35"),
  2263. cpp_dec_float("3.323069989462289682259517650700861440000000000000000000000000000000000000000000000000000000000000000e35"),
  2264. cpp_dec_float("6.646139978924579364519035301401722880000000000000000000000000000000000000000000000000000000000000000e35"),
  2265. cpp_dec_float("1.329227995784915872903807060280344576000000000000000000000000000000000000000000000000000000000000000e36"),
  2266. cpp_dec_float("2.658455991569831745807614120560689152000000000000000000000000000000000000000000000000000000000000000e36"),
  2267. cpp_dec_float("5.316911983139663491615228241121378304000000000000000000000000000000000000000000000000000000000000000e36"),
  2268. cpp_dec_float("1.063382396627932698323045648224275660800000000000000000000000000000000000000000000000000000000000000e37"),
  2269. cpp_dec_float("2.126764793255865396646091296448551321600000000000000000000000000000000000000000000000000000000000000e37"),
  2270. cpp_dec_float("4.253529586511730793292182592897102643200000000000000000000000000000000000000000000000000000000000000e37"),
  2271. cpp_dec_float("8.507059173023461586584365185794205286400000000000000000000000000000000000000000000000000000000000000e37"),
  2272. cpp_dec_float("1.701411834604692317316873037158841057280000000000000000000000000000000000000000000000000000000000000e38")
  2273. }};
  2274. if((p > static_cast<boost::long_long_type>(-128)) && (p < static_cast<boost::long_long_type>(+128)))
  2275. {
  2276. return p2_data[static_cast<std::size_t>(p + ((p2_data.size() - 1u) / 2u))];
  2277. }
  2278. else
  2279. {
  2280. // Compute and return 2^p.
  2281. if(p < static_cast<boost::long_long_type>(0))
  2282. {
  2283. return pow2(static_cast<boost::long_long_type>(-p)).calculate_inv();
  2284. }
  2285. else
  2286. {
  2287. cpp_dec_float<Digits10, ExponentType, Allocator> t;
  2288. default_ops::detail::pow_imp(t, two(), p, mpl::true_());
  2289. return t;
  2290. }
  2291. }
  2292. }
  2293. template <unsigned Digits10, class ExponentType, class Allocator>
  2294. inline void eval_add(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& o)
  2295. {
  2296. result += o;
  2297. }
  2298. template <unsigned Digits10, class ExponentType, class Allocator>
  2299. inline void eval_subtract(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& o)
  2300. {
  2301. result -= o;
  2302. }
  2303. template <unsigned Digits10, class ExponentType, class Allocator>
  2304. inline void eval_multiply(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& o)
  2305. {
  2306. result *= o;
  2307. }
  2308. template <unsigned Digits10, class ExponentType, class Allocator>
  2309. inline void eval_divide(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& o)
  2310. {
  2311. result /= o;
  2312. }
  2313. template <unsigned Digits10, class ExponentType, class Allocator>
  2314. inline void eval_add(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const boost::ulong_long_type& o)
  2315. {
  2316. result.add_unsigned_long_long(o);
  2317. }
  2318. template <unsigned Digits10, class ExponentType, class Allocator>
  2319. inline void eval_subtract(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const boost::ulong_long_type& o)
  2320. {
  2321. result.sub_unsigned_long_long(o);
  2322. }
  2323. template <unsigned Digits10, class ExponentType, class Allocator>
  2324. inline void eval_multiply(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const boost::ulong_long_type& o)
  2325. {
  2326. result.mul_unsigned_long_long(o);
  2327. }
  2328. template <unsigned Digits10, class ExponentType, class Allocator>
  2329. inline void eval_divide(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const boost::ulong_long_type& o)
  2330. {
  2331. result.div_unsigned_long_long(o);
  2332. }
  2333. template <unsigned Digits10, class ExponentType, class Allocator>
  2334. inline void eval_add(cpp_dec_float<Digits10, ExponentType, Allocator>& result, boost::long_long_type o)
  2335. {
  2336. if(o < 0)
  2337. result.sub_unsigned_long_long(boost::multiprecision::detail::unsigned_abs(o));
  2338. else
  2339. result.add_unsigned_long_long(o);
  2340. }
  2341. template <unsigned Digits10, class ExponentType, class Allocator>
  2342. inline void eval_subtract(cpp_dec_float<Digits10, ExponentType, Allocator>& result, boost::long_long_type o)
  2343. {
  2344. if(o < 0)
  2345. result.add_unsigned_long_long(boost::multiprecision::detail::unsigned_abs(o));
  2346. else
  2347. result.sub_unsigned_long_long(o);
  2348. }
  2349. template <unsigned Digits10, class ExponentType, class Allocator>
  2350. inline void eval_multiply(cpp_dec_float<Digits10, ExponentType, Allocator>& result, boost::long_long_type o)
  2351. {
  2352. if(o < 0)
  2353. {
  2354. result.mul_unsigned_long_long(boost::multiprecision::detail::unsigned_abs(o));
  2355. result.negate();
  2356. }
  2357. else
  2358. result.mul_unsigned_long_long(o);
  2359. }
  2360. template <unsigned Digits10, class ExponentType, class Allocator>
  2361. inline void eval_divide(cpp_dec_float<Digits10, ExponentType, Allocator>& result, boost::long_long_type o)
  2362. {
  2363. if(o < 0)
  2364. {
  2365. result.div_unsigned_long_long(boost::multiprecision::detail::unsigned_abs(o));
  2366. result.negate();
  2367. }
  2368. else
  2369. result.div_unsigned_long_long(o);
  2370. }
  2371. template <unsigned Digits10, class ExponentType, class Allocator>
  2372. inline void eval_convert_to(boost::ulong_long_type* result, const cpp_dec_float<Digits10, ExponentType, Allocator>& val)
  2373. {
  2374. *result = val.extract_unsigned_long_long();
  2375. }
  2376. template <unsigned Digits10, class ExponentType, class Allocator>
  2377. inline void eval_convert_to(boost::long_long_type* result, const cpp_dec_float<Digits10, ExponentType, Allocator>& val)
  2378. {
  2379. *result = val.extract_signed_long_long();
  2380. }
  2381. template <unsigned Digits10, class ExponentType, class Allocator>
  2382. inline void eval_convert_to(long double* result, cpp_dec_float<Digits10, ExponentType, Allocator>& val)
  2383. {
  2384. *result = val.extract_long_double();
  2385. }
  2386. //
  2387. // Non member function support:
  2388. //
  2389. template <unsigned Digits10, class ExponentType, class Allocator>
  2390. inline int eval_fpclassify(const cpp_dec_float<Digits10, ExponentType, Allocator>& x)
  2391. {
  2392. if((x.isinf)())
  2393. return FP_INFINITE;
  2394. if((x.isnan)())
  2395. return FP_NAN;
  2396. if(x.iszero())
  2397. return FP_ZERO;
  2398. return FP_NORMAL;
  2399. }
  2400. template <unsigned Digits10, class ExponentType, class Allocator>
  2401. inline void eval_abs(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x)
  2402. {
  2403. result = x;
  2404. if(x.isneg())
  2405. result.negate();
  2406. }
  2407. template <unsigned Digits10, class ExponentType, class Allocator>
  2408. inline void eval_fabs(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x)
  2409. {
  2410. result = x;
  2411. if(x.isneg())
  2412. result.negate();
  2413. }
  2414. template <unsigned Digits10, class ExponentType, class Allocator>
  2415. inline void eval_sqrt(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x)
  2416. {
  2417. result = x;
  2418. result.calculate_sqrt();
  2419. }
  2420. template <unsigned Digits10, class ExponentType, class Allocator>
  2421. inline void eval_floor(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x)
  2422. {
  2423. result = x;
  2424. if(!(x.isfinite)() || x.isint())
  2425. {
  2426. if((x.isnan)())
  2427. errno = EDOM;
  2428. return;
  2429. }
  2430. if(x.isneg())
  2431. result -= cpp_dec_float<Digits10, ExponentType, Allocator>::one();
  2432. result = result.extract_integer_part();
  2433. }
  2434. template <unsigned Digits10, class ExponentType, class Allocator>
  2435. inline void eval_ceil(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x)
  2436. {
  2437. result = x;
  2438. if(!(x.isfinite)() || x.isint())
  2439. {
  2440. if((x.isnan)())
  2441. errno = EDOM;
  2442. return;
  2443. }
  2444. if(!x.isneg())
  2445. result += cpp_dec_float<Digits10, ExponentType, Allocator>::one();
  2446. result = result.extract_integer_part();
  2447. }
  2448. template <unsigned Digits10, class ExponentType, class Allocator>
  2449. inline void eval_trunc(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x)
  2450. {
  2451. if(x.isint() || !(x.isfinite)())
  2452. {
  2453. result = x;
  2454. if((x.isnan)())
  2455. errno = EDOM;
  2456. return;
  2457. }
  2458. result = x.extract_integer_part();
  2459. }
  2460. template <unsigned Digits10, class ExponentType, class Allocator>
  2461. inline ExponentType eval_ilogb(const cpp_dec_float<Digits10, ExponentType, Allocator>& val)
  2462. {
  2463. if(val.iszero())
  2464. return (std::numeric_limits<ExponentType>::min)();
  2465. if((val.isinf)())
  2466. return INT_MAX;
  2467. if((val.isnan)())
  2468. #ifdef FP_ILOGBNAN
  2469. return FP_ILOGBNAN;
  2470. #else
  2471. return INT_MAX;
  2472. #endif
  2473. // Set result, to the exponent of val:
  2474. return val.order();
  2475. }
  2476. template <unsigned Digits10, class ExponentType, class Allocator, class ArgType>
  2477. inline void eval_scalbn(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& val, ArgType e_)
  2478. {
  2479. using default_ops::eval_multiply;
  2480. const ExponentType e = static_cast<ExponentType>(e_);
  2481. cpp_dec_float<Digits10, ExponentType, Allocator> t(1.0, e);
  2482. eval_multiply(result, val, t);
  2483. }
  2484. template <unsigned Digits10, class ExponentType, class Allocator, class ArgType>
  2485. inline void eval_ldexp(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x, ArgType e)
  2486. {
  2487. const boost::long_long_type the_exp = static_cast<boost::long_long_type>(e);
  2488. if((the_exp > (std::numeric_limits<ExponentType>::max)()) || (the_exp < (std::numeric_limits<ExponentType>::min)()))
  2489. BOOST_THROW_EXCEPTION(std::runtime_error(std::string("Exponent value is out of range.")));
  2490. result = x;
  2491. if ((the_exp > static_cast<boost::long_long_type>(-std::numeric_limits<boost::long_long_type>::digits)) && (the_exp < static_cast<boost::long_long_type>(0)))
  2492. result.div_unsigned_long_long(1ULL << static_cast<boost::long_long_type>(-the_exp));
  2493. else if((the_exp < static_cast<boost::long_long_type>( std::numeric_limits<boost::long_long_type>::digits)) && (the_exp > static_cast<boost::long_long_type>(0)))
  2494. result.mul_unsigned_long_long(1ULL << the_exp);
  2495. else if(the_exp != static_cast<boost::long_long_type>(0))
  2496. result *= cpp_dec_float<Digits10, ExponentType, Allocator>::pow2(e);
  2497. }
  2498. template <unsigned Digits10, class ExponentType, class Allocator>
  2499. inline void eval_frexp(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x, ExponentType* e)
  2500. {
  2501. result = x;
  2502. if(result.iszero() || (result.isinf)() || (result.isnan)())
  2503. {
  2504. *e = 0;
  2505. return;
  2506. }
  2507. if(result.isneg())
  2508. result.negate();
  2509. ExponentType t = result.order();
  2510. BOOST_MP_USING_ABS
  2511. if(abs(t) < ((std::numeric_limits<ExponentType>::max)() / 1000))
  2512. {
  2513. t *= 1000;
  2514. t /= 301;
  2515. }
  2516. else
  2517. {
  2518. t /= 301;
  2519. t *= 1000;
  2520. }
  2521. result *= cpp_dec_float<Digits10, ExponentType, Allocator>::pow2(-t);
  2522. if(result.iszero() || (result.isinf)() || (result.isnan)())
  2523. {
  2524. // pow2 overflowed, slip the calculation up:
  2525. result = x;
  2526. if(result.isneg())
  2527. result.negate();
  2528. t /= 2;
  2529. result *= cpp_dec_float<Digits10, ExponentType, Allocator>::pow2(-t);
  2530. }
  2531. BOOST_MP_USING_ABS
  2532. if(abs(result.order()) > 5)
  2533. {
  2534. // If our first estimate doesn't get close enough then try recursion until we do:
  2535. ExponentType e2;
  2536. cpp_dec_float<Digits10, ExponentType, Allocator> r2;
  2537. eval_frexp(r2, result, &e2);
  2538. // overflow protection:
  2539. if((t > 0) && (e2 > 0) && (t > (std::numeric_limits<ExponentType>::max)() - e2))
  2540. BOOST_THROW_EXCEPTION(std::runtime_error("Exponent is too large to be represented as a power of 2."));
  2541. if((t < 0) && (e2 < 0) && (t < (std::numeric_limits<ExponentType>::min)() - e2))
  2542. BOOST_THROW_EXCEPTION(std::runtime_error("Exponent is too large to be represented as a power of 2."));
  2543. t += e2;
  2544. result = r2;
  2545. }
  2546. while(result.compare(cpp_dec_float<Digits10, ExponentType, Allocator>::one()) >= 0)
  2547. {
  2548. result /= cpp_dec_float<Digits10, ExponentType, Allocator>::two();
  2549. ++t;
  2550. }
  2551. while(result.compare(cpp_dec_float<Digits10, ExponentType, Allocator>::half()) < 0)
  2552. {
  2553. result *= cpp_dec_float<Digits10, ExponentType, Allocator>::two();
  2554. --t;
  2555. }
  2556. *e = t;
  2557. if(x.isneg())
  2558. result.negate();
  2559. }
  2560. template <unsigned Digits10, class ExponentType, class Allocator>
  2561. inline typename disable_if<is_same<ExponentType, int> >::type eval_frexp(cpp_dec_float<Digits10, ExponentType, Allocator>& result, const cpp_dec_float<Digits10, ExponentType, Allocator>& x, int* e)
  2562. {
  2563. ExponentType t;
  2564. eval_frexp(result, x, &t);
  2565. if((t > (std::numeric_limits<int>::max)()) || (t < (std::numeric_limits<int>::min)()))
  2566. BOOST_THROW_EXCEPTION(std::runtime_error("Exponent is outside the range of an int"));
  2567. *e = static_cast<int>(t);
  2568. }
  2569. template <unsigned Digits10, class ExponentType, class Allocator>
  2570. inline bool eval_is_zero(const cpp_dec_float<Digits10, ExponentType, Allocator>& val)
  2571. {
  2572. return val.iszero();
  2573. }
  2574. template <unsigned Digits10, class ExponentType, class Allocator>
  2575. inline int eval_get_sign(const cpp_dec_float<Digits10, ExponentType, Allocator>& val)
  2576. {
  2577. return val.iszero() ? 0 : val.isneg() ? -1 : 1;
  2578. }
  2579. template <unsigned Digits10, class ExponentType, class Allocator>
  2580. inline std::size_t hash_value(const cpp_dec_float<Digits10, ExponentType, Allocator>& val)
  2581. {
  2582. return val.hash();
  2583. }
  2584. } // namespace backends
  2585. using boost::multiprecision::backends::cpp_dec_float;
  2586. typedef number<cpp_dec_float<50> > cpp_dec_float_50;
  2587. typedef number<cpp_dec_float<100> > cpp_dec_float_100;
  2588. #ifdef BOOST_NO_SFINAE_EXPR
  2589. namespace detail{
  2590. template<unsigned D1, class E1, class A1, unsigned D2, class E2, class A2>
  2591. struct is_explicitly_convertible<cpp_dec_float<D1, E1, A1>, cpp_dec_float<D2, E2, A2> > : public mpl::true_ {};
  2592. }
  2593. #endif
  2594. }}
  2595. namespace std
  2596. {
  2597. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2598. class numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >
  2599. {
  2600. public:
  2601. BOOST_STATIC_CONSTEXPR bool is_specialized = true;
  2602. BOOST_STATIC_CONSTEXPR bool is_signed = true;
  2603. BOOST_STATIC_CONSTEXPR bool is_integer = false;
  2604. BOOST_STATIC_CONSTEXPR bool is_exact = false;
  2605. BOOST_STATIC_CONSTEXPR bool is_bounded = true;
  2606. BOOST_STATIC_CONSTEXPR bool is_modulo = false;
  2607. BOOST_STATIC_CONSTEXPR bool is_iec559 = false;
  2608. BOOST_STATIC_CONSTEXPR int digits = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10;
  2609. BOOST_STATIC_CONSTEXPR int digits10 = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10;
  2610. BOOST_STATIC_CONSTEXPR int max_digits10 = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_total_digits10;
  2611. BOOST_STATIC_CONSTEXPR ExponentType min_exponent = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp; // Type differs from int.
  2612. BOOST_STATIC_CONSTEXPR ExponentType min_exponent10 = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_min_exp10; // Type differs from int.
  2613. BOOST_STATIC_CONSTEXPR ExponentType max_exponent = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp; // Type differs from int.
  2614. BOOST_STATIC_CONSTEXPR ExponentType max_exponent10 = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_max_exp10; // Type differs from int.
  2615. BOOST_STATIC_CONSTEXPR int radix = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_radix;
  2616. BOOST_STATIC_CONSTEXPR std::float_round_style round_style = std::round_indeterminate;
  2617. BOOST_STATIC_CONSTEXPR bool has_infinity = true;
  2618. BOOST_STATIC_CONSTEXPR bool has_quiet_NaN = true;
  2619. BOOST_STATIC_CONSTEXPR bool has_signaling_NaN = false;
  2620. BOOST_STATIC_CONSTEXPR std::float_denorm_style has_denorm = std::denorm_absent;
  2621. BOOST_STATIC_CONSTEXPR bool has_denorm_loss = false;
  2622. BOOST_STATIC_CONSTEXPR bool traps = false;
  2623. BOOST_STATIC_CONSTEXPR bool tinyness_before = false;
  2624. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> (min) () { return (boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::min)(); }
  2625. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> (max) () { return (boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::max)(); }
  2626. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> lowest () { return boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::zero(); }
  2627. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> epsilon () { return boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::eps(); }
  2628. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> round_error () { return 0.5L; }
  2629. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> infinity () { return boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::inf(); }
  2630. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> quiet_NaN () { return boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::nan(); }
  2631. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> signaling_NaN() { return boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::zero(); }
  2632. BOOST_STATIC_CONSTEXPR boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> denorm_min () { return boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::zero(); }
  2633. };
  2634. #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
  2635. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2636. BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::digits;
  2637. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2638. BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::digits10;
  2639. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2640. BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::max_digits10;
  2641. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2642. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::is_signed;
  2643. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2644. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::is_integer;
  2645. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2646. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::is_exact;
  2647. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2648. BOOST_CONSTEXPR_OR_CONST int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::radix;
  2649. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2650. BOOST_CONSTEXPR_OR_CONST ExponentType numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::min_exponent;
  2651. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2652. BOOST_CONSTEXPR_OR_CONST ExponentType numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::min_exponent10;
  2653. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2654. BOOST_CONSTEXPR_OR_CONST ExponentType numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::max_exponent;
  2655. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2656. BOOST_CONSTEXPR_OR_CONST ExponentType numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::max_exponent10;
  2657. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2658. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::has_infinity;
  2659. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2660. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::has_quiet_NaN;
  2661. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2662. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::has_signaling_NaN;
  2663. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2664. BOOST_CONSTEXPR_OR_CONST float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::has_denorm;
  2665. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2666. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::has_denorm_loss;
  2667. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2668. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::is_iec559;
  2669. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2670. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::is_bounded;
  2671. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2672. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::is_modulo;
  2673. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2674. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::traps;
  2675. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2676. BOOST_CONSTEXPR_OR_CONST bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::tinyness_before;
  2677. template <unsigned Digits10, class ExponentType, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  2678. BOOST_CONSTEXPR_OR_CONST float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates> >::round_style;
  2679. #endif
  2680. }
  2681. namespace boost{ namespace math{
  2682. namespace policies{
  2683. template <unsigned Digits10, class ExponentType, class Allocator, class Policy, boost::multiprecision::expression_template_option ExpressionTemplates>
  2684. struct precision< boost::multiprecision::number<boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>, ExpressionTemplates>, Policy>
  2685. {
  2686. // Define a local copy of cpp_dec_float_digits10 because it might differ
  2687. // from the template parameter Digits10 for small or large digit counts.
  2688. static const boost::int32_t cpp_dec_float_digits10 = boost::multiprecision::cpp_dec_float<Digits10, ExponentType, Allocator>::cpp_dec_float_digits10;
  2689. typedef typename Policy::precision_type precision_type;
  2690. typedef digits2<((cpp_dec_float_digits10 + 1LL) * 1000LL) / 301LL> digits_2;
  2691. typedef typename mpl::if_c<
  2692. ((digits_2::value <= precision_type::value)
  2693. || (Policy::precision_type::value <= 0)),
  2694. // Default case, full precision for RealType:
  2695. digits_2,
  2696. // User customized precision:
  2697. precision_type
  2698. >::type type;
  2699. };
  2700. } // namespace policies
  2701. }} // namespaces boost::math
  2702. #ifdef BOOST_MSVC
  2703. #pragma warning(pop)
  2704. #endif
  2705. #endif