cstdfloat_cmath.hpp 37 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // Copyright Christopher Kormanyos 2014.
  3. // Copyright John Maddock 2014.
  4. // Copyright Paul Bristow 2014.
  5. // Distributed under the Boost Software License,
  6. // Version 1.0. (See accompanying file LICENSE_1_0.txt
  7. // or copy at http://www.boost.org/LICENSE_1_0.txt)
  8. //
  9. // Implement quadruple-precision <cmath> support.
  10. #ifndef _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_
  11. #define _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_
  12. #include <boost/math/cstdfloat/cstdfloat_types.hpp>
  13. #include <boost/math/cstdfloat/cstdfloat_limits.hpp>
  14. #if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT)
  15. #include <cmath>
  16. #include <stdexcept>
  17. #include <boost/cstdint.hpp>
  18. #include <boost/static_assert.hpp>
  19. #include <boost/throw_exception.hpp>
  20. #if defined(_WIN32) && defined(__GNUC__)
  21. // Several versions of Mingw and probably cygwin too have broken
  22. // libquadmath implementations that segfault as soon as you call
  23. // expq or any function that depends on it.
  24. #define BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS
  25. #endif
  26. // Here is a helper function used for raising the value of a given
  27. // floating-point type to the power of n, where n has integral type.
  28. namespace boost { namespace math { namespace cstdfloat { namespace detail {
  29. template<class float_type, class integer_type>
  30. inline float_type pown(const float_type& x, const integer_type p)
  31. {
  32. const bool isneg = (x < 0);
  33. const bool isnan = (x != x);
  34. const bool isinf = ((!isneg) ? bool(+x > (std::numeric_limits<float_type>::max)())
  35. : bool(-x > (std::numeric_limits<float_type>::max)()));
  36. if(isnan) { return x; }
  37. if(isinf) { return std::numeric_limits<float_type>::quiet_NaN(); }
  38. const bool x_is_neg = (x < 0);
  39. const float_type abs_x = (x_is_neg ? -x : x);
  40. if(p < static_cast<integer_type>(0))
  41. {
  42. if(abs_x < (std::numeric_limits<float_type>::min)())
  43. {
  44. return (x_is_neg ? -std::numeric_limits<float_type>::infinity()
  45. : +std::numeric_limits<float_type>::infinity());
  46. }
  47. else
  48. {
  49. return float_type(1) / pown(x, static_cast<integer_type>(-p));
  50. }
  51. }
  52. if(p == static_cast<integer_type>(0))
  53. {
  54. return float_type(1);
  55. }
  56. else
  57. {
  58. if(p == static_cast<integer_type>(1)) { return x; }
  59. if(abs_x > (std::numeric_limits<float_type>::max)())
  60. {
  61. return (x_is_neg ? -std::numeric_limits<float_type>::infinity()
  62. : +std::numeric_limits<float_type>::infinity());
  63. }
  64. if (p == static_cast<integer_type>(2)) { return (x * x); }
  65. else if(p == static_cast<integer_type>(3)) { return ((x * x) * x); }
  66. else if(p == static_cast<integer_type>(4)) { const float_type x2 = (x * x); return (x2 * x2); }
  67. else
  68. {
  69. // The variable xn stores the binary powers of x.
  70. float_type result(((p % integer_type(2)) != integer_type(0)) ? x : float_type(1));
  71. float_type xn (x);
  72. integer_type p2 = p;
  73. while(integer_type(p2 /= 2) != integer_type(0))
  74. {
  75. // Square xn for each binary power.
  76. xn *= xn;
  77. const bool has_binary_power = (integer_type(p2 % integer_type(2)) != integer_type(0));
  78. if(has_binary_power)
  79. {
  80. // Multiply the result with each binary power contained in the exponent.
  81. result *= xn;
  82. }
  83. }
  84. return result;
  85. }
  86. }
  87. }
  88. } } } } // boost::math::cstdfloat::detail
  89. // We will now define preprocessor symbols representing quadruple-precision <cmath> functions.
  90. #if defined(BOOST_INTEL)
  91. #define BOOST_CSTDFLOAT_FLOAT128_LDEXP __ldexpq
  92. #define BOOST_CSTDFLOAT_FLOAT128_FREXP __frexpq
  93. #define BOOST_CSTDFLOAT_FLOAT128_FABS __fabsq
  94. #define BOOST_CSTDFLOAT_FLOAT128_FLOOR __floorq
  95. #define BOOST_CSTDFLOAT_FLOAT128_CEIL __ceilq
  96. #if !defined(BOOST_CSTDFLOAT_FLOAT128_SQRT)
  97. #define BOOST_CSTDFLOAT_FLOAT128_SQRT __sqrtq
  98. #endif
  99. #define BOOST_CSTDFLOAT_FLOAT128_TRUNC __truncq
  100. #define BOOST_CSTDFLOAT_FLOAT128_EXP __expq
  101. #define BOOST_CSTDFLOAT_FLOAT128_EXPM1 __expm1q
  102. #define BOOST_CSTDFLOAT_FLOAT128_POW __powq
  103. #define BOOST_CSTDFLOAT_FLOAT128_LOG __logq
  104. #define BOOST_CSTDFLOAT_FLOAT128_LOG10 __log10q
  105. #define BOOST_CSTDFLOAT_FLOAT128_SIN __sinq
  106. #define BOOST_CSTDFLOAT_FLOAT128_COS __cosq
  107. #define BOOST_CSTDFLOAT_FLOAT128_TAN __tanq
  108. #define BOOST_CSTDFLOAT_FLOAT128_ASIN __asinq
  109. #define BOOST_CSTDFLOAT_FLOAT128_ACOS __acosq
  110. #define BOOST_CSTDFLOAT_FLOAT128_ATAN __atanq
  111. #define BOOST_CSTDFLOAT_FLOAT128_SINH __sinhq
  112. #define BOOST_CSTDFLOAT_FLOAT128_COSH __coshq
  113. #define BOOST_CSTDFLOAT_FLOAT128_TANH __tanhq
  114. #define BOOST_CSTDFLOAT_FLOAT128_ASINH __asinhq
  115. #define BOOST_CSTDFLOAT_FLOAT128_ACOSH __acoshq
  116. #define BOOST_CSTDFLOAT_FLOAT128_ATANH __atanhq
  117. #define BOOST_CSTDFLOAT_FLOAT128_FMOD __fmodq
  118. #define BOOST_CSTDFLOAT_FLOAT128_ATAN2 __atan2q
  119. #define BOOST_CSTDFLOAT_FLOAT128_LGAMMA __lgammaq
  120. #define BOOST_CSTDFLOAT_FLOAT128_TGAMMA __tgammaq
  121. #elif defined(__GNUC__)
  122. #define BOOST_CSTDFLOAT_FLOAT128_LDEXP ldexpq
  123. #define BOOST_CSTDFLOAT_FLOAT128_FREXP frexpq
  124. #define BOOST_CSTDFLOAT_FLOAT128_FABS fabsq
  125. #define BOOST_CSTDFLOAT_FLOAT128_FLOOR floorq
  126. #define BOOST_CSTDFLOAT_FLOAT128_CEIL ceilq
  127. #if !defined(BOOST_CSTDFLOAT_FLOAT128_SQRT)
  128. #define BOOST_CSTDFLOAT_FLOAT128_SQRT sqrtq
  129. #endif
  130. #define BOOST_CSTDFLOAT_FLOAT128_TRUNC truncq
  131. #define BOOST_CSTDFLOAT_FLOAT128_POW powq
  132. #define BOOST_CSTDFLOAT_FLOAT128_LOG logq
  133. #define BOOST_CSTDFLOAT_FLOAT128_LOG10 log10q
  134. #define BOOST_CSTDFLOAT_FLOAT128_SIN sinq
  135. #define BOOST_CSTDFLOAT_FLOAT128_COS cosq
  136. #define BOOST_CSTDFLOAT_FLOAT128_TAN tanq
  137. #define BOOST_CSTDFLOAT_FLOAT128_ASIN asinq
  138. #define BOOST_CSTDFLOAT_FLOAT128_ACOS acosq
  139. #define BOOST_CSTDFLOAT_FLOAT128_ATAN atanq
  140. #define BOOST_CSTDFLOAT_FLOAT128_FMOD fmodq
  141. #define BOOST_CSTDFLOAT_FLOAT128_ATAN2 atan2q
  142. #define BOOST_CSTDFLOAT_FLOAT128_LGAMMA lgammaq
  143. #if !defined(BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS)
  144. #define BOOST_CSTDFLOAT_FLOAT128_EXP expq
  145. #define BOOST_CSTDFLOAT_FLOAT128_EXPM1 expm1q_internal
  146. #define BOOST_CSTDFLOAT_FLOAT128_SINH sinhq
  147. #define BOOST_CSTDFLOAT_FLOAT128_COSH coshq
  148. #define BOOST_CSTDFLOAT_FLOAT128_TANH tanhq
  149. #define BOOST_CSTDFLOAT_FLOAT128_ASINH asinhq
  150. #define BOOST_CSTDFLOAT_FLOAT128_ACOSH acoshq
  151. #define BOOST_CSTDFLOAT_FLOAT128_ATANH atanhq
  152. #define BOOST_CSTDFLOAT_FLOAT128_TGAMMA tgammaq
  153. #else // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS
  154. #define BOOST_CSTDFLOAT_FLOAT128_EXP expq_patch
  155. #define BOOST_CSTDFLOAT_FLOAT128_SINH sinhq_patch
  156. #define BOOST_CSTDFLOAT_FLOAT128_COSH coshq_patch
  157. #define BOOST_CSTDFLOAT_FLOAT128_TANH tanhq_patch
  158. #define BOOST_CSTDFLOAT_FLOAT128_ASINH asinhq_patch
  159. #define BOOST_CSTDFLOAT_FLOAT128_ACOSH acoshq_patch
  160. #define BOOST_CSTDFLOAT_FLOAT128_ATANH atanhq_patch
  161. #define BOOST_CSTDFLOAT_FLOAT128_TGAMMA tgammaq_patch
  162. #endif // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS
  163. #endif
  164. // Implement quadruple-precision <cmath> functions in the namespace
  165. // boost::math::cstdfloat::detail. Subsequently inject these into the
  166. // std namespace via *using* directive.
  167. // Begin with some forward function declarations. Also implement patches
  168. // for compilers that have broken float128 exponential functions.
  169. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LDEXP (boost::math::cstdfloat::detail::float_internal128_t, int) throw();
  170. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FREXP (boost::math::cstdfloat::detail::float_internal128_t, int*) throw();
  171. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FABS (boost::math::cstdfloat::detail::float_internal128_t) throw();
  172. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FLOOR (boost::math::cstdfloat::detail::float_internal128_t) throw();
  173. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_CEIL (boost::math::cstdfloat::detail::float_internal128_t) throw();
  174. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SQRT (boost::math::cstdfloat::detail::float_internal128_t) throw();
  175. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TRUNC (boost::math::cstdfloat::detail::float_internal128_t) throw();
  176. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_POW (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw();
  177. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG (boost::math::cstdfloat::detail::float_internal128_t) throw();
  178. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG10 (boost::math::cstdfloat::detail::float_internal128_t) throw();
  179. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SIN (boost::math::cstdfloat::detail::float_internal128_t) throw();
  180. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COS (boost::math::cstdfloat::detail::float_internal128_t) throw();
  181. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TAN (boost::math::cstdfloat::detail::float_internal128_t) throw();
  182. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASIN (boost::math::cstdfloat::detail::float_internal128_t) throw();
  183. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOS (boost::math::cstdfloat::detail::float_internal128_t) throw();
  184. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN (boost::math::cstdfloat::detail::float_internal128_t) throw();
  185. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FMOD (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw();
  186. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN2 (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw();
  187. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LGAMMA(boost::math::cstdfloat::detail::float_internal128_t) throw();
  188. #if !defined(BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS)
  189. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::math::cstdfloat::detail::float_internal128_t x) throw();
  190. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH (boost::math::cstdfloat::detail::float_internal128_t x) throw();
  191. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH (boost::math::cstdfloat::detail::float_internal128_t x) throw();
  192. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH (boost::math::cstdfloat::detail::float_internal128_t x) throw();
  193. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASINH (boost::math::cstdfloat::detail::float_internal128_t x) throw();
  194. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOSH (boost::math::cstdfloat::detail::float_internal128_t x) throw();
  195. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATANH (boost::math::cstdfloat::detail::float_internal128_t x) throw();
  196. extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::math::cstdfloat::detail::float_internal128_t x) throw();
  197. #else // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS
  198. // Forward declaration of the patched exponent function, exp(x).
  199. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::math::cstdfloat::detail::float_internal128_t x);
  200. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXPM1 (boost::math::cstdfloat::detail::float_internal128_t x)
  201. {
  202. // Compute exp(x) - 1 for x small.
  203. // Use an order-36 polynomial approximation of the exponential function
  204. // in the range of (-ln2 < x < ln2). Scale the argument to this range
  205. // and subsequently multiply the result by 2^n accordingly.
  206. // Derive the polynomial coefficients with Mathematica(R) by generating
  207. // a table of high-precision values of exp(x) in the range (-ln2 < x < ln2)
  208. // and subsequently applying the built-in *Fit* function.
  209. // Table[{x, Exp[x] - 1}, {x, -Log[2], Log[2], 1/180}]
  210. // N[%, 120]
  211. // Fit[%, {x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10, x^11, x^12,
  212. // x^13, x^14, x^15, x^16, x^17, x^18, x^19, x^20, x^21, x^22,
  213. // x^23, x^24, x^25, x^26, x^27, x^28, x^29, x^30, x^31, x^32,
  214. // x^33, x^34, x^35, x^36}, x]
  215. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  216. float_type sum;
  217. if(x > BOOST_FLOAT128_C(0.693147180559945309417232121458176568075500134360255))
  218. {
  219. sum = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x) - float_type(1);
  220. }
  221. else
  222. {
  223. // Compute the polynomial approximation of exp(alpha).
  224. sum = (((((((((((((((((((((((((((((((((((( float_type(BOOST_FLOAT128_C(2.69291698127774166063293705964720493864630783729857438187365E-42)) * x
  225. + float_type(BOOST_FLOAT128_C(9.70937085471487654794114679403710456028986572118859594614033E-41))) * x
  226. + float_type(BOOST_FLOAT128_C(3.38715585158055097155585505318085512156885389014410753080500E-39))) * x
  227. + float_type(BOOST_FLOAT128_C(1.15162718532861050809222658798662695267019717760563645440433E-37))) * x
  228. + float_type(BOOST_FLOAT128_C(3.80039074689434663295873584133017767349635602413675471702393E-36))) * x
  229. + float_type(BOOST_FLOAT128_C(1.21612504934087520075905434734158045947460467096773246215239E-34))) * x
  230. + float_type(BOOST_FLOAT128_C(3.76998762883139753126119821241037824830069851253295480396224E-33))) * x
  231. + float_type(BOOST_FLOAT128_C(1.13099628863830344684998293828608215735777107850991029729440E-31))) * x
  232. + float_type(BOOST_FLOAT128_C(3.27988923706982293204067897468714277771890104022419696770352E-30))) * x
  233. + float_type(BOOST_FLOAT128_C(9.18368986379558482800593745627556950089950023355628325088207E-29))) * x
  234. + float_type(BOOST_FLOAT128_C(2.47959626322479746949155352659617642905315302382639380521497E-27))) * x
  235. + float_type(BOOST_FLOAT128_C(6.44695028438447337900255966737803112935639344283098705091949E-26))) * x
  236. + float_type(BOOST_FLOAT128_C(1.61173757109611834904452725462599961406036904573072897122957E-24))) * x
  237. + float_type(BOOST_FLOAT128_C(3.86817017063068403772269360016918092488847584660382953555804E-23))) * x
  238. + float_type(BOOST_FLOAT128_C(8.89679139245057328674891109315654704307721758924206107351744E-22))) * x
  239. + float_type(BOOST_FLOAT128_C(1.95729410633912612308475595397946731738088422488032228717097E-20))) * x
  240. + float_type(BOOST_FLOAT128_C(4.11031762331216485847799061511674191805055663711439605760231E-19))) * x
  241. + float_type(BOOST_FLOAT128_C(8.22063524662432971695598123977873600603370758794431071426640E-18))) * x
  242. + float_type(BOOST_FLOAT128_C(1.56192069685862264622163643500633782667263448653185159383285E-16))) * x
  243. + float_type(BOOST_FLOAT128_C(2.81145725434552076319894558300988749849555291507956994126835E-15))) * x
  244. + float_type(BOOST_FLOAT128_C(4.77947733238738529743820749111754320727153728139716409114011E-14))) * x
  245. + float_type(BOOST_FLOAT128_C(7.64716373181981647590113198578807092707697416852226691068627E-13))) * x
  246. + float_type(BOOST_FLOAT128_C(1.14707455977297247138516979786821056670509688396295740818677E-11))) * x
  247. + float_type(BOOST_FLOAT128_C(1.60590438368216145993923771701549479323291461578567184216302E-10))) * x
  248. + float_type(BOOST_FLOAT128_C(2.08767569878680989792100903212014323125428376052986408239620E-09))) * x
  249. + float_type(BOOST_FLOAT128_C(2.50521083854417187750521083854417187750523408006206780016659E-08))) * x
  250. + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573195144226062684604E-07))) * x
  251. + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573191310049321957902E-06))) * x
  252. + float_type(BOOST_FLOAT128_C(0.00002480158730158730158730158730158730158730158730149317774))) * x
  253. + float_type(BOOST_FLOAT128_C(0.00019841269841269841269841269841269841269841269841293575920))) * x
  254. + float_type(BOOST_FLOAT128_C(0.00138888888888888888888888888888888888888888888888889071045))) * x
  255. + float_type(BOOST_FLOAT128_C(0.00833333333333333333333333333333333333333333333333332986595))) * x
  256. + float_type(BOOST_FLOAT128_C(0.04166666666666666666666666666666666666666666666666666664876))) * x
  257. + float_type(BOOST_FLOAT128_C(0.16666666666666666666666666666666666666666666666666666669048))) * x
  258. + float_type(BOOST_FLOAT128_C(0.50000000000000000000000000000000000000000000000000000000006))) * x
  259. + float_type(BOOST_FLOAT128_C(0.99999999999999999999999999999999999999999999999999999999995))) * x);
  260. }
  261. return sum;
  262. }
  263. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::math::cstdfloat::detail::float_internal128_t x)
  264. {
  265. // Patch the expq() function for a subset of broken GCC compilers
  266. // like GCC 4.7, 4.8 on MinGW.
  267. // Use an order-36 polynomial approximation of the exponential function
  268. // in the range of (-ln2 < x < ln2). Scale the argument to this range
  269. // and subsequently multiply the result by 2^n accordingly.
  270. // Derive the polynomial coefficients with Mathematica(R) by generating
  271. // a table of high-precision values of exp(x) in the range (-ln2 < x < ln2)
  272. // and subsequently applying the built-in *Fit* function.
  273. // Table[{x, Exp[x] - 1}, {x, -Log[2], Log[2], 1/180}]
  274. // N[%, 120]
  275. // Fit[%, {x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10, x^11, x^12,
  276. // x^13, x^14, x^15, x^16, x^17, x^18, x^19, x^20, x^21, x^22,
  277. // x^23, x^24, x^25, x^26, x^27, x^28, x^29, x^30, x^31, x^32,
  278. // x^33, x^34, x^35, x^36}, x]
  279. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  280. // Scale the argument x to the range (-ln2 < x < ln2).
  281. BOOST_CONSTEXPR_OR_CONST float_type one_over_ln2 = float_type(BOOST_FLOAT128_C(1.44269504088896340735992468100189213742664595415299));
  282. const float_type x_over_ln2 = x * one_over_ln2;
  283. boost::int_fast32_t n;
  284. if(x != x)
  285. {
  286. // The argument is NaN.
  287. return std::numeric_limits<float_type>::quiet_NaN();
  288. }
  289. else if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x) > BOOST_FLOAT128_C(+0.693147180559945309417232121458176568075500134360255))
  290. {
  291. // The absolute value of the argument exceeds ln2.
  292. n = static_cast<boost::int_fast32_t>(::BOOST_CSTDFLOAT_FLOAT128_FLOOR(x_over_ln2));
  293. }
  294. else if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x) < BOOST_FLOAT128_C(+0.693147180559945309417232121458176568075500134360255))
  295. {
  296. // The absolute value of the argument is less than ln2.
  297. n = static_cast<boost::int_fast32_t>(0);
  298. }
  299. else
  300. {
  301. // The absolute value of the argument is exactly equal to ln2 (in the sense of floating-point equality).
  302. return float_type(2);
  303. }
  304. // Check if the argument is very near an integer.
  305. const float_type floor_of_x = ::BOOST_CSTDFLOAT_FLOAT128_FLOOR(x);
  306. if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x - floor_of_x) < float_type(BOOST_CSTDFLOAT_FLOAT128_EPS))
  307. {
  308. // Return e^n for arguments very near an integer.
  309. return boost::math::cstdfloat::detail::pown(BOOST_FLOAT128_C(2.71828182845904523536028747135266249775724709369996), static_cast<boost::int_fast32_t>(floor_of_x));
  310. }
  311. // Compute the scaled argument alpha.
  312. const float_type alpha = x - (n * BOOST_FLOAT128_C(0.693147180559945309417232121458176568075500134360255));
  313. // Compute the polynomial approximation of expm1(alpha) and add to it
  314. // in order to obtain the scaled result.
  315. const float_type scaled_result = ::BOOST_CSTDFLOAT_FLOAT128_EXPM1(alpha) + float_type(1);
  316. // Rescale the result and return it.
  317. return scaled_result * boost::math::cstdfloat::detail::pown(float_type(2), n);
  318. }
  319. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH (boost::math::cstdfloat::detail::float_internal128_t x)
  320. {
  321. // Patch the sinhq() function for a subset of broken GCC compilers
  322. // like GCC 4.7, 4.8 on MinGW.
  323. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  324. // Here, we use the following:
  325. // Set: ex = exp(x)
  326. // Set: em1 = expm1(x)
  327. // Then
  328. // sinh(x) = (ex - 1/ex) / 2 ; for |x| >= 1
  329. // sinh(x) = (2em1 + em1^2) / (2ex) ; for |x| < 1
  330. const float_type ex = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x);
  331. if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x) < float_type(+1))
  332. {
  333. const float_type em1 = ::BOOST_CSTDFLOAT_FLOAT128_EXPM1(x);
  334. return ((em1 * 2) + (em1 * em1)) / (ex * 2);
  335. }
  336. else
  337. {
  338. return (ex - (float_type(1) / ex)) / 2;
  339. }
  340. }
  341. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH (boost::math::cstdfloat::detail::float_internal128_t x)
  342. {
  343. // Patch the coshq() function for a subset of broken GCC compilers
  344. // like GCC 4.7, 4.8 on MinGW.
  345. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  346. const float_type ex = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x);
  347. return (ex + (float_type(1) / ex)) / 2;
  348. }
  349. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH (boost::math::cstdfloat::detail::float_internal128_t x)
  350. {
  351. // Patch the tanhq() function for a subset of broken GCC compilers
  352. // like GCC 4.7, 4.8 on MinGW.
  353. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  354. const float_type ex_plus = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x);
  355. const float_type ex_minus = (float_type(1) / ex_plus);
  356. return (ex_plus - ex_minus) / (ex_plus + ex_minus);
  357. }
  358. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASINH(boost::math::cstdfloat::detail::float_internal128_t x) throw()
  359. {
  360. // Patch the asinh() function since quadmath does not have it.
  361. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  362. return ::BOOST_CSTDFLOAT_FLOAT128_LOG(x + ::BOOST_CSTDFLOAT_FLOAT128_SQRT((x * x) + float_type(1)));
  363. }
  364. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOSH(boost::math::cstdfloat::detail::float_internal128_t x) throw()
  365. {
  366. // Patch the acosh() function since quadmath does not have it.
  367. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  368. const float_type zp(x + float_type(1));
  369. const float_type zm(x - float_type(1));
  370. return ::BOOST_CSTDFLOAT_FLOAT128_LOG(x + (zp * ::BOOST_CSTDFLOAT_FLOAT128_SQRT(zm / zp)));
  371. }
  372. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATANH(boost::math::cstdfloat::detail::float_internal128_t x) throw()
  373. {
  374. // Patch the atanh() function since quadmath does not have it.
  375. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  376. return ( ::BOOST_CSTDFLOAT_FLOAT128_LOG(float_type(1) + x)
  377. - ::BOOST_CSTDFLOAT_FLOAT128_LOG(float_type(1) - x)) / 2;
  378. }
  379. inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::math::cstdfloat::detail::float_internal128_t x) throw()
  380. {
  381. // Patch the tgammaq() function for a subset of broken GCC compilers
  382. // like GCC 4.7, 4.8 on MinGW.
  383. typedef boost::math::cstdfloat::detail::float_internal128_t float_type;
  384. if(x > float_type(0))
  385. {
  386. return ::BOOST_CSTDFLOAT_FLOAT128_EXP(::BOOST_CSTDFLOAT_FLOAT128_LGAMMA(x));
  387. }
  388. else if(x < float_type(0))
  389. {
  390. // For x < 0, compute tgamma(-x) and use the reflection formula.
  391. const float_type positive_x = -x;
  392. float_type gamma_value = ::BOOST_CSTDFLOAT_FLOAT128_TGAMMA(positive_x);
  393. const float_type floor_of_positive_x = ::BOOST_CSTDFLOAT_FLOAT128_FLOOR (positive_x);
  394. // Take the reflection checks (slightly adapted) from <boost/math/gamma.hpp>.
  395. const bool floor_of_z_is_equal_to_z = (positive_x == ::BOOST_CSTDFLOAT_FLOAT128_FLOOR(positive_x));
  396. BOOST_CONSTEXPR_OR_CONST float_type my_pi = BOOST_FLOAT128_C(3.14159265358979323846264338327950288419716939937511);
  397. if(floor_of_z_is_equal_to_z)
  398. {
  399. const bool is_odd = ((boost::int32_t(floor_of_positive_x) % boost::int32_t(2)) != boost::int32_t(0));
  400. return (is_odd ? -std::numeric_limits<float_type>::infinity()
  401. : +std::numeric_limits<float_type>::infinity());
  402. }
  403. const float_type sinpx_value = x * ::BOOST_CSTDFLOAT_FLOAT128_SIN(my_pi * x);
  404. gamma_value *= sinpx_value;
  405. const bool result_is_too_large_to_represent = ( (::BOOST_CSTDFLOAT_FLOAT128_FABS(gamma_value) < float_type(1))
  406. && (((std::numeric_limits<float_type>::max)() * ::BOOST_CSTDFLOAT_FLOAT128_FABS(gamma_value)) < my_pi));
  407. if(result_is_too_large_to_represent)
  408. {
  409. const bool is_odd = ((boost::int32_t(floor_of_positive_x) % boost::int32_t(2)) != boost::int32_t(0));
  410. return (is_odd ? -std::numeric_limits<float_type>::infinity()
  411. : +std::numeric_limits<float_type>::infinity());
  412. }
  413. gamma_value = -my_pi / gamma_value;
  414. if((gamma_value > float_type(0)) || (gamma_value < float_type(0)))
  415. {
  416. return gamma_value;
  417. }
  418. else
  419. {
  420. // The value of gamma is too small to represent. Return 0.0 here.
  421. return float_type(0);
  422. }
  423. }
  424. else
  425. {
  426. // Gamma of zero is complex infinity. Return NaN here.
  427. return std::numeric_limits<float_type>::quiet_NaN();
  428. }
  429. }
  430. #endif // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS
  431. // Define the quadruple-precision <cmath> functions in the namespace boost::math::cstdfloat::detail.
  432. namespace boost { namespace math { namespace cstdfloat { namespace detail {
  433. inline boost::math::cstdfloat::detail::float_internal128_t ldexp (boost::math::cstdfloat::detail::float_internal128_t x, int n) { return ::BOOST_CSTDFLOAT_FLOAT128_LDEXP (x, n); }
  434. inline boost::math::cstdfloat::detail::float_internal128_t frexp (boost::math::cstdfloat::detail::float_internal128_t x, int* pn) { return ::BOOST_CSTDFLOAT_FLOAT128_FREXP (x, pn); }
  435. inline boost::math::cstdfloat::detail::float_internal128_t fabs (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FABS (x); }
  436. inline boost::math::cstdfloat::detail::float_internal128_t abs (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FABS (x); }
  437. inline boost::math::cstdfloat::detail::float_internal128_t floor (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FLOOR (x); }
  438. inline boost::math::cstdfloat::detail::float_internal128_t ceil (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_CEIL (x); }
  439. inline boost::math::cstdfloat::detail::float_internal128_t sqrt (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SQRT (x); }
  440. inline boost::math::cstdfloat::detail::float_internal128_t trunc (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TRUNC (x); }
  441. inline boost::math::cstdfloat::detail::float_internal128_t exp (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_EXP (x); }
  442. inline boost::math::cstdfloat::detail::float_internal128_t pow (boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t a) { return ::BOOST_CSTDFLOAT_FLOAT128_POW (x, a); }
  443. inline boost::math::cstdfloat::detail::float_internal128_t pow (boost::math::cstdfloat::detail::float_internal128_t x, int a) { return ::BOOST_CSTDFLOAT_FLOAT128_POW (x, boost::math::cstdfloat::detail::float_internal128_t(a)); }
  444. inline boost::math::cstdfloat::detail::float_internal128_t log (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG (x); }
  445. inline boost::math::cstdfloat::detail::float_internal128_t log10 (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG10 (x); }
  446. inline boost::math::cstdfloat::detail::float_internal128_t sin (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SIN (x); }
  447. inline boost::math::cstdfloat::detail::float_internal128_t cos (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_COS (x); }
  448. inline boost::math::cstdfloat::detail::float_internal128_t tan (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TAN (x); }
  449. inline boost::math::cstdfloat::detail::float_internal128_t asin (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ASIN (x); }
  450. inline boost::math::cstdfloat::detail::float_internal128_t acos (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ACOS (x); }
  451. inline boost::math::cstdfloat::detail::float_internal128_t atan (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ATAN (x); }
  452. inline boost::math::cstdfloat::detail::float_internal128_t sinh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SINH (x); }
  453. inline boost::math::cstdfloat::detail::float_internal128_t cosh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_COSH (x); }
  454. inline boost::math::cstdfloat::detail::float_internal128_t tanh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TANH (x); }
  455. inline boost::math::cstdfloat::detail::float_internal128_t asinh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ASINH (x); }
  456. inline boost::math::cstdfloat::detail::float_internal128_t acosh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ACOSH (x); }
  457. inline boost::math::cstdfloat::detail::float_internal128_t atanh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ATANH (x); }
  458. inline boost::math::cstdfloat::detail::float_internal128_t fmod (boost::math::cstdfloat::detail::float_internal128_t a, boost::math::cstdfloat::detail::float_internal128_t b) { return ::BOOST_CSTDFLOAT_FLOAT128_FMOD (a, b); }
  459. inline boost::math::cstdfloat::detail::float_internal128_t atan2 (boost::math::cstdfloat::detail::float_internal128_t y, boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ATAN2 (y, x); }
  460. inline boost::math::cstdfloat::detail::float_internal128_t lgamma(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LGAMMA(x); }
  461. inline boost::math::cstdfloat::detail::float_internal128_t tgamma(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TGAMMA(x); }
  462. } } } } // boost::math::cstdfloat::detail
  463. // We will now inject the quadruple-precision <cmath> functions
  464. // into the std namespace. This is done via *using* directive.
  465. namespace std
  466. {
  467. using boost::math::cstdfloat::detail::ldexp;
  468. using boost::math::cstdfloat::detail::frexp;
  469. using boost::math::cstdfloat::detail::fabs;
  470. using boost::math::cstdfloat::detail::abs;
  471. using boost::math::cstdfloat::detail::floor;
  472. using boost::math::cstdfloat::detail::ceil;
  473. using boost::math::cstdfloat::detail::sqrt;
  474. using boost::math::cstdfloat::detail::trunc;
  475. using boost::math::cstdfloat::detail::exp;
  476. using boost::math::cstdfloat::detail::pow;
  477. using boost::math::cstdfloat::detail::log;
  478. using boost::math::cstdfloat::detail::log10;
  479. using boost::math::cstdfloat::detail::sin;
  480. using boost::math::cstdfloat::detail::cos;
  481. using boost::math::cstdfloat::detail::tan;
  482. using boost::math::cstdfloat::detail::asin;
  483. using boost::math::cstdfloat::detail::acos;
  484. using boost::math::cstdfloat::detail::atan;
  485. using boost::math::cstdfloat::detail::sinh;
  486. using boost::math::cstdfloat::detail::cosh;
  487. using boost::math::cstdfloat::detail::tanh;
  488. using boost::math::cstdfloat::detail::asinh;
  489. using boost::math::cstdfloat::detail::acosh;
  490. using boost::math::cstdfloat::detail::atanh;
  491. using boost::math::cstdfloat::detail::fmod;
  492. using boost::math::cstdfloat::detail::atan2;
  493. using boost::math::cstdfloat::detail::lgamma;
  494. using boost::math::cstdfloat::detail::tgamma;
  495. } // namespace std
  496. // We will now remove the preprocessor symbols representing quadruple-precision <cmath>
  497. // functions from the preprocessor.
  498. #undef BOOST_CSTDFLOAT_FLOAT128_LDEXP
  499. #undef BOOST_CSTDFLOAT_FLOAT128_FREXP
  500. #undef BOOST_CSTDFLOAT_FLOAT128_FABS
  501. #undef BOOST_CSTDFLOAT_FLOAT128_FLOOR
  502. #undef BOOST_CSTDFLOAT_FLOAT128_CEIL
  503. #undef BOOST_CSTDFLOAT_FLOAT128_SQRT
  504. #undef BOOST_CSTDFLOAT_FLOAT128_TRUNC
  505. #undef BOOST_CSTDFLOAT_FLOAT128_EXP
  506. #undef BOOST_CSTDFLOAT_FLOAT128_EXPM1
  507. #undef BOOST_CSTDFLOAT_FLOAT128_POW
  508. #undef BOOST_CSTDFLOAT_FLOAT128_LOG
  509. #undef BOOST_CSTDFLOAT_FLOAT128_LOG10
  510. #undef BOOST_CSTDFLOAT_FLOAT128_SIN
  511. #undef BOOST_CSTDFLOAT_FLOAT128_COS
  512. #undef BOOST_CSTDFLOAT_FLOAT128_TAN
  513. #undef BOOST_CSTDFLOAT_FLOAT128_ASIN
  514. #undef BOOST_CSTDFLOAT_FLOAT128_ACOS
  515. #undef BOOST_CSTDFLOAT_FLOAT128_ATAN
  516. #undef BOOST_CSTDFLOAT_FLOAT128_SINH
  517. #undef BOOST_CSTDFLOAT_FLOAT128_COSH
  518. #undef BOOST_CSTDFLOAT_FLOAT128_TANH
  519. #undef BOOST_CSTDFLOAT_FLOAT128_ASINH
  520. #undef BOOST_CSTDFLOAT_FLOAT128_ACOSH
  521. #undef BOOST_CSTDFLOAT_FLOAT128_ATANH
  522. #undef BOOST_CSTDFLOAT_FLOAT128_FMOD
  523. #undef BOOST_CSTDFLOAT_FLOAT128_ATAN2
  524. #undef BOOST_CSTDFLOAT_FLOAT128_LGAMMA
  525. #undef BOOST_CSTDFLOAT_FLOAT128_TGAMMA
  526. #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT (i.e., the user would like to have libquadmath support)
  527. #endif // _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_