cpp_bin_float.hpp 91 KB

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