cpp_bin_float.hpp 106 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365
  1. ////////////////////////////////////////////////////////////////
  2. // Copyright 2013 - 2025 John Maddock.
  3. // Copyright 2022 - 2025 Christopher Kormanyos.
  4. // Distributed under the Boost Software License,
  5. // Version 1.0. (See accompanying file LICENSE_1_0.txt
  6. // or copy at https://www.boost.org/LICENSE_1_0.txt)
  7. #ifndef BOOST_MP_CPP_BIN_FLOAT_HPP
  8. #define BOOST_MP_CPP_BIN_FLOAT_HPP
  9. #include <cmath>
  10. #include <cstdint>
  11. #include <limits>
  12. #include <type_traits>
  13. #include <boost/multiprecision/cpp_int.hpp>
  14. #include <boost/multiprecision/integer.hpp>
  15. #include <boost/multiprecision/detail/standalone_config.hpp>
  16. #include <boost/multiprecision/detail/fpclassify.hpp>
  17. #include <boost/multiprecision/detail/float_string_cvt.hpp>
  18. #include <boost/multiprecision/traits/max_digits10.hpp>
  19. #include <boost/multiprecision/detail/hash.hpp>
  20. #include <boost/multiprecision/detail/no_exceptions_support.hpp>
  21. #include <boost/multiprecision/detail/assert.hpp>
  22. #include <boost/multiprecision/detail/float128_functions.hpp>
  23. #include <boost/multiprecision/detail/functions/trunc.hpp>
  24. //
  25. // Some includes we need from Boost.Math, since we rely on that library to provide these functions:
  26. //
  27. #ifdef BOOST_MP_MATH_AVAILABLE
  28. #include <boost/math/special_functions/asinh.hpp>
  29. #include <boost/math/special_functions/acosh.hpp>
  30. #include <boost/math/special_functions/atanh.hpp>
  31. #include <boost/math/special_functions/cbrt.hpp>
  32. #include <boost/math/special_functions/expm1.hpp>
  33. #include <boost/math/special_functions/gamma.hpp>
  34. #endif
  35. #ifdef BOOST_HAS_FLOAT128
  36. # if __has_include(<quadmath.h>)
  37. # include <quadmath.h>
  38. # define BOOST_MP_HAS_FLOAT128_SUPPORT
  39. # endif
  40. #endif
  41. namespace boost {
  42. namespace multiprecision {
  43. namespace backends {
  44. #ifdef BOOST_MSVC
  45. #pragma warning(push)
  46. #pragma warning(disable : 4522 6326) // multiple assignment operators specified, comparison of two constants
  47. #endif
  48. namespace detail {
  49. template <class U>
  50. inline typename std::enable_if<boost::multiprecision::detail::is_unsigned<U>::value, bool>::type is_negative(U) { return false; }
  51. template <class S>
  52. inline typename std::enable_if< !boost::multiprecision::detail::is_unsigned<S>::value, bool>::type is_negative(S s) { return s < 0; }
  53. template <class Float, std::ptrdiff_t, bool = number_category<Float>::value == number_kind_floating_point>
  54. struct is_cpp_bin_float_implicitly_constructible_from_type
  55. {
  56. static constexpr bool value = false;
  57. };
  58. template <class Float, std::ptrdiff_t bit_count>
  59. struct is_cpp_bin_float_implicitly_constructible_from_type<Float, bit_count, true>
  60. {
  61. static constexpr bool value = (std::numeric_limits<Float>::digits <= static_cast<int>(bit_count)) && (std::numeric_limits<Float>::radix == 2) && std::numeric_limits<Float>::is_specialized
  62. #ifdef BOOST_MP_HAS_FLOAT128_SUPPORT
  63. && !std::is_same<Float, float128_type>::value
  64. #endif
  65. && (std::is_floating_point<Float>::value || is_number<Float>::value);
  66. };
  67. template <class Float, std::ptrdiff_t, bool = number_category<Float>::value == number_kind_floating_point>
  68. struct is_cpp_bin_float_explicitly_constructible_from_type
  69. {
  70. static constexpr bool value = false;
  71. };
  72. template <class Float, std::ptrdiff_t bit_count>
  73. struct is_cpp_bin_float_explicitly_constructible_from_type<Float, bit_count, true>
  74. {
  75. static constexpr bool value = (std::numeric_limits<Float>::digits > static_cast<int>(bit_count)) && (std::numeric_limits<Float>::radix == 2) && std::numeric_limits<Float>::is_specialized
  76. #ifdef BOOST_MP_HAS_FLOAT128_SUPPORT
  77. && !std::is_same<Float, float128_type>::value
  78. #endif
  79. ;
  80. };
  81. } // namespace detail
  82. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinExponent, Exponent MaxExponent>
  83. class cpp_bin_float
  84. {
  85. public:
  86. static constexpr unsigned bit_count = DigitBase == digit_base_2 ? Digits : (Digits * 1000uL) / 301uL + (((Digits * 1000uL) % 301) ? 2u : 1u);
  87. using rep_type = cpp_int_backend<std::is_void<Allocator>::value ? bit_count : 0, bit_count, std::is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator>;
  88. using double_rep_type = cpp_int_backend<std::is_void<Allocator>::value ? 2 * bit_count : 0, 2 * bit_count, std::is_void<Allocator>::value ? unsigned_magnitude : signed_magnitude, unchecked, Allocator>;
  89. using signed_types = typename rep_type::signed_types;
  90. using unsigned_types = typename rep_type::unsigned_types;
  91. using float_types = std::tuple<float, double, long double>;
  92. using exponent_type = Exponent;
  93. static constexpr exponent_type max_exponent_limit = (std::numeric_limits<exponent_type>::max)()- 2 * static_cast<exponent_type>(bit_count);
  94. static constexpr exponent_type min_exponent_limit = (std::numeric_limits<exponent_type>::min)() + 2 * static_cast<exponent_type>(bit_count);
  95. static_assert(MinExponent >= min_exponent_limit, "Template parameter MinExponent is too negative for our internal logic to function correctly, sorry!");
  96. static_assert(MaxExponent <= max_exponent_limit, "Template parameter MaxExponent is too large for our internal logic to function correctly, sorry!");
  97. static_assert(MinExponent <= 0, "Template parameter MinExponent can not be positive!");
  98. static_assert(MaxExponent >= 0, "Template parameter MaxExponent can not be negative!");
  99. static constexpr exponent_type max_exponent = MaxExponent == 0 ? max_exponent_limit : MaxExponent;
  100. static constexpr exponent_type min_exponent = MinExponent == 0 ? min_exponent_limit : MinExponent;
  101. static constexpr exponent_type exponent_zero = max_exponent + 1;
  102. static constexpr exponent_type exponent_infinity = max_exponent + 2;
  103. static constexpr exponent_type exponent_nan = max_exponent + 3;
  104. private:
  105. rep_type m_data;
  106. exponent_type m_exponent;
  107. bool m_sign;
  108. public:
  109. cpp_bin_float() noexcept(noexcept(rep_type())) : m_data(), m_exponent(exponent_zero), m_sign(false) {}
  110. cpp_bin_float(const cpp_bin_float& o) noexcept(noexcept(rep_type(std::declval<const rep_type&>())))
  111. : m_data(o.m_data), m_exponent(o.m_exponent), m_sign(o.m_sign) {}
  112. template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
  113. cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE>& o, typename std::enable_if<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = nullptr)
  114. {
  115. *this = o;
  116. }
  117. template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
  118. explicit cpp_bin_float(const cpp_bin_float<D, B, A, E, MinE, MaxE>& o, typename std::enable_if< !(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = nullptr)
  119. : m_exponent(o.exponent()), m_sign(o.sign())
  120. {
  121. *this = o;
  122. }
  123. // rvalue copy:
  124. template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
  125. cpp_bin_float(cpp_bin_float<D, B, A, E, MinE, MaxE>&& o, typename std::enable_if<(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = nullptr)noexcept(noexcept(rep_type(std::declval<rep_type&&>())))
  126. {
  127. *this = std::move(o);
  128. }
  129. template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
  130. explicit cpp_bin_float(cpp_bin_float<D, B, A, E, MinE, MaxE>&& o, typename std::enable_if< !(bit_count >= cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count)>::type const* = nullptr) noexcept(noexcept(rep_type(std::declval<rep_type&&>())))
  131. : m_exponent(o.exponent()), m_sign(o.sign())
  132. {
  133. *this = std::move(o);
  134. }
  135. template <class Float>
  136. cpp_bin_float(const Float& f,
  137. typename std::enable_if<detail::is_cpp_bin_float_implicitly_constructible_from_type<Float, static_cast<std::ptrdiff_t>(bit_count)>::value>::type const* = nullptr)
  138. : m_data(), m_exponent(0), m_sign(false)
  139. {
  140. this->assign_float(f);
  141. }
  142. template <class Float>
  143. explicit cpp_bin_float(const Float& f,
  144. typename std::enable_if<detail::is_cpp_bin_float_explicitly_constructible_from_type<Float, static_cast<std::ptrdiff_t>(bit_count)>::value>::type const* = nullptr)
  145. : m_data(), m_exponent(0), m_sign(false)
  146. {
  147. this->assign_float(f);
  148. }
  149. #ifdef BOOST_MP_HAS_FLOAT128_SUPPORT
  150. template <class Float>
  151. cpp_bin_float(const Float& f,
  152. typename std::enable_if<
  153. std::is_same<Float, float128_type>::value && (static_cast<int>(bit_count) >= 113)>::type const* = nullptr)
  154. : m_data(), m_exponent(0), m_sign(false)
  155. {
  156. this->assign_float(f);
  157. }
  158. template <class Float>
  159. explicit cpp_bin_float(const Float& f,
  160. typename std::enable_if<
  161. std::is_same<Float, float128_type>::value && (static_cast<int>(bit_count) < 113)>::type const* = nullptr)
  162. : m_data(), m_exponent(0), m_sign(false)
  163. {
  164. this->assign_float(f);
  165. }
  166. #endif
  167. cpp_bin_float& operator=(const cpp_bin_float& o) noexcept(noexcept(std::declval<rep_type&>() = std::declval<const rep_type&>()))
  168. {
  169. m_data = o.m_data;
  170. m_exponent = o.m_exponent;
  171. m_sign = o.m_sign;
  172. return *this;
  173. }
  174. template <class A, class E, E MinE, E MaxE>
  175. cpp_bin_float& operator=(const cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>& o) noexcept(noexcept(std::declval<rep_type&>() = std::declval<const rep_type&>()))
  176. {
  177. m_data = o.bits();
  178. m_sign = o.sign();
  179. if (o.exponent() == cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>::exponent_zero)
  180. m_exponent = exponent_zero;
  181. else if (o.exponent() == cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>::exponent_nan)
  182. m_exponent = exponent_nan;
  183. else if (o.exponent() == cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>::exponent_infinity)
  184. m_exponent = exponent_infinity;
  185. else if (o.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
  186. {
  187. // Overflow:
  188. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
  189. bits() = static_cast<limb_type>(0u);
  190. }
  191. else if (o.exponent() < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
  192. {
  193. // Underflow:
  194. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  195. bits() = static_cast<limb_type>(0u);
  196. }
  197. else
  198. m_exponent = o.exponent();
  199. return *this;
  200. }
  201. // rvalue copy:
  202. template <class A, class E, E MinE, E MaxE>
  203. cpp_bin_float& operator=(cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>&& o) noexcept(noexcept(std::declval<rep_type&>() = std::declval<rep_type&&>()))
  204. {
  205. m_data = std::move(o.bits());
  206. m_sign = o.sign();
  207. if (o.exponent() == cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>::exponent_zero)
  208. m_exponent = exponent_zero;
  209. else if (o.exponent() == cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>::exponent_nan)
  210. m_exponent = exponent_nan;
  211. else if (o.exponent() == cpp_bin_float<Digits, DigitBase, A, E, MinE, MaxE>::exponent_infinity)
  212. m_exponent = exponent_infinity;
  213. else if (o.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
  214. {
  215. // Overflow:
  216. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
  217. bits() = static_cast<limb_type>(0u);
  218. }
  219. else if (o.exponent() < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
  220. {
  221. // Underflow:
  222. exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  223. bits() = static_cast<limb_type>(0u);
  224. }
  225. else
  226. m_exponent = o.exponent();
  227. return *this;
  228. }
  229. template <unsigned D, digit_base_type B, class A, class E, E MinE, E MaxE>
  230. cpp_bin_float& operator=(const cpp_bin_float<D, B, A, E, MinE, MaxE>& f)
  231. {
  232. switch (eval_fpclassify(f))
  233. {
  234. case FP_ZERO:
  235. m_data = limb_type(0);
  236. m_sign = f.sign();
  237. m_exponent = exponent_zero;
  238. break;
  239. case FP_NAN:
  240. m_data = limb_type(0);
  241. m_sign = false;
  242. m_exponent = exponent_nan;
  243. break;
  244. ;
  245. case FP_INFINITE:
  246. m_data = limb_type(0);
  247. m_sign = f.sign();
  248. m_exponent = exponent_infinity;
  249. break;
  250. default:
  251. typename cpp_bin_float<D, B, A, E, MinE, MaxE>::rep_type b(f.bits());
  252. this->exponent() = f.exponent() + (E)bit_count - (E)cpp_bin_float<D, B, A, E, MinE, MaxE>::bit_count;
  253. this->sign() = f.sign();
  254. copy_and_round(*this, b);
  255. }
  256. return *this;
  257. }
  258. #ifdef BOOST_MP_HAS_FLOAT128_SUPPORT
  259. template <class Float>
  260. typename std::enable_if<
  261. (number_category<Float>::value == number_kind_floating_point)
  262. //&& (std::numeric_limits<Float>::digits <= static_cast<int>(bit_count))
  263. && ((std::numeric_limits<Float>::radix == 2) || (std::is_same<Float, float128_type>::value)),
  264. cpp_bin_float&>::type
  265. operator=(const Float& f)
  266. #else
  267. template <class Float>
  268. typename std::enable_if<
  269. (number_category<Float>::value == number_kind_floating_point)
  270. //&& (std::numeric_limits<Float>::digits <= static_cast<int>(bit_count))
  271. && (std::numeric_limits<Float>::radix == 2),
  272. cpp_bin_float&>::type
  273. operator=(const Float& f)
  274. #endif
  275. {
  276. return assign_float(f);
  277. }
  278. #ifdef BOOST_MP_HAS_FLOAT128_SUPPORT
  279. template <class Float>
  280. typename std::enable_if<std::is_same<Float, float128_type>::value && (std::numeric_limits<Float>::digits > Digits), cpp_bin_float&>::type assign_float(Float f)
  281. {
  282. cpp_bin_float<113, DigitBase, Allocator, Exponent, MinExponent, MaxExponent> bf(f);
  283. return *this = bf;
  284. }
  285. template <class Float>
  286. typename std::enable_if<std::is_same<Float, float128_type>::value && (std::numeric_limits<Float>::digits <= Digits), cpp_bin_float&>::type assign_float(Float f)
  287. {
  288. using default_ops::eval_add;
  289. using bf_int_type = typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type;
  290. if (f == 0)
  291. {
  292. m_data = limb_type(0);
  293. m_sign = (signbitq(f) > 0);
  294. m_exponent = exponent_zero;
  295. return *this;
  296. }
  297. else if (isnanq(f))
  298. {
  299. m_data = limb_type(0);
  300. m_sign = false;
  301. m_exponent = exponent_nan;
  302. return *this;
  303. }
  304. else if (isinfq(f))
  305. {
  306. m_data = limb_type(0);
  307. m_sign = (f < 0);
  308. m_exponent = exponent_infinity;
  309. return *this;
  310. }
  311. if (f < 0)
  312. {
  313. *this = -f;
  314. this->negate();
  315. return *this;
  316. }
  317. using ui_type = typename std::tuple_element<0, unsigned_types>::type;
  318. m_data = static_cast<ui_type>(0u);
  319. m_sign = false;
  320. m_exponent = 0;
  321. constexpr std::ptrdiff_t bits = static_cast<Exponent>(sizeof(int) * CHAR_BIT - 1) < MaxExponent - 1 ? sizeof(int) * CHAR_BIT - 1 : 3;
  322. int e;
  323. f = frexpq(f, &e);
  324. while (f)
  325. {
  326. f = ldexpq(f, bits);
  327. e -= bits;
  328. int ipart = static_cast<int>(truncq(f));
  329. f -= ipart;
  330. m_exponent += bits;
  331. cpp_bin_float t;
  332. t = static_cast<bf_int_type>(ipart);
  333. eval_add(*this, t);
  334. }
  335. m_exponent += static_cast<Exponent>(e);
  336. if (m_exponent > max_exponent)
  337. {
  338. m_exponent = exponent_infinity;
  339. m_data = static_cast<ui_type>(0u);
  340. }
  341. else if (m_exponent < min_exponent)
  342. {
  343. m_exponent = exponent_zero;
  344. m_data = static_cast<ui_type>(0u);
  345. }
  346. return *this;
  347. }
  348. #endif
  349. #ifdef BOOST_MP_HAS_FLOAT128_SUPPORT
  350. template <class Float>
  351. typename std::enable_if<std::is_floating_point<Float>::value && !std::is_same<Float, float128_type>::value && (std::numeric_limits<Float>::digits > Digits), cpp_bin_float&>::type assign_float(Float f)
  352. #else
  353. template <class Float>
  354. typename std::enable_if<std::is_floating_point<Float>::value && (std::numeric_limits<Float>::digits > Digits), cpp_bin_float&>::type assign_float(Float f)
  355. #endif
  356. {
  357. cpp_bin_float<std::numeric_limits<Float>::digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent> bf(f);
  358. return *this = bf;
  359. }
  360. #ifdef BOOST_MP_HAS_FLOAT128_SUPPORT
  361. template <class Float>
  362. typename std::enable_if<std::is_floating_point<Float>::value && !std::is_same<Float, float128_type>::value && (std::numeric_limits<Float>::digits <= Digits), cpp_bin_float&>::type assign_float(Float f)
  363. #else
  364. template <class Float>
  365. typename std::enable_if<std::is_floating_point<Float>::value && (std::numeric_limits<Float>::digits <= Digits), cpp_bin_float&>::type assign_float(Float f)
  366. #endif
  367. {
  368. using std::frexp;
  369. using std::ldexp;
  370. using std::signbit;
  371. using default_ops::eval_add;
  372. using bf_int_type = typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type;
  373. switch (BOOST_MP_FPCLASSIFY(f))
  374. {
  375. case FP_ZERO:
  376. m_data = limb_type(0);
  377. m_sign = ((signbit)(f));
  378. m_exponent = exponent_zero;
  379. return *this;
  380. case FP_NAN:
  381. m_data = limb_type(0);
  382. m_sign = false;
  383. m_exponent = exponent_nan;
  384. return *this;
  385. case FP_INFINITE:
  386. m_data = limb_type(0);
  387. m_sign = (f < 0);
  388. m_exponent = exponent_infinity;
  389. return *this;
  390. default:
  391. break;
  392. }
  393. if (f < 0)
  394. {
  395. *this = -f;
  396. this->negate();
  397. return *this;
  398. }
  399. using ui_type = typename std::tuple_element<0, unsigned_types>::type;
  400. m_data = static_cast<ui_type>(0u);
  401. m_sign = false;
  402. m_exponent = 0;
  403. //
  404. // This code picks off the bits in f a few at a time and injects them into *this.
  405. // It does not do roundingm so we must have more digits precision in *this than
  406. // in the floating point value (the normal situation, unless we're emulating another
  407. // type like float16_t).
  408. //
  409. constexpr std::ptrdiff_t bits = static_cast<std::ptrdiff_t>(sizeof(int) * CHAR_BIT - 1) < static_cast<std::ptrdiff_t>(MaxExponent - 1) ? static_cast<std::ptrdiff_t>(sizeof(int) * CHAR_BIT - 1) : 3;
  410. int e;
  411. f = frexp(f, &e);
  412. while (f != static_cast<Float>(0.0F))
  413. {
  414. f = ldexp(f, bits);
  415. e -= static_cast<int>(bits);
  416. int ipart = boost::multiprecision::detail::itrunc(f);
  417. f -= static_cast<Float>(ipart);
  418. m_exponent += static_cast<exponent_type>(bits);
  419. cpp_bin_float t;
  420. t = static_cast<bf_int_type>(ipart);
  421. eval_add(*this, t);
  422. }
  423. m_exponent += static_cast<Exponent>(e);
  424. if (m_exponent > max_exponent)
  425. {
  426. m_exponent = exponent_infinity;
  427. m_data = static_cast<ui_type>(0u);
  428. }
  429. else if(m_exponent < min_exponent)
  430. {
  431. m_exponent = exponent_zero;
  432. m_data = static_cast<ui_type>(0u);
  433. }
  434. return *this;
  435. }
  436. template <class Float>
  437. typename std::enable_if<
  438. (number_category<Float>::value == number_kind_floating_point) && !std::is_floating_point<Float>::value && (number_category<Float>::value == number_kind_floating_point),
  439. cpp_bin_float&>::type
  440. assign_float(Float f)
  441. {
  442. using default_ops::eval_add;
  443. using default_ops::eval_convert_to;
  444. using default_ops::eval_get_sign;
  445. using default_ops::eval_subtract;
  446. using f_int_type = typename boost::multiprecision::detail::canonical<int, Float>::type ;
  447. using bf_int_type = typename boost::multiprecision::detail::canonical<int, cpp_bin_float>::type;
  448. switch (eval_fpclassify(f))
  449. {
  450. case FP_ZERO:
  451. m_data = limb_type(0);
  452. m_sign = (eval_get_sign(f) > 0);
  453. m_exponent = exponent_zero;
  454. return *this;
  455. case FP_NAN:
  456. m_data = limb_type(0);
  457. m_sign = false;
  458. m_exponent = exponent_nan;
  459. return *this;
  460. case FP_INFINITE:
  461. m_data = limb_type(0);
  462. m_sign = eval_get_sign(f) < 0;
  463. m_exponent = exponent_infinity;
  464. return *this;
  465. default:
  466. break;
  467. }
  468. if (eval_get_sign(f) < 0)
  469. {
  470. f.negate();
  471. assign_float(f);
  472. this->negate();
  473. return *this;
  474. }
  475. using ui_type = typename std::tuple_element<0, unsigned_types>::type;
  476. m_data = static_cast<ui_type>(0u);
  477. m_sign = false;
  478. m_exponent = 0;
  479. constexpr std::ptrdiff_t bits = sizeof(int) * CHAR_BIT - 1;
  480. int e;
  481. eval_frexp(f, f, &e);
  482. while (eval_get_sign(f) != 0)
  483. {
  484. eval_ldexp(f, f, bits);
  485. e -= bits;
  486. int ipart;
  487. eval_convert_to(&ipart, f);
  488. eval_subtract(f, static_cast<f_int_type>(ipart));
  489. m_exponent += bits;
  490. eval_add(*this, static_cast<bf_int_type>(ipart));
  491. }
  492. m_exponent += e;
  493. if (m_exponent > max_exponent)
  494. m_exponent = exponent_infinity;
  495. if (m_exponent < min_exponent)
  496. {
  497. m_data = limb_type(0u);
  498. m_exponent = exponent_zero;
  499. m_sign = (eval_get_sign(f) > 0);
  500. }
  501. else if (eval_get_sign(m_data) == 0)
  502. {
  503. m_exponent = exponent_zero;
  504. m_sign = (eval_get_sign(f) > 0);
  505. }
  506. return *this;
  507. }
  508. template <class B, expression_template_option et>
  509. cpp_bin_float& assign_float(const number<B, et>& f)
  510. {
  511. return assign_float(f.backend());
  512. }
  513. template <class I>
  514. typename std::enable_if<boost::multiprecision::detail::is_integral<I>::value, cpp_bin_float&>::type operator=(const I& i)
  515. {
  516. using default_ops::eval_bit_test;
  517. if (!i)
  518. {
  519. m_data = static_cast<limb_type>(0);
  520. m_exponent = exponent_zero;
  521. m_sign = false;
  522. }
  523. else
  524. {
  525. using ui_type = typename boost::multiprecision::detail::make_unsigned<I>::type ;
  526. ui_type fi = static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(i));
  527. using ar_type = typename boost::multiprecision::detail::canonical<ui_type, rep_type>::type;
  528. m_data = static_cast<ar_type>(fi);
  529. std::size_t shift = msb(fi);
  530. if (shift > max_exponent)
  531. {
  532. m_exponent = exponent_infinity;
  533. m_data = static_cast<limb_type>(0);
  534. }
  535. else if (shift >= bit_count)
  536. {
  537. m_exponent = static_cast<Exponent>(shift);
  538. m_data = static_cast<ar_type>(fi >> (shift + 1 - bit_count));
  539. }
  540. else
  541. {
  542. m_exponent = static_cast<Exponent>(shift);
  543. eval_left_shift(m_data, bit_count - shift - 1);
  544. }
  545. BOOST_MP_ASSERT((m_exponent == exponent_infinity) || eval_bit_test(m_data, bit_count - 1));
  546. m_sign = detail::is_negative(i);
  547. }
  548. return *this;
  549. }
  550. cpp_bin_float& operator=(const char* s);
  551. void swap(cpp_bin_float& o) noexcept
  552. {
  553. m_data.swap(o.m_data);
  554. std::swap(m_exponent, o.m_exponent);
  555. std::swap(m_sign, o.m_sign);
  556. }
  557. std::string str(std::streamsize dig, std::ios_base::fmtflags f) const;
  558. void negate()
  559. {
  560. if (m_exponent != exponent_nan)
  561. m_sign = !m_sign;
  562. }
  563. int compare(const cpp_bin_float& o) const noexcept
  564. {
  565. if (m_sign != o.m_sign)
  566. return (m_exponent == exponent_zero) && (m_exponent == o.m_exponent) ? 0 : m_sign ? -1 : 1;
  567. int result;
  568. if (m_exponent == exponent_nan)
  569. return -1;
  570. else if (m_exponent != o.m_exponent)
  571. {
  572. if (m_exponent == exponent_zero)
  573. result = -1;
  574. else if (o.m_exponent == exponent_zero)
  575. result = 1;
  576. else
  577. result = m_exponent > o.m_exponent ? 1 : -1;
  578. }
  579. else
  580. result = m_data.compare(o.m_data);
  581. if (m_sign)
  582. result = -result;
  583. return result;
  584. }
  585. template <class A>
  586. int compare(const A& o) const noexcept
  587. {
  588. cpp_bin_float b;
  589. b = o;
  590. return compare(b);
  591. }
  592. rep_type& bits() { return m_data; }
  593. const rep_type& bits() const { return m_data; }
  594. exponent_type& exponent() { return m_exponent; }
  595. const exponent_type& exponent() const { return m_exponent; }
  596. bool& sign() { return m_sign; }
  597. const bool& sign() const { return m_sign; }
  598. void check_invariants()
  599. {
  600. using default_ops::eval_bit_test;
  601. using default_ops::eval_is_zero;
  602. if ((m_exponent <= max_exponent) && (m_exponent >= min_exponent))
  603. {
  604. BOOST_MP_ASSERT(eval_bit_test(m_data, bit_count - 1));
  605. }
  606. else
  607. {
  608. BOOST_MP_ASSERT(m_exponent > max_exponent);
  609. BOOST_MP_ASSERT(m_exponent <= exponent_nan);
  610. BOOST_MP_ASSERT(eval_is_zero(m_data));
  611. }
  612. }
  613. #ifndef BOOST_MP_STANDALONE
  614. template <class Archive>
  615. void serialize(Archive& ar, const unsigned int /*version*/)
  616. {
  617. ar& boost::make_nvp("data", m_data);
  618. ar& boost::make_nvp("exponent", m_exponent);
  619. ar& boost::make_nvp("sign", m_sign);
  620. }
  621. #endif
  622. };
  623. #ifdef BOOST_MSVC
  624. #pragma warning(pop)
  625. #endif
  626. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class Int>
  627. inline void copy_and_round(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, Int& arg, std::ptrdiff_t bits_to_keep = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
  628. {
  629. // Precondition: exponent of res must have been set before this function is called
  630. // as we may need to adjust it based on how many bits_to_keep in arg are set.
  631. using default_ops::eval_bit_test;
  632. using default_ops::eval_get_sign;
  633. using default_ops::eval_increment;
  634. using default_ops::eval_left_shift;
  635. using default_ops::eval_lsb;
  636. using default_ops::eval_msb;
  637. using default_ops::eval_right_shift;
  638. // cancellation may have resulted in arg being all zeros:
  639. if (eval_get_sign(arg) == 0)
  640. {
  641. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  642. res.sign() = false;
  643. res.bits() = static_cast<limb_type>(0u);
  644. return;
  645. }
  646. std::ptrdiff_t msb = static_cast<std::ptrdiff_t>(eval_msb(arg));
  647. if (static_cast<std::ptrdiff_t >(bits_to_keep) > msb + 1)
  648. {
  649. // Must have had cancellation in subtraction,
  650. // or be converting from a narrower type, so shift left:
  651. res.bits() = arg;
  652. eval_left_shift(res.bits(), static_cast<double_limb_type>(bits_to_keep - msb - 1));
  653. res.exponent() -= static_cast<Exponent>(bits_to_keep - msb - 1);
  654. }
  655. else if (static_cast<std::ptrdiff_t >(bits_to_keep) < msb + 1)
  656. {
  657. // We have more bits_to_keep than we need, so round as required,
  658. // first get the rounding bit:
  659. bool roundup = eval_bit_test(arg, static_cast<std::size_t>(msb - bits_to_keep));
  660. // Then check for a tie:
  661. if (roundup && (msb - bits_to_keep == static_cast<std::ptrdiff_t>(eval_lsb(arg))))
  662. {
  663. // Ties round towards even:
  664. if (!eval_bit_test(arg, static_cast<std::size_t>(msb - bits_to_keep + 1)))
  665. roundup = false;
  666. }
  667. // Shift off the bits_to_keep we don't need:
  668. eval_right_shift(arg, static_cast<double_limb_type>(msb - bits_to_keep + 1));
  669. res.exponent() += static_cast<Exponent>(msb - bits_to_keep + 1);
  670. if (roundup)
  671. {
  672. eval_increment(arg);
  673. if (bits_to_keep)
  674. {
  675. if (eval_bit_test(arg, static_cast<std::size_t>(bits_to_keep)))
  676. {
  677. // This happens very very rairly, all the bits left after
  678. // truncation must be 1's and we're rounding up an order of magnitude:
  679. eval_right_shift(arg, 1u);
  680. ++res.exponent();
  681. }
  682. }
  683. else
  684. {
  685. // We get here when bits_to_keep is zero but we're rounding up,
  686. // as a result we end up with a single digit that is a 1:
  687. ++bits_to_keep;
  688. }
  689. }
  690. if (bits_to_keep != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
  691. {
  692. // Normalize result when we're rounding to fewer bits than we can hold, only happens in conversions
  693. // to narrower types:
  694. eval_left_shift(arg, static_cast<double_limb_type>(static_cast<std::ptrdiff_t>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - bits_to_keep));
  695. res.exponent() -= static_cast<Exponent>(static_cast<std::ptrdiff_t>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - bits_to_keep);
  696. }
  697. res.bits() = arg;
  698. }
  699. else
  700. {
  701. res.bits() = arg;
  702. }
  703. if (!bits_to_keep && !res.bits().limbs()[0])
  704. {
  705. // We're keeping zero bits and did not round up, so result is zero:
  706. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  707. return;
  708. }
  709. // Result must be normalized:
  710. BOOST_MP_ASSERT(((std::ptrdiff_t )eval_msb(res.bits()) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
  711. if (res.exponent() > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent)
  712. {
  713. // Overflow:
  714. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
  715. res.bits() = static_cast<limb_type>(0u);
  716. }
  717. else if (res.exponent() < cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent)
  718. {
  719. // Underflow:
  720. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  721. res.bits() = static_cast<limb_type>(0u);
  722. }
  723. }
  724. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class BinFloat2, class BinFloat3>
  725. inline void do_eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  726. const BinFloat2& a, const BinFloat3& b)
  727. {
  728. if (a.exponent() < b.exponent())
  729. {
  730. bool s = a.sign();
  731. do_eval_add(res, b, a);
  732. if (res.sign() != s)
  733. res.negate();
  734. return;
  735. }
  736. using default_ops::eval_add;
  737. using default_ops::eval_bit_test;
  738. using exponent_type = typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type;
  739. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
  740. // Special cases first:
  741. switch (a.exponent())
  742. {
  743. case BinFloat2::exponent_zero:
  744. {
  745. bool s = a.sign();
  746. res = b;
  747. res.sign() = s;
  748. return;
  749. }
  750. case BinFloat2::exponent_infinity:
  751. if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan)
  752. res = b;
  753. else
  754. res = a;
  755. return; // result is still infinite.
  756. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  757. res = a;
  758. return; // result is still a NaN.
  759. default:
  760. break;
  761. }
  762. switch (b.exponent())
  763. {
  764. case BinFloat3::exponent_zero:
  765. res = a;
  766. return;
  767. case BinFloat3::exponent_infinity:
  768. res = b;
  769. if (res.sign())
  770. res.negate();
  771. return; // result is infinite.
  772. case BinFloat3::exponent_nan:
  773. res = b;
  774. return; // result is a NaN.
  775. default:
  776. break;
  777. }
  778. static_assert((std::numeric_limits<exponent_type>::max)() - cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count > cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent, "Exponent range check failed");
  779. bool s = a.sign();
  780. dt = a.bits();
  781. if (a.exponent() > (std::ptrdiff_t )cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + b.exponent())
  782. {
  783. res.exponent() = a.exponent();
  784. }
  785. else
  786. {
  787. exponent_type e_diff = a.exponent() - b.exponent();
  788. BOOST_MP_ASSERT(e_diff >= 0);
  789. eval_left_shift(dt, static_cast<double_limb_type>(e_diff));
  790. res.exponent() = a.exponent() - e_diff;
  791. eval_add(dt, b.bits());
  792. }
  793. copy_and_round(res, dt);
  794. res.check_invariants();
  795. if (res.sign() != s)
  796. res.negate();
  797. }
  798. template <class BinFloat1, class BinFloat2, class BinFloat3>
  799. inline void do_eval_subtract(BinFloat1& res, const BinFloat2& a, const BinFloat3& b)
  800. {
  801. using default_ops::eval_bit_test;
  802. using default_ops::eval_decrement;
  803. using default_ops::eval_subtract;
  804. typename BinFloat1::double_rep_type dt;
  805. // Special cases first:
  806. switch (a.exponent())
  807. {
  808. case BinFloat2::exponent_zero:
  809. if (b.exponent() == BinFloat3::exponent_nan)
  810. res = std::numeric_limits<number<BinFloat1> >::quiet_NaN().backend();
  811. else
  812. {
  813. bool s = a.sign();
  814. res = b;
  815. if (res.exponent() == BinFloat1::exponent_zero)
  816. res.sign() = false;
  817. else if (res.sign() == s)
  818. res.negate();
  819. }
  820. return;
  821. case BinFloat2::exponent_infinity:
  822. if ((b.exponent() == BinFloat3::exponent_nan) || (b.exponent() == BinFloat3::exponent_infinity))
  823. res = std::numeric_limits<number<BinFloat1> >::quiet_NaN().backend();
  824. else
  825. res = a;
  826. return;
  827. case BinFloat2::exponent_nan:
  828. res = a;
  829. return; // result is still a NaN.
  830. default:
  831. break;
  832. }
  833. switch (b.exponent())
  834. {
  835. case BinFloat3::exponent_zero:
  836. res = a;
  837. return;
  838. case BinFloat3::exponent_infinity:
  839. res.exponent() = BinFloat1::exponent_infinity;
  840. res.sign() = !a.sign();
  841. res.bits() = static_cast<limb_type>(0u);
  842. return; // result is a NaN.
  843. case BinFloat3::exponent_nan:
  844. res = b;
  845. return; // result is still a NaN.
  846. default:
  847. break;
  848. }
  849. bool s = a.sign();
  850. if ((a.exponent() > b.exponent()) || ((a.exponent() == b.exponent()) && a.bits().compare(b.bits()) >= 0))
  851. {
  852. dt = a.bits();
  853. if (a.exponent() <= (std::ptrdiff_t )BinFloat1::bit_count + b.exponent())
  854. {
  855. typename BinFloat1::exponent_type e_diff = a.exponent() - b.exponent();
  856. eval_left_shift(dt, static_cast<double_limb_type>(e_diff));
  857. res.exponent() = a.exponent() - e_diff;
  858. eval_subtract(dt, b.bits());
  859. }
  860. else if (a.exponent() == (std::ptrdiff_t )BinFloat1::bit_count + b.exponent() + 1)
  861. {
  862. if ((eval_lsb(a.bits()) == BinFloat1::bit_count - 1)
  863. && (eval_lsb(b.bits()) != BinFloat1::bit_count - 1))
  864. {
  865. eval_left_shift(dt, 1);
  866. eval_decrement(dt);
  867. res.exponent() = a.exponent() - 1;
  868. }
  869. else
  870. res.exponent() = a.exponent();
  871. }
  872. else
  873. res.exponent() = a.exponent();
  874. }
  875. else
  876. {
  877. dt = b.bits();
  878. if (b.exponent() <= (std::ptrdiff_t )BinFloat1::bit_count + a.exponent())
  879. {
  880. typename BinFloat1::exponent_type e_diff = a.exponent() - b.exponent();
  881. eval_left_shift(dt, static_cast<double_limb_type>(-e_diff));
  882. res.exponent() = b.exponent() + e_diff;
  883. eval_subtract(dt, a.bits());
  884. }
  885. else if (b.exponent() == (std::ptrdiff_t )BinFloat1::bit_count + a.exponent() + 1)
  886. {
  887. if ((eval_lsb(a.bits()) != BinFloat1::bit_count - 1)
  888. && eval_lsb(b.bits()))
  889. {
  890. eval_left_shift(dt, 1);
  891. eval_decrement(dt);
  892. res.exponent() = b.exponent() - 1;
  893. }
  894. else
  895. res.exponent() = b.exponent();
  896. }
  897. else
  898. res.exponent() = b.exponent();
  899. s = !s;
  900. }
  901. copy_and_round(res, dt);
  902. if (res.exponent() == BinFloat1::exponent_zero)
  903. res.sign() = false;
  904. else if (res.sign() != s)
  905. res.negate();
  906. res.check_invariants();
  907. }
  908. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  909. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2,
  910. class Allocator3, class Exponent3, Exponent MinE3, Exponent MaxE3>
  911. inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  912. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a,
  913. const cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>& b)
  914. {
  915. if (a.sign() == b.sign())
  916. do_eval_add(res, a, b);
  917. else
  918. do_eval_subtract(res, a, b);
  919. }
  920. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  921. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2>
  922. inline void eval_add(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  923. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a)
  924. {
  925. return eval_add(res, res, a);
  926. }
  927. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  928. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2,
  929. class Allocator3, class Exponent3, Exponent MinE3, Exponent MaxE3>
  930. inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  931. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a,
  932. const cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>& b)
  933. {
  934. if (a.sign() != b.sign())
  935. do_eval_add(res, a, b);
  936. else
  937. do_eval_subtract(res, a, b);
  938. }
  939. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  940. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2>
  941. inline void eval_subtract(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  942. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a)
  943. {
  944. return eval_subtract(res, res, a);
  945. }
  946. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  947. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2,
  948. class Allocator3, class Exponent3, Exponent MinE3, Exponent MaxE3>
  949. inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  950. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a,
  951. const cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>& b)
  952. {
  953. using default_ops::eval_bit_test;
  954. using default_ops::eval_multiply;
  955. // Special cases first:
  956. switch (a.exponent())
  957. {
  958. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_zero:
  959. {
  960. if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_nan)
  961. res = b;
  962. else if (b.exponent() == cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_infinity)
  963. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  964. else
  965. {
  966. bool s = a.sign() != b.sign();
  967. res = a;
  968. res.sign() = s;
  969. }
  970. return;
  971. }
  972. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_infinity:
  973. switch (b.exponent())
  974. {
  975. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_zero:
  976. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  977. break;
  978. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_nan:
  979. res = b;
  980. break;
  981. default:
  982. bool s = a.sign() != b.sign();
  983. res = a;
  984. res.sign() = s;
  985. break;
  986. }
  987. return;
  988. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_nan:
  989. res = a;
  990. return;
  991. default:
  992. break;
  993. }
  994. if (b.exponent() > cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::max_exponent)
  995. {
  996. bool s = a.sign() != b.sign();
  997. res = b;
  998. res.sign() = s;
  999. return;
  1000. }
  1001. if ((a.exponent() > 0) && (b.exponent() > 0))
  1002. {
  1003. if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + 2 - a.exponent() < b.exponent())
  1004. {
  1005. // We will certainly overflow:
  1006. bool s = a.sign() != b.sign();
  1007. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
  1008. res.sign() = s;
  1009. res.bits() = static_cast<limb_type>(0u);
  1010. return;
  1011. }
  1012. }
  1013. if ((a.exponent() < 0) && (b.exponent() < 0))
  1014. {
  1015. if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - 2 - a.exponent() > b.exponent())
  1016. {
  1017. // We will certainly underflow:
  1018. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  1019. res.sign() = a.sign() != b.sign();
  1020. res.bits() = static_cast<limb_type>(0u);
  1021. return;
  1022. }
  1023. }
  1024. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
  1025. eval_multiply(dt, a.bits(), b.bits());
  1026. res.exponent() = a.exponent() + b.exponent() - (Exponent)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count + 1;
  1027. copy_and_round(res, dt);
  1028. res.check_invariants();
  1029. res.sign() = a.sign() != b.sign();
  1030. }
  1031. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1032. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2>
  1033. inline void eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  1034. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a)
  1035. {
  1036. eval_multiply(res, res, a);
  1037. }
  1038. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1039. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2, class U>
  1040. inline typename std::enable_if<boost::multiprecision::detail::is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  1041. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a, const U& b)
  1042. {
  1043. using default_ops::eval_bit_test;
  1044. using default_ops::eval_multiply;
  1045. bool s = a.sign(); // saved for later in case a and res are the same object.
  1046. // Special cases first:
  1047. switch (a.exponent())
  1048. {
  1049. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_zero:
  1050. {
  1051. res = a;
  1052. res.sign() = s;
  1053. return;
  1054. }
  1055. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_infinity:
  1056. if (b == 0)
  1057. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1058. else
  1059. res = a;
  1060. return;
  1061. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_nan:
  1062. res = a;
  1063. return;
  1064. default:
  1065. break;
  1066. }
  1067. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type dt;
  1068. using canon_ui_type = typename boost::multiprecision::detail::canonical<U, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::type;
  1069. eval_multiply(dt, a.bits(), static_cast<canon_ui_type>(b));
  1070. res.exponent() = a.exponent();
  1071. copy_and_round(res, dt);
  1072. res.check_invariants();
  1073. res.sign() = s;
  1074. }
  1075. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
  1076. inline typename std::enable_if<boost::multiprecision::detail::is_unsigned<U>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const U& b)
  1077. {
  1078. eval_multiply(res, res, b);
  1079. }
  1080. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1081. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2, class S>
  1082. inline typename std::enable_if<boost::multiprecision::detail::is_signed<S>::value && boost::multiprecision::detail::is_integral<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  1083. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& a, const S& b)
  1084. {
  1085. using ui_type = typename boost::multiprecision::detail::make_unsigned<S>::type;
  1086. eval_multiply(res, a, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(b)));
  1087. if (b < 0)
  1088. res.negate();
  1089. }
  1090. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
  1091. inline typename std::enable_if<boost::multiprecision::detail::is_signed<S>::value && boost::multiprecision::detail::is_integral<S>::value>::type eval_multiply(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const S& b)
  1092. {
  1093. eval_multiply(res, res, b);
  1094. }
  1095. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1096. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2,
  1097. class Allocator3, class Exponent3, Exponent MinE3, Exponent MaxE3>
  1098. inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  1099. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& u,
  1100. const cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>& v)
  1101. {
  1102. #ifdef BOOST_MSVC
  1103. #pragma warning(push)
  1104. #pragma warning(disable : 6326) // comparison of two constants
  1105. #endif
  1106. using default_ops::eval_bit_test;
  1107. using default_ops::eval_get_sign;
  1108. using default_ops::eval_increment;
  1109. using default_ops::eval_qr;
  1110. using default_ops::eval_subtract;
  1111. //
  1112. // Special cases first:
  1113. //
  1114. switch (u.exponent())
  1115. {
  1116. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_zero:
  1117. {
  1118. switch (v.exponent())
  1119. {
  1120. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_zero:
  1121. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_nan:
  1122. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1123. return;
  1124. default:
  1125. break;
  1126. }
  1127. bool s = u.sign() != v.sign();
  1128. res = u;
  1129. res.sign() = s;
  1130. return;
  1131. }
  1132. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_infinity:
  1133. {
  1134. switch (v.exponent())
  1135. {
  1136. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_infinity:
  1137. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_nan:
  1138. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1139. return;
  1140. default:
  1141. break;
  1142. }
  1143. bool s = u.sign() != v.sign();
  1144. res = u;
  1145. res.sign() = s;
  1146. return;
  1147. }
  1148. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_nan:
  1149. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1150. return;
  1151. default:
  1152. break;
  1153. }
  1154. switch (v.exponent())
  1155. {
  1156. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_zero:
  1157. {
  1158. bool s = u.sign() != v.sign();
  1159. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
  1160. res.sign() = s;
  1161. return;
  1162. }
  1163. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_infinity:
  1164. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  1165. res.bits() = limb_type(0);
  1166. res.sign() = u.sign() != v.sign();
  1167. return;
  1168. case cpp_bin_float<Digits, DigitBase, Allocator3, Exponent3, MinE3, MaxE3>::exponent_nan:
  1169. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1170. return;
  1171. default:
  1172. break;
  1173. }
  1174. // We can scale u and v so that both are integers, then perform integer
  1175. // division to obtain quotient q and remainder r, such that:
  1176. //
  1177. // q * v + r = u
  1178. //
  1179. // and hence:
  1180. //
  1181. // q + r/v = u/v
  1182. //
  1183. // From this, assuming q has cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count
  1184. // bits we only need to determine whether
  1185. // r/v is less than, equal to, or greater than 0.5 to determine rounding -
  1186. // this we can do with a shift and comparison.
  1187. //
  1188. // We can set the exponent and sign of the result up front:
  1189. //
  1190. if ((v.exponent() < 0) && (u.exponent() > 0))
  1191. {
  1192. // Check for overflow:
  1193. if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent + v.exponent() < u.exponent() - 1)
  1194. {
  1195. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
  1196. res.sign() = u.sign() != v.sign();
  1197. res.bits() = static_cast<limb_type>(0u);
  1198. return;
  1199. }
  1200. }
  1201. else if ((v.exponent() > 0) && (u.exponent() < 0))
  1202. {
  1203. // Check for underflow:
  1204. if (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent + v.exponent() > u.exponent())
  1205. {
  1206. // We will certainly underflow:
  1207. res.exponent() = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  1208. res.sign() = u.sign() != v.sign();
  1209. res.bits() = static_cast<limb_type>(0u);
  1210. return;
  1211. }
  1212. }
  1213. res.exponent() = u.exponent() - v.exponent() - 1;
  1214. res.sign() = u.sign() != v.sign();
  1215. //
  1216. // Now get the quotient and remainder:
  1217. //
  1218. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), t2(v.bits()), q, r;
  1219. eval_left_shift(t, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count);
  1220. eval_qr(t, t2, q, r);
  1221. //
  1222. // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count"
  1223. // or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant
  1224. // bits in q.
  1225. //
  1226. constexpr unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  1227. if (eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
  1228. {
  1229. //
  1230. // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 bits,
  1231. // so we already have rounding info,
  1232. // we just need to changes things if the last bit is 1 and either the
  1233. // remainder is non-zero (ie we do not have a tie) or the quotient would
  1234. // be odd if it were shifted to the correct number of bits (ie a tiebreak).
  1235. //
  1236. BOOST_MP_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count));
  1237. if ((q.limbs()[0] & 1u) && (eval_get_sign(r) || (q.limbs()[0] & 2u)))
  1238. {
  1239. eval_increment(q);
  1240. }
  1241. }
  1242. else
  1243. {
  1244. //
  1245. // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" bits in q.
  1246. // Get rounding info, which we can get by comparing 2r with v.
  1247. // We want to call copy_and_round to handle rounding and general cleanup,
  1248. // so we'll left shift q and add some fake digits on the end to represent
  1249. // how we'll be rounding.
  1250. //
  1251. using local_exponent_type = typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type;
  1252. BOOST_MP_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
  1253. constexpr unsigned lshift = (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits) ? 2 : limb_bits;
  1254. eval_left_shift(q, lshift);
  1255. res.exponent() -= static_cast<local_exponent_type>(lshift);
  1256. eval_left_shift(r, 1u);
  1257. int c = r.compare(v.bits());
  1258. if (c == 0)
  1259. q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1);
  1260. else if (c > 0)
  1261. q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u);
  1262. }
  1263. copy_and_round(res, q);
  1264. #ifdef BOOST_MSVC
  1265. #pragma warning(pop)
  1266. #endif
  1267. }
  1268. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1269. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2>
  1270. inline void eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  1271. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& arg)
  1272. {
  1273. eval_divide(res, res, arg);
  1274. }
  1275. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1276. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2, class U>
  1277. inline typename std::enable_if<boost::multiprecision::detail::is_unsigned<U>::value && (std::numeric_limits<U>::digits <= Digits)>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  1278. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& u, const U& v)
  1279. {
  1280. #ifdef BOOST_MSVC
  1281. #pragma warning(push)
  1282. #pragma warning(disable : 6326) // comparison of two constants
  1283. #endif
  1284. using default_ops::eval_bit_test;
  1285. using default_ops::eval_get_sign;
  1286. using default_ops::eval_increment;
  1287. using default_ops::eval_qr;
  1288. using default_ops::eval_subtract;
  1289. //
  1290. // Special cases first:
  1291. //
  1292. switch (u.exponent())
  1293. {
  1294. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_zero:
  1295. {
  1296. if (v == 0)
  1297. {
  1298. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1299. return;
  1300. }
  1301. bool s = u.sign() != (v < 0);
  1302. res = u;
  1303. res.sign() = s;
  1304. return;
  1305. }
  1306. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_infinity:
  1307. res = u;
  1308. return;
  1309. case cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>::exponent_nan:
  1310. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1311. return;
  1312. default:
  1313. break;
  1314. }
  1315. if (v == 0)
  1316. {
  1317. bool s = u.sign();
  1318. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
  1319. res.sign() = s;
  1320. return;
  1321. }
  1322. // We can scale u and v so that both are integers, then perform integer
  1323. // division to obtain quotient q and remainder r, such that:
  1324. //
  1325. // q * v + r = u
  1326. //
  1327. // and hence:
  1328. //
  1329. // q + r/v = u/v
  1330. //
  1331. // From this, assuming q has "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, we only need to determine whether
  1332. // r/v is less than, equal to, or greater than 0.5 to determine rounding -
  1333. // this we can do with a shift and comparison.
  1334. //
  1335. // We can set the exponent and sign of the result up front:
  1336. //
  1337. std::ptrdiff_t gb = static_cast<std::ptrdiff_t>(msb(v));
  1338. res.exponent() = u.exponent() - static_cast<Exponent>(gb) - static_cast<Exponent>(1);
  1339. res.sign() = u.sign();
  1340. //
  1341. // Now get the quotient and remainder:
  1342. //
  1343. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(u.bits()), q, r;
  1344. eval_left_shift(t, static_cast<double_limb_type>(gb + 1));
  1345. eval_qr(t, number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v), q, r);
  1346. //
  1347. // We now have either "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" or "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1" significant cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q.
  1348. //
  1349. constexpr unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  1350. if (eval_bit_test(q, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
  1351. {
  1352. //
  1353. // OK we have cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count+1 cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count, so we already have rounding info,
  1354. // we just need to changes things if the last bit is 1 and the
  1355. // remainder is non-zero (ie we do not have a tie).
  1356. //
  1357. BOOST_MP_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count));
  1358. if ((q.limbs()[0] & 1u) && eval_get_sign(r))
  1359. {
  1360. eval_increment(q);
  1361. }
  1362. }
  1363. else
  1364. {
  1365. //
  1366. // We have exactly "cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count" cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in q.
  1367. // Get rounding info, which we can get by comparing 2r with v.
  1368. // We want to call copy_and_round to handle rounding and general cleanup,
  1369. // so we'll left shift q and add some fake cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count on the end to represent
  1370. // how we'll be rounding.
  1371. //
  1372. using local_exponent_type = typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type;
  1373. BOOST_MP_ASSERT((eval_msb(q) == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1));
  1374. constexpr unsigned lshift = cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count < limb_bits ? 2 : limb_bits;
  1375. eval_left_shift(q, lshift);
  1376. res.exponent() -= static_cast<local_exponent_type>(lshift);
  1377. eval_left_shift(r, 1u);
  1378. int c = r.compare(number<typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type>::canonical_value(v));
  1379. if (c == 0)
  1380. q.limbs()[0] |= static_cast<limb_type>(1u) << (lshift - 1);
  1381. else if (c > 0)
  1382. q.limbs()[0] |= (static_cast<limb_type>(1u) << (lshift - 1)) + static_cast<limb_type>(1u);
  1383. }
  1384. copy_and_round(res, q);
  1385. #ifdef BOOST_MSVC
  1386. #pragma warning(pop)
  1387. #endif
  1388. }
  1389. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class U>
  1390. inline typename std::enable_if<boost::multiprecision::detail::is_unsigned<U>::value && (std::numeric_limits<U>::digits <= Digits)>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const U& v)
  1391. {
  1392. eval_divide(res, res, v);
  1393. }
  1394. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1395. class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2, class S>
  1396. inline typename std::enable_if<boost::multiprecision::detail::is_signed<S>::value && boost::multiprecision::detail::is_integral<S>::value && (std::numeric_limits<S>::digits <= Digits)>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res,
  1397. const cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2>& u, const S& v)
  1398. {
  1399. using ui_type = typename boost::multiprecision::detail::make_unsigned<S>::type;
  1400. eval_divide(res, u, static_cast<ui_type>(boost::multiprecision::detail::unsigned_abs(v)));
  1401. if (v < 0)
  1402. res.negate();
  1403. }
  1404. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class S>
  1405. inline typename std::enable_if<boost::multiprecision::detail::is_signed<S>::value && boost::multiprecision::detail::is_integral<S>::value && (std::numeric_limits<S>::digits <= Digits)>::type eval_divide(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const S& v)
  1406. {
  1407. eval_divide(res, res, v);
  1408. }
  1409. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1410. inline int eval_get_sign(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1411. {
  1412. return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero ? 0 : arg.sign() ? -1 : 1;
  1413. }
  1414. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1415. inline bool eval_is_zero(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1416. {
  1417. return arg.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero;
  1418. }
  1419. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1420. inline bool eval_eq(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& a, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& b)
  1421. {
  1422. if (a.exponent() == b.exponent())
  1423. {
  1424. if (a.exponent() == cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero)
  1425. return true;
  1426. return (a.sign() == b.sign()) && (a.bits().compare(b.bits()) == 0) && (a.exponent() != cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan);
  1427. }
  1428. return false;
  1429. }
  1430. template <class I, unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1431. inline void convert_to_signed_int(I* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1432. {
  1433. static constexpr int digits = std::numeric_limits<I>::is_specialized ? std::numeric_limits<I>::digits : sizeof(I) * CHAR_BIT - 1;
  1434. static constexpr I max_val = std::numeric_limits<I>::is_specialized ? (std::numeric_limits<I>::max)() : (((I(1) << (sizeof(I) * CHAR_BIT - 2)) - 1) << 1) + 1;
  1435. static constexpr I min_val = std::numeric_limits<I>::is_specialized ? (std::numeric_limits<I>::min)() : -max_val - 1;
  1436. switch (arg.exponent())
  1437. {
  1438. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1439. *res = 0;
  1440. return;
  1441. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1442. BOOST_MP_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer."));
  1443. return;
  1444. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1445. *res = max_val;
  1446. if (arg.sign())
  1447. *res = -*res;
  1448. return;
  1449. default:
  1450. break;
  1451. }
  1452. using shift_type = typename std::conditional<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type;
  1453. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type man(arg.bits());
  1454. shift_type shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent();
  1455. if (shift > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
  1456. {
  1457. *res = 0;
  1458. return;
  1459. }
  1460. if (arg.sign() && (arg.compare(min_val) <= 0))
  1461. {
  1462. *res = min_val;
  1463. return;
  1464. }
  1465. else if (!arg.sign() && (arg.compare(max_val) >= 0))
  1466. {
  1467. *res = max_val;
  1468. return;
  1469. }
  1470. if (shift < 0)
  1471. {
  1472. if (static_cast<int>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - static_cast<int>(shift) <= digits)
  1473. {
  1474. // We have more bits in long_long_type than the float, so it's OK to left shift:
  1475. eval_convert_to(res, man);
  1476. *res <<= -shift;
  1477. }
  1478. else
  1479. {
  1480. *res = (std::numeric_limits<I>::max)();
  1481. return;
  1482. }
  1483. }
  1484. else
  1485. {
  1486. eval_right_shift(man, static_cast<double_limb_type>(shift));
  1487. eval_convert_to(res, man);
  1488. }
  1489. if (arg.sign())
  1490. {
  1491. *res = -*res;
  1492. }
  1493. }
  1494. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1495. inline void eval_convert_to(long long* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1496. {
  1497. convert_to_signed_int(res, arg);
  1498. }
  1499. #ifdef BOOST_HAS_INT128
  1500. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1501. inline void eval_convert_to(int128_type* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1502. {
  1503. convert_to_signed_int(res, arg);
  1504. }
  1505. #endif
  1506. template <class I, unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1507. inline void convert_to_unsigned_int(I* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1508. {
  1509. static constexpr int digits = std::numeric_limits<I>::is_specialized ? std::numeric_limits<I>::digits : sizeof(I) * CHAR_BIT;
  1510. static constexpr I max_val = std::numeric_limits<I>::is_specialized ? (std::numeric_limits<I>::max)() : ~static_cast<I>(0);
  1511. switch (arg.exponent())
  1512. {
  1513. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1514. *res = 0;
  1515. return;
  1516. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1517. BOOST_MP_THROW_EXCEPTION(std::runtime_error("Could not convert NaN to integer."));
  1518. #if defined(BOOST_NO_EXCEPTIONS)
  1519. *res = max_val;
  1520. return;
  1521. #endif
  1522. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1523. *res = max_val;
  1524. return;
  1525. default:
  1526. break;
  1527. }
  1528. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type man(arg.bits());
  1529. using shift_type = typename std::conditional<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type;
  1530. const shift_type shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1 - arg.exponent();
  1531. if (arg.sign())
  1532. {
  1533. using si_type = typename boost::multiprecision::detail::make_signed<I>::type;
  1534. si_type val;
  1535. convert_to_signed_int(&val, arg);
  1536. *res = static_cast<I>(val);
  1537. return;
  1538. }
  1539. else if (shift > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1)
  1540. {
  1541. *res = 0;
  1542. return;
  1543. }
  1544. else if (arg.compare(max_val) >= 0)
  1545. {
  1546. *res = max_val;
  1547. return;
  1548. }
  1549. else if (shift < 0)
  1550. {
  1551. if ((shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - shift <= (shift_type)digits)
  1552. {
  1553. // We have more bits in ulong_long_type than the float, so it's OK to left shift:
  1554. eval_convert_to(res, man);
  1555. *res <<= -shift;
  1556. return;
  1557. }
  1558. *res = max_val;
  1559. return;
  1560. }
  1561. eval_right_shift(man, static_cast<unsigned>(shift));
  1562. eval_convert_to(res, man);
  1563. }
  1564. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1565. inline void eval_convert_to(unsigned long long* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1566. {
  1567. convert_to_unsigned_int(res, arg);
  1568. }
  1569. #ifdef BOOST_HAS_INT128
  1570. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1571. inline void eval_convert_to(uint128_type* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1572. {
  1573. convert_to_unsigned_int(res, arg);
  1574. }
  1575. #endif
  1576. template <class Float, unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1577. inline typename std::enable_if<std::is_floating_point<Float>::value>::type eval_convert_to(Float* res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& original_arg)
  1578. {
  1579. using conv_type = cpp_bin_float<std::numeric_limits<Float>::digits, digit_base_2, void, Exponent, MinE, MaxE>;
  1580. using common_exp_type = typename std::common_type<typename conv_type::exponent_type, int>::type;
  1581. static constexpr int float_digits = boost::multiprecision::detail::is_float128<Float>::value ? 113 : std::numeric_limits<Float>::digits;
  1582. BOOST_MP_FLOAT128_USING using std::ldexp;
  1583. //
  1584. // Special cases first:
  1585. //
  1586. switch (original_arg.exponent())
  1587. {
  1588. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1589. *res = 0;
  1590. if (original_arg.sign())
  1591. *res = -*res;
  1592. return;
  1593. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1594. BOOST_IF_CONSTEXPR(boost::multiprecision::detail::is_float128<Float>::value)
  1595. {
  1596. *res = static_cast<Float>(std::numeric_limits<double>::quiet_NaN());
  1597. }
  1598. else
  1599. {
  1600. *res = std::numeric_limits<Float>::quiet_NaN();
  1601. }
  1602. return;
  1603. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1604. BOOST_IF_CONSTEXPR(boost::multiprecision::detail::is_float128<Float>::value)
  1605. {
  1606. *res = static_cast<Float>((std::numeric_limits<double>::infinity)());
  1607. }
  1608. else
  1609. {
  1610. *res = (std::numeric_limits<Float>::infinity)();
  1611. }
  1612. if (original_arg.sign())
  1613. *res = -*res;
  1614. return;
  1615. default:
  1616. break;
  1617. }
  1618. //
  1619. // Check for super large exponent that must be converted to infinity:
  1620. //
  1621. if (original_arg.exponent() > (boost::multiprecision::detail::is_float128<Float>::value ? 16384 : std::numeric_limits<Float>::max_exponent))
  1622. {
  1623. BOOST_IF_CONSTEXPR(boost::multiprecision::detail::is_float128<Float>::value)
  1624. {
  1625. *res = static_cast<Float>(std::numeric_limits<double>::infinity());
  1626. }
  1627. else
  1628. {
  1629. *res = std::numeric_limits<Float>::has_infinity ? std::numeric_limits<Float>::infinity() : (std::numeric_limits<Float>::max)();
  1630. }
  1631. if (original_arg.sign())
  1632. *res = -*res;
  1633. return;
  1634. }
  1635. //
  1636. // Figure out how many digits we will have in our result,
  1637. // allowing for a possibly denormalized result:
  1638. //
  1639. common_exp_type digits_to_round_to = float_digits;
  1640. if (original_arg.exponent() < std::numeric_limits<Float>::min_exponent - 1)
  1641. {
  1642. common_exp_type diff = original_arg.exponent();
  1643. diff -= boost::multiprecision::detail::is_float128<Float>::value ? -16382 : std::numeric_limits<Float>::min_exponent - 1;
  1644. digits_to_round_to += diff;
  1645. }
  1646. if (digits_to_round_to < 0)
  1647. {
  1648. // Result must be zero:
  1649. *res = 0;
  1650. if (original_arg.sign())
  1651. *res = -*res;
  1652. return;
  1653. }
  1654. //
  1655. // Perform rounding first, then afterwards extract the digits:
  1656. //
  1657. cpp_bin_float<static_cast<unsigned>(float_digits), digit_base_2, Allocator, Exponent, 0, 0> arg;
  1658. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::rep_type bits(original_arg.bits());
  1659. arg.exponent() = original_arg.exponent();
  1660. copy_and_round(arg, bits, (std::ptrdiff_t)digits_to_round_to);
  1661. common_exp_type e = arg.exponent();
  1662. e -= static_cast<common_exp_type>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1;
  1663. constexpr std::size_t limbs_needed = static_cast<std::size_t>(float_digits) / (sizeof(*arg.bits().limbs()) * CHAR_BIT) + (static_cast<std::size_t>(float_digits) % (sizeof(*arg.bits().limbs()) * CHAR_BIT) ? 1 : 0);
  1664. std::size_t first_limb_needed = arg.bits().size() - limbs_needed;
  1665. *res = 0;
  1666. e += static_cast<common_exp_type>(first_limb_needed * sizeof(*arg.bits().limbs()) * CHAR_BIT);
  1667. while (first_limb_needed < arg.bits().size())
  1668. {
  1669. *res += ldexp(static_cast<Float>(arg.bits().limbs()[first_limb_needed]), static_cast<int>(e));
  1670. ++first_limb_needed;
  1671. e += static_cast<common_exp_type>(sizeof(*arg.bits().limbs()) * CHAR_BIT);
  1672. }
  1673. if (original_arg.sign())
  1674. *res = -*res;
  1675. }
  1676. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1677. inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, Exponent* e)
  1678. {
  1679. switch (arg.exponent())
  1680. {
  1681. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1682. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1683. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1684. *e = 0;
  1685. res = arg;
  1686. return;
  1687. default:
  1688. break;
  1689. }
  1690. res = arg;
  1691. *e = arg.exponent() + 1;
  1692. res.exponent() = -1;
  1693. }
  1694. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
  1695. inline void eval_frexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I* pe)
  1696. {
  1697. Exponent e;
  1698. eval_frexp(res, arg, &e);
  1699. if ((e > (std::numeric_limits<I>::max)()) || (e < (std::numeric_limits<I>::min)()))
  1700. {
  1701. BOOST_MP_THROW_EXCEPTION(std::runtime_error("Exponent was outside of the range of the argument type to frexp."));
  1702. }
  1703. *pe = static_cast<I>(e);
  1704. }
  1705. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1706. inline void eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, Exponent e)
  1707. {
  1708. switch (arg.exponent())
  1709. {
  1710. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1711. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1712. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1713. res = arg;
  1714. return;
  1715. default:
  1716. break;
  1717. }
  1718. if ((e > 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent - e < arg.exponent()))
  1719. {
  1720. // Overflow:
  1721. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
  1722. res.sign() = arg.sign();
  1723. }
  1724. else if ((e < 0) && (cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent - e > arg.exponent()))
  1725. {
  1726. // Underflow:
  1727. res = limb_type(0);
  1728. }
  1729. else
  1730. {
  1731. res = arg;
  1732. res.exponent() += e;
  1733. }
  1734. }
  1735. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
  1736. inline typename std::enable_if<boost::multiprecision::detail::is_unsigned<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I e)
  1737. {
  1738. using si_type = typename boost::multiprecision::detail::make_signed<I>::type;
  1739. if (e > static_cast<I>((std::numeric_limits<si_type>::max)()))
  1740. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
  1741. else
  1742. eval_ldexp(res, arg, static_cast<si_type>(e));
  1743. }
  1744. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class I>
  1745. inline typename std::enable_if<boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value>::type eval_ldexp(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg, I e)
  1746. {
  1747. if ((e > (std::numeric_limits<Exponent>::max)()) || (e < (std::numeric_limits<Exponent>::min)()))
  1748. {
  1749. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::infinity().backend();
  1750. if (e < 0)
  1751. res.negate();
  1752. }
  1753. else
  1754. eval_ldexp(res, arg, static_cast<Exponent>(e));
  1755. }
  1756. /*
  1757. * Sign manipulation
  1758. */
  1759. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1760. unsigned Digits2, digit_base_type DigitBase2, class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2>
  1761. inline void eval_abs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits2, DigitBase2, Allocator2, Exponent2, MinE2, MaxE2>& arg)
  1762. {
  1763. res = arg;
  1764. res.sign() = false;
  1765. }
  1766. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1767. inline void eval_abs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1768. {
  1769. res = arg;
  1770. res.sign() = false;
  1771. }
  1772. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE,
  1773. unsigned Digits2, digit_base_type DigitBase2, class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2>
  1774. inline void eval_fabs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits2, DigitBase2, Allocator2, Exponent2, MinE2, MaxE2>& arg)
  1775. {
  1776. res = arg;
  1777. res.sign() = false;
  1778. }
  1779. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1780. inline void eval_fabs(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1781. {
  1782. res = arg;
  1783. res.sign() = false;
  1784. }
  1785. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1786. inline int eval_fpclassify(const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1787. {
  1788. switch (arg.exponent())
  1789. {
  1790. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1791. return FP_ZERO;
  1792. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1793. return FP_INFINITE;
  1794. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1795. return FP_NAN;
  1796. default:
  1797. break;
  1798. }
  1799. return FP_NORMAL;
  1800. }
  1801. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1802. inline void eval_sqrt(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1803. {
  1804. using default_ops::eval_bit_test;
  1805. using default_ops::eval_increment;
  1806. using default_ops::eval_integer_sqrt;
  1807. switch (arg.exponent())
  1808. {
  1809. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1810. errno = EDOM;
  1811. BOOST_MP_FALLTHROUGH;
  1812. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1813. res = arg;
  1814. return;
  1815. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1816. if (arg.sign())
  1817. {
  1818. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1819. errno = EDOM;
  1820. }
  1821. else
  1822. res = arg;
  1823. return;
  1824. default:
  1825. break;
  1826. }
  1827. if (arg.sign())
  1828. {
  1829. res = std::numeric_limits<number<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > >::quiet_NaN().backend();
  1830. errno = EDOM;
  1831. return;
  1832. }
  1833. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::double_rep_type t(arg.bits()), r, s;
  1834. eval_left_shift(t, arg.exponent() & 1 ? cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count : cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1);
  1835. eval_integer_sqrt(s, r, t);
  1836. if (!eval_bit_test(s, cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count))
  1837. {
  1838. // We have exactly the right number of cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count in the result, round as required:
  1839. if (s.compare(r) < 0)
  1840. {
  1841. eval_increment(s);
  1842. }
  1843. }
  1844. typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type ae = arg.exponent();
  1845. res.exponent() = ae / 2;
  1846. res.sign() = false;
  1847. if ((ae & 1) && (ae < 0))
  1848. --res.exponent();
  1849. copy_and_round(res, s);
  1850. }
  1851. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1852. inline void eval_floor(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1853. {
  1854. using default_ops::eval_increment;
  1855. switch (arg.exponent())
  1856. {
  1857. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1858. errno = EDOM;
  1859. BOOST_MP_FALLTHROUGH;
  1860. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1861. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1862. res = arg;
  1863. return;
  1864. default:
  1865. break;
  1866. }
  1867. using shift_type = typename std::conditional<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type;
  1868. shift_type shift =
  1869. (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1;
  1870. if ((arg.exponent() > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0))
  1871. {
  1872. // Either arg is already an integer, or a special value:
  1873. res = arg;
  1874. return;
  1875. }
  1876. if (shift >= (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
  1877. {
  1878. res = static_cast<signed_limb_type>(arg.sign() ? -1 : 0);
  1879. return;
  1880. }
  1881. bool fractional = (shift_type)eval_lsb(arg.bits()) < shift;
  1882. res = arg;
  1883. eval_right_shift(res.bits(), static_cast<double_limb_type>(shift));
  1884. if (fractional && res.sign())
  1885. {
  1886. eval_increment(res.bits());
  1887. const std::ptrdiff_t shift_check =
  1888. static_cast<std::ptrdiff_t>(static_cast<std::ptrdiff_t>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1 - static_cast<std::ptrdiff_t>(shift));
  1889. if (static_cast<std::ptrdiff_t>(eval_msb(res.bits())) != shift_check)
  1890. {
  1891. // Must have extended result by one bit in the increment:
  1892. --shift;
  1893. ++res.exponent();
  1894. }
  1895. }
  1896. eval_left_shift(res.bits(), static_cast<double_limb_type>(shift));
  1897. }
  1898. template <unsigned Digits, digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE>
  1899. inline void eval_ceil(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& res, const cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>& arg)
  1900. {
  1901. using default_ops::eval_increment;
  1902. switch (arg.exponent())
  1903. {
  1904. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity:
  1905. errno = EDOM;
  1906. BOOST_MP_FALLTHROUGH;
  1907. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_zero:
  1908. case cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan:
  1909. res = arg;
  1910. return;
  1911. default:
  1912. break;
  1913. }
  1914. using shift_type = typename std::conditional<sizeof(typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type) < sizeof(int), int, typename cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type>::type;
  1915. shift_type shift = (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - arg.exponent() - 1;
  1916. if ((arg.exponent() > (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent) || (shift <= 0))
  1917. {
  1918. // Either arg is already an integer, or a special value:
  1919. res = arg;
  1920. return;
  1921. }
  1922. if (shift >= (shift_type)cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count)
  1923. {
  1924. bool s = arg.sign(); // takes care of signed zeros
  1925. res = static_cast<signed_limb_type>(arg.sign() ? 0 : 1);
  1926. res.sign() = s;
  1927. return;
  1928. }
  1929. bool fractional = (shift_type)eval_lsb(arg.bits()) < shift;
  1930. res = arg;
  1931. eval_right_shift(res.bits(), static_cast<unsigned>(shift));
  1932. if (fractional && !res.sign())
  1933. {
  1934. eval_increment(res.bits());
  1935. if ((std::ptrdiff_t)eval_msb(res.bits()) != (std::ptrdiff_t)(static_cast<shift_type>(cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count) - 1 - shift))
  1936. {
  1937. // Must have extended result by one bit in the increment:
  1938. --shift;
  1939. ++res.exponent();
  1940. }
  1941. }
  1942. eval_left_shift(res.bits(), static_cast<unsigned>(shift));
  1943. }
  1944. template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2>
  1945. int eval_signbit(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val)
  1946. {
  1947. return val.sign();
  1948. }
  1949. template <unsigned D1, backends::digit_base_type B1, class A1, class E1, E1 M1, E1 M2>
  1950. inline std::size_t hash_value(const cpp_bin_float<D1, B1, A1, E1, M1, M2>& val)
  1951. {
  1952. std::size_t result = hash_value(val.bits());
  1953. boost::multiprecision::detail::hash_combine(result, val.exponent(), val.sign());
  1954. return result;
  1955. }
  1956. } // namespace backends
  1957. namespace detail {
  1958. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinExponent, Exponent MaxExponent>
  1959. struct transcendental_reduction_type<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent> >
  1960. {
  1961. //
  1962. // The type used for trigonometric reduction needs 3 times the precision of the base type.
  1963. // This is double the precision of the original type, plus the largest exponent supported.
  1964. // As a practical measure the largest argument supported is 1/eps, as supporting larger
  1965. // arguments requires the division of argument by PI/2 to also be done at higher precision,
  1966. // otherwise the result (an integer) can not be represented exactly.
  1967. //
  1968. // See ARGUMENT REDUCTION FOR HUGE ARGUMENTS. K C Ng.
  1969. //
  1970. using type = boost::multiprecision::backends::cpp_bin_float<
  1971. boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent>::bit_count * 3,
  1972. boost::multiprecision::backends::digit_base_2,
  1973. Allocator, Exponent, MinExponent, MaxExponent>;
  1974. };
  1975. #ifdef BOOST_HAS_INT128
  1976. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinExponent, Exponent MaxExponent>
  1977. struct is_convertible_arithmetic<int128_type, boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent> > : public std::true_type
  1978. {};
  1979. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinExponent, Exponent MaxExponent>
  1980. struct is_convertible_arithmetic<uint128_type, boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent> > : public std::true_type
  1981. {};
  1982. #endif
  1983. } // namespace detail
  1984. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator, boost::multiprecision::expression_template_option ExpressionTemplates>
  1985. inline boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>
  1986. copysign BOOST_PREVENT_MACRO_SUBSTITUTION(
  1987. const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& a,
  1988. const boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>& b)
  1989. {
  1990. boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> res(a);
  1991. res.backend().sign() = b.backend().sign();
  1992. return res;
  1993. }
  1994. template <unsigned Digits, backends::digit_base_type DigitBase, class Exponent, Exponent MinE, Exponent MaxE, class Allocator>
  1995. struct number_category<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE> > : public std::integral_constant<int, boost::multiprecision::number_kind_floating_point>
  1996. {};
  1997. template <unsigned Digits, backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, class Allocator2, class Exponent2, Exponent MinE2, Exponent MaxE2>
  1998. struct is_equivalent_number_type<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, cpp_bin_float<Digits, DigitBase, Allocator2, Exponent2, MinE2, MaxE2> >
  1999. : public std::integral_constant<bool, true> {};
  2000. } // namespace multiprecision
  2001. namespace math {
  2002. using boost::multiprecision::copysign;
  2003. using boost::multiprecision::signbit;
  2004. } // namespace math
  2005. } // namespace boost
  2006. #include <boost/multiprecision/cpp_bin_float/io.hpp>
  2007. #include <boost/multiprecision/cpp_bin_float/transcendental.hpp>
  2008. namespace std {
  2009. //
  2010. // numeric_limits [partial] specializations for the types declared in this header:
  2011. //
  2012. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2013. class numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >
  2014. {
  2015. using number_type = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates>;
  2016. private:
  2017. //
  2018. // Functions to calculate cached values stored in static values:
  2019. //
  2020. static number_type get_min()
  2021. {
  2022. using ui_type = typename std::tuple_element<0, typename number_type::backend_type::unsigned_types>::type;
  2023. number_type value(ui_type(1u));
  2024. value.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
  2025. return value;
  2026. }
  2027. #ifdef BOOST_MSVC
  2028. #pragma warning(push)
  2029. #pragma warning(disable : 4127) // conditional expression is constant
  2030. #endif
  2031. static number_type get_max()
  2032. {
  2033. number_type value;
  2034. BOOST_IF_CONSTEXPR(std::is_void<Allocator>::value)
  2035. eval_complement(value.backend().bits(), value.backend().bits());
  2036. else
  2037. {
  2038. // We jump through hoops here using the backend type directly just to keep VC12 happy
  2039. // (ie compiler workaround, for very strange compiler bug):
  2040. using boost::multiprecision::default_ops::eval_add;
  2041. using boost::multiprecision::default_ops::eval_decrement;
  2042. using boost::multiprecision::default_ops::eval_left_shift;
  2043. using int_backend_type = typename number_type::backend_type::rep_type;
  2044. using ui_type = typename std::tuple_element<0, typename int_backend_type::unsigned_types>::type;
  2045. int_backend_type i;
  2046. i = ui_type(1u);
  2047. eval_left_shift(i, boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count - 1);
  2048. int_backend_type j(i);
  2049. eval_decrement(i);
  2050. eval_add(j, i);
  2051. value.backend().bits() = j;
  2052. }
  2053. value.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
  2054. return value;
  2055. }
  2056. #ifdef BOOST_MSVC
  2057. #pragma warning(pop)
  2058. #endif
  2059. static number_type get_epsilon()
  2060. {
  2061. using ui_type = typename std::tuple_element<0, typename number_type::backend_type::unsigned_types>::type;
  2062. number_type value(ui_type(1u));
  2063. return ldexp(value, 1 - static_cast<int>(digits));
  2064. }
  2065. // What value should this be????
  2066. static number_type get_round_error()
  2067. {
  2068. // returns 0.5
  2069. return ldexp(number_type(1u), -1);
  2070. }
  2071. static number_type get_infinity()
  2072. {
  2073. number_type value;
  2074. value.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_infinity;
  2075. return value;
  2076. }
  2077. static number_type get_quiet_NaN()
  2078. {
  2079. number_type value;
  2080. value.backend().exponent() = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_nan;
  2081. return value;
  2082. }
  2083. public:
  2084. static constexpr bool is_specialized = true;
  2085. static number_type(min)()
  2086. {
  2087. // C++11 thread safe static initialization:
  2088. static number_type value = get_min();
  2089. return value;
  2090. }
  2091. static number_type(max)()
  2092. {
  2093. // C++11 thread safe static initialization:
  2094. static number_type value = get_max();
  2095. return value;
  2096. }
  2097. static constexpr number_type lowest()
  2098. {
  2099. return -(max)();
  2100. }
  2101. static constexpr int digits = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::bit_count;
  2102. static constexpr int digits10 = boost::multiprecision::detail::calc_digits10<static_cast<unsigned>(digits)>::value;
  2103. // Is this really correct???
  2104. static constexpr int max_digits10 = boost::multiprecision::detail::calc_max_digits10<static_cast<unsigned>(digits)>::value;
  2105. static constexpr bool is_signed = true;
  2106. static constexpr bool is_integer = false;
  2107. static constexpr bool is_exact = false;
  2108. static constexpr int radix = 2;
  2109. static number_type epsilon()
  2110. {
  2111. // C++11 thread safe static initialization:
  2112. static number_type value = get_epsilon();
  2113. return value;
  2114. }
  2115. // What value should this be????
  2116. static number_type round_error()
  2117. {
  2118. // returns 0.5
  2119. // C++11 thread safe static initialization:
  2120. static number_type value = get_round_error();
  2121. return value;
  2122. }
  2123. static constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::min_exponent;
  2124. static constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type min_exponent10 = (min_exponent / 1000) * 301L;
  2125. static constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent = boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::max_exponent;
  2126. static constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type max_exponent10 = (max_exponent / 1000) * 301L;
  2127. static constexpr bool has_infinity = true;
  2128. static constexpr bool has_quiet_NaN = true;
  2129. static constexpr bool has_signaling_NaN = false;
  2130. #ifdef _MSC_VER
  2131. #pragma warning(push)
  2132. #pragma warning(disable:4996)
  2133. #endif
  2134. static constexpr float_denorm_style has_denorm = denorm_absent;
  2135. #ifdef _MSC_VER
  2136. #pragma warning(pop)
  2137. #endif
  2138. static constexpr bool has_denorm_loss = false;
  2139. static number_type infinity()
  2140. {
  2141. // C++11 thread safe static initialization:
  2142. static number_type value = get_infinity();
  2143. return value;
  2144. }
  2145. static number_type quiet_NaN()
  2146. {
  2147. // C++11 thread safe static initialization:
  2148. static number_type value = get_quiet_NaN();
  2149. return value;
  2150. }
  2151. static constexpr number_type signaling_NaN()
  2152. {
  2153. return number_type(0);
  2154. }
  2155. static constexpr number_type denorm_min() { return get_min(); }
  2156. static constexpr bool is_iec559 = false;
  2157. static constexpr bool is_bounded = true;
  2158. static constexpr bool is_modulo = false;
  2159. static constexpr bool traps = true;
  2160. static constexpr bool tinyness_before = false;
  2161. static constexpr float_round_style round_style = round_to_nearest;
  2162. };
  2163. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2164. constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits;
  2165. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2166. constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::digits10;
  2167. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2168. constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_digits10;
  2169. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2170. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_signed;
  2171. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2172. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_integer;
  2173. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2174. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_exact;
  2175. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2176. constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::radix;
  2177. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2178. constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent;
  2179. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2180. constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::min_exponent10;
  2181. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2182. constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent;
  2183. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2184. constexpr typename boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>::exponent_type numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::max_exponent10;
  2185. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2186. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_infinity;
  2187. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2188. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_quiet_NaN;
  2189. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2190. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_signaling_NaN;
  2191. #ifdef _MSC_VER
  2192. #pragma warning(push)
  2193. #pragma warning(disable:4996)
  2194. #endif
  2195. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2196. constexpr float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm;
  2197. #ifdef _MSC_VER
  2198. #pragma warning(pop)
  2199. #endif
  2200. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2201. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::has_denorm_loss;
  2202. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2203. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_iec559;
  2204. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2205. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_bounded;
  2206. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2207. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::is_modulo;
  2208. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2209. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::traps;
  2210. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2211. constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::tinyness_before;
  2212. template <unsigned Digits, boost::multiprecision::backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinE, Exponent MaxE, boost::multiprecision::expression_template_option ExpressionTemplates>
  2213. constexpr float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinE, MaxE>, ExpressionTemplates> >::round_style;
  2214. } // namespace std
  2215. #endif