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