ceph/src/boost/boost/math/cstdfloat/cstdfloat_cmath.hpp

   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>
  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
 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   }
 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   }
 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   }
 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   }
 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   }
 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   }
 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   }
 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   }
 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 {
 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); }
 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); }
 503   inline   boost::math::cstdfloat::detail::float_internal128_t fabs  (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_FABS  (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); }
 505   inline   boost::math::cstdfloat::detail::float_internal128_t floor (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_FLOOR (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); }
 507   inline   boost::math::cstdfloat::detail::float_internal128_t sqrt  (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_SQRT  (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); }
 509   inline   boost::math::cstdfloat::detail::float_internal128_t exp   (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_EXP   (x); }
 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); }
 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)); }
 512   inline   boost::math::cstdfloat::detail::float_internal128_t log   (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_LOG   (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); }
 514   inline   boost::math::cstdfloat::detail::float_internal128_t sin   (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_SIN   (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); }
 516   inline   boost::math::cstdfloat::detail::float_internal128_t tan   (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_TAN   (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); }
 518   inline   boost::math::cstdfloat::detail::float_internal128_t acos  (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_ACOS  (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); }
 520   inline   boost::math::cstdfloat::detail::float_internal128_t sinh  (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_SINH  (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); }
 522   inline   boost::math::cstdfloat::detail::float_internal128_t tanh  (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_TANH  (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); }
 524   inline   boost::math::cstdfloat::detail::float_internal128_t acosh (boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_ACOSH (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); }
 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); }
 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); }
 528   inline   boost::math::cstdfloat::detail::float_internal128_t lgamma(boost::math::cstdfloat::detail::float_internal128_t x)                                                        { return ::BOOST_CSTDFLOAT_FLOAT128_LGAMMA(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 #if !(defined(_GLIBCXX_USE_FLOAT128) && defined(__GNUC__) && (__GNUC__ >= 7))
 540     using boost::math::cstdfloat::detail::abs;
 541 #endif
 542     using boost::math::cstdfloat::detail::floor;
 543     using boost::math::cstdfloat::detail::ceil;
 544     using boost::math::cstdfloat::detail::sqrt;
 545     using boost::math::cstdfloat::detail::trunc;
 546     using boost::math::cstdfloat::detail::exp;
 547     using boost::math::cstdfloat::detail::pow;
 548     using boost::math::cstdfloat::detail::log;
 549     using boost::math::cstdfloat::detail::log10;
 550     using boost::math::cstdfloat::detail::sin;
 551     using boost::math::cstdfloat::detail::cos;
 552     using boost::math::cstdfloat::detail::tan;
 553     using boost::math::cstdfloat::detail::asin;
 554     using boost::math::cstdfloat::detail::acos;
 555     using boost::math::cstdfloat::detail::atan;
 556     using boost::math::cstdfloat::detail::sinh;
 557     using boost::math::cstdfloat::detail::cosh;
 558     using boost::math::cstdfloat::detail::tanh;
 559     using boost::math::cstdfloat::detail::asinh;
 560     using boost::math::cstdfloat::detail::acosh;
 561     using boost::math::cstdfloat::detail::atanh;
 562     using boost::math::cstdfloat::detail::fmod;
 563     using boost::math::cstdfloat::detail::atan2;
 564     using boost::math::cstdfloat::detail::lgamma;
 565     using boost::math::cstdfloat::detail::tgamma;
 566   } // namespace std
 567
 568   // We will now remove the preprocessor symbols representing quadruple-precision <cmath>
 569   // functions from the preprocessor.
 570
 571   #undef BOOST_CSTDFLOAT_FLOAT128_LDEXP
 572   #undef BOOST_CSTDFLOAT_FLOAT128_FREXP
 573   #undef BOOST_CSTDFLOAT_FLOAT128_FABS
 574   #undef BOOST_CSTDFLOAT_FLOAT128_FLOOR
 575   #undef BOOST_CSTDFLOAT_FLOAT128_CEIL
 576   #undef BOOST_CSTDFLOAT_FLOAT128_SQRT
 577   #undef BOOST_CSTDFLOAT_FLOAT128_TRUNC
 578   #undef BOOST_CSTDFLOAT_FLOAT128_EXP
 579   #undef BOOST_CSTDFLOAT_FLOAT128_EXPM1
 580   #undef BOOST_CSTDFLOAT_FLOAT128_POW
 581   #undef BOOST_CSTDFLOAT_FLOAT128_LOG
 582   #undef BOOST_CSTDFLOAT_FLOAT128_LOG10
 583   #undef BOOST_CSTDFLOAT_FLOAT128_SIN
 584   #undef BOOST_CSTDFLOAT_FLOAT128_COS
 585   #undef BOOST_CSTDFLOAT_FLOAT128_TAN
 586   #undef BOOST_CSTDFLOAT_FLOAT128_ASIN
 587   #undef BOOST_CSTDFLOAT_FLOAT128_ACOS
 588   #undef BOOST_CSTDFLOAT_FLOAT128_ATAN
 589   #undef BOOST_CSTDFLOAT_FLOAT128_SINH
 590   #undef BOOST_CSTDFLOAT_FLOAT128_COSH
 591   #undef BOOST_CSTDFLOAT_FLOAT128_TANH
 592   #undef BOOST_CSTDFLOAT_FLOAT128_ASINH
 593   #undef BOOST_CSTDFLOAT_FLOAT128_ACOSH
 594   #undef BOOST_CSTDFLOAT_FLOAT128_ATANH
 595   #undef BOOST_CSTDFLOAT_FLOAT128_FMOD
 596   #undef BOOST_CSTDFLOAT_FLOAT128_ATAN2
 597   #undef BOOST_CSTDFLOAT_FLOAT128_LGAMMA
 598   #undef BOOST_CSTDFLOAT_FLOAT128_TGAMMA
 599
 600   #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT (i.e., the user would like to have libquadmath support)
 601
 602 #endif // _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_