From f9b8cd6bedfde49825a06bb3fe189f73cee2ecdd Mon Sep 17 00:00:00 2001 From: ckormanyos Date: Wed, 12 Feb 2014 01:00:28 +0100 Subject: [PATCH] In , patch exp() function for __float128. Add most of for __float128. Still needs some wrinkle ironing. --- include/boost/cstdfloat.hpp | 1323 +++++++++++++++++++++++++++++------ 1 file changed, 1128 insertions(+), 195 deletions(-) diff --git a/include/boost/cstdfloat.hpp b/include/boost/cstdfloat.hpp index 0e71622ee..5083436c0 100644 --- a/include/boost/cstdfloat.hpp +++ b/include/boost/cstdfloat.hpp @@ -200,12 +200,17 @@ #endif #endif - // Check if float_internal128_t from GCC's libquadmath or if (potentially) - // ICC's /Qlong-double flag is supported. - // TODO: Should we allow BOOST_MATH_USE_FLOAT128 for ICC? - // Here, we use the BOOST_MATH_USE_FLOAT128 pre-processor - // definition from . - #if (BOOST_CSTDFLOAT_HAS_FLOAT128_NATIVE_TYPE == 0) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_GCC_FLOAT128_SUPPORT) + // Check if float_internal128_t from GCC's libquadmath or if + // ICC's /Qlong-double flag is supported. Here, we use the + // BOOST_MATH_USE_FLOAT128 pre-processor definition from + // to query for libquadmath. + // If libquadmath is, quadruple-precision math is based on + // GCC's __float128 or ICC's _Quad. + + // What then follows are some rather long sections that optionally + // implements various parts of the C++ standard library for float128_t. + // Thise include , , I/O stream support, and . + #if (BOOST_CSTDFLOAT_HAS_FLOAT128_NATIVE_TYPE == 0) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT) namespace boost { namespace cstdfloat { namespace detail { #if defined(BOOST_INTEL) @@ -217,8 +222,8 @@ #define BOOST_CSTDFLOAT_FLOAT128_FLOOR __floorq #define BOOST_CSTDFLOAT_FLOAT128_CEIL __ceilq #define BOOST_CSTDFLOAT_FLOAT128_SQRT __sqrtq - #define BOOST_CSTDFLOAT_FLOAT128_TRUND __truncq - #define BOOST_CSTDFLOAT_FLOAT128_EXP __expq + #define BOOST_CSTDFLOAT_FLOAT128_TRUNC __truncq + #define BOOST_CSTDFLOAT_FLOAT128_EXP __expq_patch #define BOOST_CSTDFLOAT_FLOAT128_POW __powq #define BOOST_CSTDFLOAT_FLOAT128_LOG __logq #define BOOST_CSTDFLOAT_FLOAT128_LOG10 __log10q @@ -228,9 +233,9 @@ #define BOOST_CSTDFLOAT_FLOAT128_ASIN __asinq #define BOOST_CSTDFLOAT_FLOAT128_ACOS __acosq #define BOOST_CSTDFLOAT_FLOAT128_ATAN __atanq - #define BOOST_CSTDFLOAT_FLOAT128_SINH __sinhq - #define BOOST_CSTDFLOAT_FLOAT128_COSH __coshq - #define BOOST_CSTDFLOAT_FLOAT128_TANH __tanhq + #define BOOST_CSTDFLOAT_FLOAT128_SINH __sinhq_patch + #define BOOST_CSTDFLOAT_FLOAT128_COSH __coshq_patch + #define BOOST_CSTDFLOAT_FLOAT128_TANH __tanhq_patch #define BOOST_CSTDFLOAT_FLOAT128_FMOD __fmodq #define BOOST_CSTDFLOAT_FLOAT128_ATAN2 __atan2q #define BOOST_CSTDFLOAT_FLOAT128_LGAMMA __lgammaq @@ -245,8 +250,8 @@ #define BOOST_CSTDFLOAT_FLOAT128_FLOOR floorq #define BOOST_CSTDFLOAT_FLOAT128_CEIL ceilq #define BOOST_CSTDFLOAT_FLOAT128_SQRT sqrtq - #define BOOST_CSTDFLOAT_FLOAT128_TRUND truncq - #define BOOST_CSTDFLOAT_FLOAT128_EXP expq + #define BOOST_CSTDFLOAT_FLOAT128_TRUNC truncq + #define BOOST_CSTDFLOAT_FLOAT128_EXP expq_patch #define BOOST_CSTDFLOAT_FLOAT128_POW powq #define BOOST_CSTDFLOAT_FLOAT128_LOG logq #define BOOST_CSTDFLOAT_FLOAT128_LOG10 log10q @@ -256,9 +261,9 @@ #define BOOST_CSTDFLOAT_FLOAT128_ASIN asinq #define BOOST_CSTDFLOAT_FLOAT128_ACOS acosq #define BOOST_CSTDFLOAT_FLOAT128_ATAN atanq - #define BOOST_CSTDFLOAT_FLOAT128_SINH sinhq - #define BOOST_CSTDFLOAT_FLOAT128_COSH coshq - #define BOOST_CSTDFLOAT_FLOAT128_TANH tanhq + #define BOOST_CSTDFLOAT_FLOAT128_SINH sinhq_patch + #define BOOST_CSTDFLOAT_FLOAT128_COSH coshq_patch + #define BOOST_CSTDFLOAT_FLOAT128_TANH tanhq_patch #define BOOST_CSTDFLOAT_FLOAT128_FMOD fmodq #define BOOST_CSTDFLOAT_FLOAT128_ATAN2 atan2q #define BOOST_CSTDFLOAT_FLOAT128_LGAMMA lgammaq @@ -279,12 +284,12 @@ #define BOOST_CSTDFLOAT_FLOAT128_EPS 1.92592994438723585305597794258492732e-0034Q #define BOOST_FLOAT128_C(x) (x ## Q) - #if !defined(BOOST_CSTDFLOAT_NO_GCC_FLOAT128_LIMITS) + #if !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_LIMITS) - // For float_internal128_t, implement a specialization of std::numeric_limits<>. + // Implement a specialization of std::numeric_limits<> for quadruple-precision. // Forward declaration of quadruple-precision square root function. - extern "C" boost::cstdfloat::detail::float_internal128_t sqrtq(boost::cstdfloat::detail::float_internal128_t); + extern "C" boost::cstdfloat::detail::float_internal128_t sqrtq(boost::cstdfloat::detail::float_internal128_t) throw(); namespace std { @@ -326,9 +331,11 @@ BOOST_STATIC_CONSTEXPR float_round_style round_style = round_to_nearest; }; } - #endif + #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_LIMITS (i.e., has libquadmath support) - #if !defined(BOOST_CSTDFLOAT_NO_GCC_FLOAT128_CMATH) + #if !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_CMATH) + + #include // Implement quadruple-precision functions in the namespace // boost::cstdfloat::detail. Subsequently *use* these in the global @@ -337,34 +344,113 @@ // Begin with some forward function declarations. // Forward declarations of quadruple-precision string functions. - extern "C" int quadmath_snprintf(char *str, size_t size, const char *format, ...); - extern "C" boost::cstdfloat::detail::float_internal128_t strtoflt128(const char*, char **); + extern "C" int quadmath_snprintf(char *str, size_t size, const char *format, ...) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t strtoflt128(const char*, char **) throw(); // Forward declarations of quadruple-precision elementary functions. - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LDEXP (boost::cstdfloat::detail::float_internal128_t, int); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FREXP (boost::cstdfloat::detail::float_internal128_t, int*); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FABS (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FLOOR (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_CEIL (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SQRT (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TRUND (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_POW (boost::cstdfloat::detail::float_internal128_t, boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG10 (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SIN (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COS (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TAN (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASIN (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOS (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH (boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FMOD (boost::cstdfloat::detail::float_internal128_t, boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN2 (boost::cstdfloat::detail::float_internal128_t, boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LGAMMA(boost::cstdfloat::detail::float_internal128_t); - extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::cstdfloat::detail::float_internal128_t); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LDEXP (boost::cstdfloat::detail::float_internal128_t, int) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FREXP (boost::cstdfloat::detail::float_internal128_t, int*) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FABS (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FLOOR (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_CEIL (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SQRT (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TRUNC (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_POW (boost::cstdfloat::detail::float_internal128_t, boost::cstdfloat::detail::float_internal128_t) throw(); + inline boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::cstdfloat::detail::float_internal128_t x) + { + typedef boost::cstdfloat::detail::float_internal128_t float_type; + + // Use an order-36 polynomial approximation of the exponential function + // in the range of (-ln2 < x < ln2). Scale the argument to this range + // and multiply the result by 2^n accordingly. + + // Generate the polynomial coefficients with Mathematica(R). + + // Table[{x, Exp[x] - 1}, {x, -Log[2], Log[2], 1/180}] + // N[%, 120] + // Fit[%, {x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10, x^11, x^12, + // x^13, x^14, x^15, x^16, x^17, x^18, x^19, x^20, x^21, x^22, + // x^23, x^24, x^25, x^26, x^27, x^28, x^29, x^30, x^31, x^32, + // x^33, x^34, x^35, x^36}, x] + + BOOST_CONSTEXPR_OR_CONST float_type one_over_ln2 = float_type(BOOST_FLOAT128_C(1.44269504088896340735992468100189213742664595415299)); + const float_type x_over_ln2 = x * one_over_ln2; + + boost::int_fast32_t n; + + if (x < -1) { n = static_cast(BOOST_CSTDFLOAT_FLOAT128_CEIL (x_over_ln2)); } + else if(x > +1) { n = static_cast(BOOST_CSTDFLOAT_FLOAT128_FLOOR(x_over_ln2)); } + else { n = static_cast(0); } + + const float_type alpha = x - (n * float_type(BOOST_FLOAT128_C(0.693147180559945309417232121458176568075500134360255))); + + const float_type sum = + (((((((((((((((((((((((((((((((((((( float_type(BOOST_FLOAT128_C(2.69291698127774166063293705964720493864630783729857438187365E-42)) * alpha + + float_type(BOOST_FLOAT128_C(9.70937085471487654794114679403710456028986572118859594614033E-41))) * alpha + + float_type(BOOST_FLOAT128_C(3.38715585158055097155585505318085512156885389014410753080500E-39))) * alpha + + float_type(BOOST_FLOAT128_C(1.15162718532861050809222658798662695267019717760563645440433E-37))) * alpha + + float_type(BOOST_FLOAT128_C(3.80039074689434663295873584133017767349635602413675471702393E-36))) * alpha + + float_type(BOOST_FLOAT128_C(1.21612504934087520075905434734158045947460467096773246215239E-34))) * alpha + + float_type(BOOST_FLOAT128_C(3.76998762883139753126119821241037824830069851253295480396224E-33))) * alpha + + float_type(BOOST_FLOAT128_C(1.13099628863830344684998293828608215735777107850991029729440E-31))) * alpha + + float_type(BOOST_FLOAT128_C(3.27988923706982293204067897468714277771890104022419696770352E-30))) * alpha + + float_type(BOOST_FLOAT128_C(9.18368986379558482800593745627556950089950023355628325088207E-29))) * alpha + + float_type(BOOST_FLOAT128_C(2.47959626322479746949155352659617642905315302382639380521497E-27))) * alpha + + float_type(BOOST_FLOAT128_C(6.44695028438447337900255966737803112935639344283098705091949E-26))) * alpha + + float_type(BOOST_FLOAT128_C(1.61173757109611834904452725462599961406036904573072897122957E-24))) * alpha + + float_type(BOOST_FLOAT128_C(3.86817017063068403772269360016918092488847584660382953555804E-23))) * alpha + + float_type(BOOST_FLOAT128_C(8.89679139245057328674891109315654704307721758924206107351744E-22))) * alpha + + float_type(BOOST_FLOAT128_C(1.95729410633912612308475595397946731738088422488032228717097E-20))) * alpha + + float_type(BOOST_FLOAT128_C(4.11031762331216485847799061511674191805055663711439605760231E-19))) * alpha + + float_type(BOOST_FLOAT128_C(8.22063524662432971695598123977873600603370758794431071426640E-18))) * alpha + + float_type(BOOST_FLOAT128_C(1.56192069685862264622163643500633782667263448653185159383285E-16))) * alpha + + float_type(BOOST_FLOAT128_C(2.81145725434552076319894558300988749849555291507956994126835E-15))) * alpha + + float_type(BOOST_FLOAT128_C(4.77947733238738529743820749111754320727153728139716409114011E-14))) * alpha + + float_type(BOOST_FLOAT128_C(7.64716373181981647590113198578807092707697416852226691068627E-13))) * alpha + + float_type(BOOST_FLOAT128_C(1.14707455977297247138516979786821056670509688396295740818677E-11))) * alpha + + float_type(BOOST_FLOAT128_C(1.60590438368216145993923771701549479323291461578567184216302E-10))) * alpha + + float_type(BOOST_FLOAT128_C(2.08767569878680989792100903212014323125428376052986408239620E-09))) * alpha + + float_type(BOOST_FLOAT128_C(2.50521083854417187750521083854417187750523408006206780016659E-08))) * alpha + + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573195144226062684604E-07))) * alpha + + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573191310049321957902E-06))) * alpha + + float_type(BOOST_FLOAT128_C(0.00002480158730158730158730158730158730158730158730149317774))) * alpha + + float_type(BOOST_FLOAT128_C(0.00019841269841269841269841269841269841269841269841293575920))) * alpha + + float_type(BOOST_FLOAT128_C(0.00138888888888888888888888888888888888888888888888889071045))) * alpha + + float_type(BOOST_FLOAT128_C(0.00833333333333333333333333333333333333333333333333332986595))) * alpha + + float_type(BOOST_FLOAT128_C(0.04166666666666666666666666666666666666666666666666666664876))) * alpha + + float_type(BOOST_FLOAT128_C(0.16666666666666666666666666666666666666666666666666666669048))) * alpha + + float_type(BOOST_FLOAT128_C(0.50000000000000000000000000000000000000000000000000000000006))) * alpha + + float_type(BOOST_FLOAT128_C(0.99999999999999999999999999999999999999999999999999999999995))) * alpha + + float_type(BOOST_FLOAT128_C(1.0))); + + return sum * BOOST_CSTDFLOAT_FLOAT128_POW(float_type(2), float_type(n)); + } + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG10 (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SIN (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COS (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TAN (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASIN (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOS (boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN (boost::cstdfloat::detail::float_internal128_t) throw(); + inline boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH (boost::cstdfloat::detail::float_internal128_t x) + { + const boost::cstdfloat::detail::float_internal128_t ex = BOOST_CSTDFLOAT_FLOAT128_EXP(x); + return (ex - (1 / ex)) / 2; + } + inline boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH (boost::cstdfloat::detail::float_internal128_t x) + { + const boost::cstdfloat::detail::float_internal128_t ex = BOOST_CSTDFLOAT_FLOAT128_EXP(x); + return (ex + (1 / ex)) / 2; + } + inline boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH (boost::cstdfloat::detail::float_internal128_t x) + { + return BOOST_CSTDFLOAT_FLOAT128_SINH(x) / BOOST_CSTDFLOAT_FLOAT128_COSH(x); + } + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FMOD (boost::cstdfloat::detail::float_internal128_t, boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN2 (boost::cstdfloat::detail::float_internal128_t, boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LGAMMA(boost::cstdfloat::detail::float_internal128_t) throw(); + extern "C" boost::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::cstdfloat::detail::float_internal128_t) throw(); // Put the quadruple-precision functions in the namespace // boost::cstdfloat::detail. @@ -376,7 +462,7 @@ inline boost::cstdfloat::detail::float_internal128_t floor (boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FLOOR (x); } inline boost::cstdfloat::detail::float_internal128_t ceil (boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_CEIL (x); } inline boost::cstdfloat::detail::float_internal128_t sqrt (boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SQRT (x); } - inline boost::cstdfloat::detail::float_internal128_t trunc (boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TRUND (x); } + inline boost::cstdfloat::detail::float_internal128_t trunc (boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TRUNC (x); } inline boost::cstdfloat::detail::float_internal128_t exp (boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_EXP (x); } inline boost::cstdfloat::detail::float_internal128_t pow (boost::cstdfloat::detail::float_internal128_t x, boost::cstdfloat::detail::float_internal128_t a) { return ::BOOST_CSTDFLOAT_FLOAT128_POW (x, a); } inline boost::cstdfloat::detail::float_internal128_t log (boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG (x); } @@ -396,31 +482,6 @@ inline boost::cstdfloat::detail::float_internal128_t tgamma(boost::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TGAMMA(x); } } } } // boost::cstdfloat::detail - #undef BOOST_CSTDFLOAT_FLOAT128_LDEXP - #undef BOOST_CSTDFLOAT_FLOAT128_FREXP - #undef BOOST_CSTDFLOAT_FLOAT128_FABS - #undef BOOST_CSTDFLOAT_FLOAT128_FLOOR - #undef BOOST_CSTDFLOAT_FLOAT128_CEIL - #undef BOOST_CSTDFLOAT_FLOAT128_SQRT - #undef BOOST_CSTDFLOAT_FLOAT128_TRUND - #undef BOOST_CSTDFLOAT_FLOAT128_EXP - #undef BOOST_CSTDFLOAT_FLOAT128_POW - #undef BOOST_CSTDFLOAT_FLOAT128_LOG - #undef BOOST_CSTDFLOAT_FLOAT128_LOG10 - #undef BOOST_CSTDFLOAT_FLOAT128_SIN - #undef BOOST_CSTDFLOAT_FLOAT128_COS - #undef BOOST_CSTDFLOAT_FLOAT128_TAN - #undef BOOST_CSTDFLOAT_FLOAT128_ASIN - #undef BOOST_CSTDFLOAT_FLOAT128_ACOS - #undef BOOST_CSTDFLOAT_FLOAT128_ATAN - #undef BOOST_CSTDFLOAT_FLOAT128_SINH - #undef BOOST_CSTDFLOAT_FLOAT128_COSH - #undef BOOST_CSTDFLOAT_FLOAT128_TANH - #undef BOOST_CSTDFLOAT_FLOAT128_FMOD - #undef BOOST_CSTDFLOAT_FLOAT128_ATAN2 - #undef BOOST_CSTDFLOAT_FLOAT128_LGAMMA - #undef BOOST_CSTDFLOAT_FLOAT128_TGAMMA - using boost::cstdfloat::detail::ldexp; using boost::cstdfloat::detail::frexp; using boost::cstdfloat::detail::fabs; @@ -449,161 +510,1021 @@ // Implement quadruple-precision functions in the std namespace. namespace std { - inline boost::cstdfloat::detail::float_internal128_t ldexp (boost::cstdfloat::detail::float_internal128_t x, int n) { return boost::cstdfloat::detail::ldexp (x, n); } - inline boost::cstdfloat::detail::float_internal128_t frexp (boost::cstdfloat::detail::float_internal128_t x, int* pn) { return boost::cstdfloat::detail::frexp (x, pn); } - inline boost::cstdfloat::detail::float_internal128_t fabs (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::fabs (x); } - inline boost::cstdfloat::detail::float_internal128_t floor (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::floor (x); } - inline boost::cstdfloat::detail::float_internal128_t ceil (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::ceil (x); } - inline boost::cstdfloat::detail::float_internal128_t sqrt (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::sqrt (x); } - inline boost::cstdfloat::detail::float_internal128_t trunc (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::trunc (x); } - inline boost::cstdfloat::detail::float_internal128_t exp (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::exp (x); } - inline boost::cstdfloat::detail::float_internal128_t pow (boost::cstdfloat::detail::float_internal128_t x, boost::cstdfloat::detail::float_internal128_t a) { return boost::cstdfloat::detail::pow (x, a); } - inline boost::cstdfloat::detail::float_internal128_t log (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::log (x); } - inline boost::cstdfloat::detail::float_internal128_t log10 (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::log10 (x); } - inline boost::cstdfloat::detail::float_internal128_t sin (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::sin (x); } - inline boost::cstdfloat::detail::float_internal128_t cos (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::cos (x); } - inline boost::cstdfloat::detail::float_internal128_t tan (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::tan (x); } - inline boost::cstdfloat::detail::float_internal128_t asin (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::asin (x); } - inline boost::cstdfloat::detail::float_internal128_t acos (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::acos (x); } - inline boost::cstdfloat::detail::float_internal128_t atan (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::atan (x); } - inline boost::cstdfloat::detail::float_internal128_t sinh (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::sinh (x); } - inline boost::cstdfloat::detail::float_internal128_t cosh (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::cosh (x); } - inline boost::cstdfloat::detail::float_internal128_t tanh (boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::tanh (x); } - inline boost::cstdfloat::detail::float_internal128_t fmod (boost::cstdfloat::detail::float_internal128_t a, boost::cstdfloat::detail::float_internal128_t b) { return boost::cstdfloat::detail::fmod (a, b); } - inline boost::cstdfloat::detail::float_internal128_t atan2 (boost::cstdfloat::detail::float_internal128_t y, boost::cstdfloat::detail::float_internal128_t x) { return boost::cstdfloat::detail::atan2 (y, x); } - inline boost::cstdfloat::detail::float_internal128_t lgamma(boost::cstdfloat::detail::float_internal128_t x ) { return boost::cstdfloat::detail::lgamma(x); } - inline boost::cstdfloat::detail::float_internal128_t tgamma(boost::cstdfloat::detail::float_internal128_t x ) { return boost::cstdfloat::detail::tgamma(x); } + using boost::cstdfloat::detail::ldexp; + using boost::cstdfloat::detail::frexp; + using boost::cstdfloat::detail::fabs; + using boost::cstdfloat::detail::floor; + using boost::cstdfloat::detail::ceil; + using boost::cstdfloat::detail::sqrt; + using boost::cstdfloat::detail::trunc; + using boost::cstdfloat::detail::exp; + using boost::cstdfloat::detail::pow; + using boost::cstdfloat::detail::log; + using boost::cstdfloat::detail::log10; + using boost::cstdfloat::detail::sin; + using boost::cstdfloat::detail::cos; + using boost::cstdfloat::detail::tan; + using boost::cstdfloat::detail::asin; + using boost::cstdfloat::detail::acos; + using boost::cstdfloat::detail::atan; + using boost::cstdfloat::detail::sinh; + using boost::cstdfloat::detail::cosh; + using boost::cstdfloat::detail::tanh; + using boost::cstdfloat::detail::fmod; + using boost::cstdfloat::detail::atan2; + using boost::cstdfloat::detail::lgamma; + using boost::cstdfloat::detail::tgamma; } - #endif + #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_CMATH (i.e., has libquadmath support) - #if !defined(BOOST_CSTDFLOAT_NO_GCC_FLOAT128_IOSTREAM) + #if !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_IOSTREAM) // Implement quadruple-precision I/O stream operations. #include #include + #include #include - inline std::ostream& operator<<(std::ostream& os, const boost::cstdfloat::detail::float_internal128_t& x) + #if defined(__GNUC__) +// #if 0 + + namespace std { - #if defined(__GNUC__) - - char my_buffer[64U]; - - const int my_prec = static_cast(os.precision()); - const int my_digits = ((my_prec == 0) ? 36 : my_prec); - - const std::ios_base::fmtflags my_flags = os.flags(); - - char my_format_string[8U]; - - std::size_t my_format_string_index = 0U; - - my_format_string[my_format_string_index] = '%'; - ++my_format_string_index; - - if(my_flags & std::ios_base::showpos) { my_format_string[my_format_string_index] = '+'; ++my_format_string_index; } - if(my_flags & std::ios_base::showpoint) { my_format_string[my_format_string_index] = '#'; ++my_format_string_index; } - - my_format_string[my_format_string_index + 0U] = '.'; - my_format_string[my_format_string_index + 1U] = '*'; - my_format_string[my_format_string_index + 2U] = 'Q'; - - my_format_string_index += 3U; - - char the_notation_char; - - if (my_flags & std::ios_base::scientific) { the_notation_char = 'e'; } - else if(my_flags & std::ios_base::fixed) { the_notation_char = 'f'; } - else { the_notation_char = 'g'; } - - my_format_string[my_format_string_index + 0U] = the_notation_char; - my_format_string[my_format_string_index + 1U] = 0; - - const int v = ::quadmath_snprintf(my_buffer, - static_cast(sizeof(my_buffer)), - my_format_string, - my_digits, - x); - - if(v < 0) { BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed.")); } - - if(v >= static_cast(sizeof(my_buffer) - 1U)) + template + inline std::basic_ostream& operator<<(std::basic_ostream& os, const boost::cstdfloat::detail::float_internal128_t& x) { - // Evidently there is a really long floating-point string here, - // such as a small decimal representation. So we have to use - // dynamic memory allocation. + std::basic_ostringstream ostr; + ostr.flags(os.flags()); + ostr.imbue(os.getloc()); + ostr.precision(os.precision()); - char* my_buffer2 = static_cast(0U); + char my_buffer[64U]; - try + const int my_prec = static_cast(os.precision()); + const int my_digits = ((my_prec == 0) ? 36 : my_prec); + + const std::ios_base::fmtflags my_flags = os.flags(); + + char my_format_string[8U]; + + std::size_t my_format_string_index = 0U; + + my_format_string[my_format_string_index] = '%'; + ++my_format_string_index; + + if(my_flags & std::ios_base::showpos) { my_format_string[my_format_string_index] = '+'; ++my_format_string_index; } + if(my_flags & std::ios_base::showpoint) { my_format_string[my_format_string_index] = '#'; ++my_format_string_index; } + + my_format_string[my_format_string_index + 0U] = '.'; + my_format_string[my_format_string_index + 1U] = '*'; + my_format_string[my_format_string_index + 2U] = 'Q'; + + my_format_string_index += 3U; + + char the_notation_char; + + if (my_flags & std::ios_base::scientific) { the_notation_char = 'e'; } + else if(my_flags & std::ios_base::fixed) { the_notation_char = 'f'; } + else { the_notation_char = 'g'; } + + my_format_string[my_format_string_index + 0U] = the_notation_char; + my_format_string[my_format_string_index + 1U] = 0; + + const int v = ::quadmath_snprintf(my_buffer, + static_cast(sizeof(my_buffer)), + my_format_string, + my_digits, + x); + + if(v < 0) { BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed.")); } + + if(v >= static_cast(sizeof(my_buffer) - 1U)) { - my_buffer2 = new char[v + 3]; + // Evidently there is a really long floating-point string here, + // such as a small decimal representation. So we have to use + // dynamic memory allocation. + + char* my_buffer2 = static_cast(0U); + + try + { + my_buffer2 = new char[v + 3]; + } + catch(const std::bad_alloc&) + { + BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed while allocating memory.")); + } + + const int v2 = ::quadmath_snprintf(my_buffer2, + v + 3, + my_format_string, + my_digits, + x); + + if(v2 >= v + 3) + { + BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed.")); + } + + ostr << my_buffer2; + + delete [] my_buffer2; } - catch(const std::bad_alloc&) + else { - BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed while allocating memory.")); + ostr << my_buffer; } - const int v2 = ::quadmath_snprintf(my_buffer2, - v + 3, - my_format_string, - my_digits, - x); + return (os << ostr.str()); + } - if(v2 >= v + 3) + inline std::istream& operator>>(std::istream& is, boost::cstdfloat::detail::float_internal128_t& x) + { + std::string str; + + is >> str; + + char* p_end; + + x = strtoflt128(str.c_str(), &p_end); + + if(static_cast(p_end - str.c_str()) != static_cast(str.length())) { - BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of boost::float128_t failed.")); + BOOST_THROW_EXCEPTION(std::runtime_error("Unable to interpret input string as a boost::float128_t")); } - os << my_buffer2; + return is; + } + } - delete [] my_buffer2; + #elif defined(BOOST_INTEL) +// #elif defined(__GNUC__) - return os; + // The following string-extraction routines are based on the methodology + // used in Boost.Multiprecision by John Maddock and Christopher Kormanyos. + // This methodology has been slightly modified here for boost::float128_t. + + #include + + namespace boost { namespace cstdfloat { namespace detail { + + template + void format_float_string(string_type& str, + int my_exp, + int digits, + const std::ios_base::fmtflags f, + const bool iszero) + { + typedef typename string_type::size_type size_type; + + const bool scientific = ((f & std::ios_base::scientific) == std::ios_base::scientific); + const bool fixed = ((f & std::ios_base::fixed) == std::ios_base::fixed); + const bool showpoint = ((f & std::ios_base::showpoint) == std::ios_base::showpoint); + const bool showpos = ((f & std::ios_base::showpos) == std::ios_base::showpos); + + const bool neg = str.size() && (str[0] == '-'); + + if(neg) + { + str.erase(0, 1); + } + + if(digits == 0) + { + digits = static_cast((std::max)(str.size(), size_type(16))); + } + + if(iszero || str.empty() || (str.find_first_not_of('0') == string_type::npos)) + { + // We will be printing zero, even though the value might not + // actually be zero (it just may have been rounded to zero). + str = "0"; + + if(scientific || fixed) + { + str.append(1, '.'); + str.append(size_type(digits), '0'); + + if(scientific) + { + str.append("e+00"); + } + } + else + { + if(showpoint) + { + str.append(1, '.'); + if(digits > 1) + { + str.append(size_type(digits - 1), '0'); + } + } + } + if(neg) + { + str.insert(0U, 1U, '-'); + } + else if(showpos) + { + str.insert(0U, 1U, '+'); + } + + return; + } + + if(!fixed && !scientific && !showpoint) + { + // Suppress trailing zeros: + typename string_type::iterator pos = str.end(); + + while(pos != str.begin() && *--pos == '0') { ; } + + if(pos != str.end()) + { + ++pos; + } + + str.erase(pos, str.end()); + + if(str.empty()) + { + str = '0'; + } + } + else if(!fixed || (my_exp >= 0)) + { + // Pad out the end with zero's if we need to. + + int chars = static_cast(str.size()); + chars = digits - chars; + + if(scientific) + { + ++chars; + } + + if(chars > 0) + { + str.append(static_cast(chars), '0'); + } + } + + if(fixed || (!scientific && (my_exp >= -4) && (my_exp < digits))) + { + if((1 + my_exp) > static_cast(str.size())) + { + // Just pad out the end with zeros: + str.append(static_cast((1 + my_exp) - static_cast(str.size())), '0'); + + if(showpoint || fixed) + { + str.append("."); + } + } + else if(my_exp + 1 < static_cast(str.size())) + { + if(my_exp < 0) + { + str.insert(0U, static_cast(-1 - my_exp), '0'); + str.insert(0U, "0."); + } + else + { + // Insert the decimal point: + str.insert(static_cast(my_exp + 1), 1, '.'); + } + } + else if(showpoint || fixed) // we have exactly the digits we require to left of the point + { + str += "."; + } + + if(fixed) + { + // We may need to add trailing zeros. + int l = static_cast(str.find('.') + 1U); + l = digits - (static_cast(str.size()) - l); + + if(l > 0) + { + str.append(size_type(l), '0'); + } + } } else { - return (os << my_buffer); + // Scientific format: + if(showpoint || (str.size() > 1)) + { + str.insert(1U, 1U, '.'); + } + + str.append(1U, 'e'); + string_type e = boost::lexical_cast(std::abs(my_exp)); + + if(e.size() < 2U) + { + e.insert(0U, 2U - e.size(), '0'); + } + + if(my_exp < 0) + { + e.insert(0U, 1U, '-'); + } + else + { + e.insert(0U, 1U, '+'); + } + + str.append(e); } - #elif defined(BOOST_INTEL) - - BOOST_THROW_EXCEPTION(std::runtime_error(" does not yet support output stream operation for ICC quadruple-precision type")); - return os; - - #endif - } - - inline std::istream& operator>>(std::istream& is, boost::cstdfloat::detail::float_internal128_t& x) - { - #if defined(__GNUC__) - - std::string str; - - is >> str; - - char* p_end; - - x = strtoflt128(str.c_str(), &p_end); - - if(static_cast(p_end - str.c_str()) != static_cast(str.length())) + if(neg) { - BOOST_THROW_EXCEPTION(std::runtime_error("Unable to interpret input string as a boost::float128_t")); + str.insert(0U, 1U, '-'); + } + else if(showpos) + { + str.insert(0U, 1U, '+'); } - - return is; - - #elif defined(BOOST_INTEL) - - BOOST_THROW_EXCEPTION(std::runtime_error(" does not yet support input stream operation for ICC quadruple-precision type")); - return is; - - #endif } - #endif - #endif + template inline void eval_convert_to(type_a* pa, const float_type& cb) { *pa = static_cast(cb); } + template inline void eval_subtract (float_type& b, const type_a& a) { b -= a; } + template inline void eval_multiply (float_type& b, const type_a& a) { b *= a; } + template inline void eval_multiply (float_type& b, const float_type& cb, const float_type& cb2) { b = (cb * cb2); } + template inline void eval_divide (float_type& b, const type_a& a) { b /= a; } + template inline void eval_log10 (float_type& b, const float_type& cb) { b = std::log10(cb); } + template inline void eval_floor (float_type& b, const float_type& cb) { b = std::floor(cb); } + template inline void eval_pow (float_type& b, const float_type& cb, const float_type& cb2) { b = std::pow(cb, cb2); } + + inline void round_string_up_at(std::string& s, int pos, int& expon) + { + // This subroutine rounds up a string representation of a + // number at the given position pos. + + if(pos < 0) + { + s.insert(0U, 1U, '1'); + s.erase(s.size() - 1U); + ++expon; + } + else if(s[pos] == '9') + { + s[pos] = '0'; + round_string_up_at(s, pos - 1, expon); + } + else + { + if((pos == 0) && (s[pos] == '0') && (s.size() == 1)) + { + ++expon; + } + + ++s[pos]; + } + } + + template + std::string convert_to_string(float_type& x, + std::streamsize digits, + const std::ios_base::fmtflags f) + { + const bool iszero = (std::fabs(x) < (std::numeric_limits::min)()); + const bool isneg = (x < 0); + const bool isnan = (x != x); + const bool isinf = (std::fabs(x) > (std::numeric_limits::max)()); + + int expon = 0; + + if(digits <= 0) { digits = std::numeric_limits::max_digits10; } + + const int org_digits = static_cast(digits); + + std::string result; + + if(iszero) + { + result = "0"; + } + else if(isinf) + { + if(x < 0) + { + return "-inf"; + } + else + { + return ((f & std::ios_base::showpos) == std::ios_base::showpos) ? "+inf" : "inf"; + } + } + else if(isnan) + { + return "nan"; + } + else + { + // Start by figuring out the base-10 exponent. + + if(isneg) + { + x = -x; + } + + float_type t; + float_type ten = 10; + + eval_log10(t, x); + eval_floor(t, t); + eval_convert_to(&expon, t); + + if(-expon > std::numeric_limits::max_exponent10 - 3) + { + int e = -expon / 2; + + float_type t2; + + eval_pow(t2, ten, float_type(e)); + eval_multiply(t, t2, x); + eval_multiply(t, t2); + + if((expon & 1) != 0) + { + eval_multiply(t, ten); + } + } + else + { + eval_pow (t, ten, float_type(-expon)); + eval_multiply(t, x); + } + + // Make sure that the value lies between [1, 10), and adjust if not. + if(t < 1) + { + eval_multiply(t, 10); + + --expon; + } + else if(t >= 10) + { + eval_divide(t, 10); + + ++expon; + } + + float_type digit; + int cdigit; + + // Adjust the number of digits required based on formatting options. + if(((f & std::ios_base::fixed) == std::ios_base::fixed) && (expon != -1)) + { + digits += (expon + 1); + } + + if((f & std::ios_base::scientific) == std::ios_base::scientific) + { + ++digits; + } + + // Extract the base-10 digits one at a time. + for(int i = 0; i < digits; ++i) + { + eval_floor(digit, t); + eval_convert_to(&cdigit, digit); + + result += static_cast('0' + cdigit); + + eval_subtract(t, digit); + eval_multiply(t, ten); + } + + // Possibly round the result. + if(digits >= 0) + { + eval_floor(digit, t); + eval_convert_to(&cdigit, digit); + eval_subtract(t, digit); + + if((cdigit == 5) && (t == 0)) + { + // Use simple bankers rounding. + + if((static_cast(*result.rbegin() - '0') & 1) != 0) + { + round_string_up_at(result, static_cast(result.size() - 1U), expon); + } + } + else if(cdigit >= 5) + { + round_string_up_at(result, static_cast(result.size() - 1), expon); + } + } + } + + while((result.size() > static_cast(digits)) && result.size()) + { + // We may get here as a result of rounding. + + if(result.size() > 1U) + { + result.erase(result.size() - 1U); + } + else + { + if(expon > 0) + { + --expon; // so we put less padding in the result. + } + else + { + ++expon; + } + + ++digits; + } + } + + if(isneg) + { + result.insert(0U, 1U, '-'); + } + + format_float_string(result, expon, org_digits, f, iszero); + + return result; + } + +/* +template +void convert_from_string(Backend& b, const char* p) +{ + using default_ops::eval_multiply; + using default_ops::eval_add; + using default_ops::eval_pow; + using default_ops::eval_divide; + + typedef typename mpl::front::type ui_type; + b = ui_type(0); + if(!p || (*p == 0)) + return; + + bool is_neg = false; + bool is_neg_expon = false; + static const ui_type ten = ui_type(10); + typename Backend::exponent_type expon = 0; + int digits_seen = 0; + typedef std::numeric_limits > limits; + static const int max_digits = limits::is_specialized ? limits::max_digits10 + 1 : INT_MAX; + + if(*p == '+') ++p; + else if(*p == '-') + { + is_neg = true; + ++p; + } + if((std::strcmp(p, "nan") == 0) || (std::strcmp(p, "NaN") == 0) || (std::strcmp(p, "NAN") == 0)) + { + eval_divide(b, ui_type(0)); + if(is_neg) + b.negate(); + return; + } + if((std::strcmp(p, "inf") == 0) || (std::strcmp(p, "Inf") == 0) || (std::strcmp(p, "INF") == 0)) + { + b = ui_type(1); + eval_divide(b, ui_type(0)); + if(is_neg) + b.negate(); + return; + } + // + // Grab all the leading digits before the decimal point: + // + while(std::isdigit(*p)) + { + eval_multiply(b, ten); + eval_add(b, ui_type(*p - '0')); + ++p; + ++digits_seen; + } + if(*p == '.') + { + // + // Grab everything after the point, stop when we've seen + // enough digits, even if there are actually more available: + // + ++p; + while(std::isdigit(*p)) + { + eval_multiply(b, ten); + eval_add(b, ui_type(*p - '0')); + ++p; + --expon; + if(++digits_seen > max_digits) + break; + } + while(std::isdigit(*p)) + ++p; + } + // + // Parse the exponent: + // + if((*p == 'e') || (*p == 'E')) + { + ++p; + if(*p == '+') ++p; + else if(*p == '-') + { + is_neg_expon = true; + ++p; + } + typename Backend::exponent_type e2 = 0; + while(std::isdigit(*p)) + { + e2 *= 10; + e2 += (*p - '0'); + ++p; + } + if(is_neg_expon) + e2 = -e2; + expon += e2; + } + if(expon) + { + // Scale by 10^expon, note that 10^expon can be + // outside the range of our number type, even though the + // result is within range, if that looks likely, then split + // the calculation in two: + Backend t; + t = ten; + if(expon > limits::min_exponent10 + 2) + { + eval_pow(t, t, expon); + eval_multiply(b, t); + } + else + { + eval_pow(t, t, expon + digits_seen + 1); + eval_multiply(b, t); + t = ten; + eval_pow(t, t, -digits_seen - 1); + eval_multiply(b, t); + } + } + if(is_neg) + b.negate(); + if(*p) + { + // Unexpected input in string: + BOOST_THROW_EXCEPTION(std::runtime_error("Unexpected characters in string being interpreted as a float128.")); + } +} +*/ + + } } } // boost::cstdfloat::detail + + namespace std + { + template + inline std::basic_ostream& operator<<(std::basic_ostream& os, const boost::cstdfloat::detail::float_internal128_t& x) + { + boost::cstdfloat::detail::float_internal128_t non_const_x = x; + + const std::string str = boost::cstdfloat::detail::convert_to_string(non_const_x, + os.precision(), + os.flags()); + + return (os << str); + } + + inline std::istream& operator>>(std::istream& is, boost::cstdfloat::detail::float_internal128_t& x) + { + BOOST_THROW_EXCEPTION(std::runtime_error(" does not yet support input stream operation for ICC quadruple-precision type")); + return is; + } + } + + #endif // Use __GNUC__ or BOOST_INTEL libquadmath + + #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_IOSTREAM (i.e., has libquadmath I/O stream support) + + #if !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_COMPLEX) + + #include + + // Implement a specialization of std::complex<> for quadruple-precision. + namespace std + { + // Forward template function definitions. + #if defined(BOOST_NO_CXX11_CONSTEXPR) + template<> boost::cstdfloat::detail::float_internal128_t& real(complex&); + template<> boost::cstdfloat::detail::float_internal128_t& imag(complex&); + #else + template<> BOOST_CONSTEXPR boost::cstdfloat::detail::float_internal128_t real(const complex&); + template<> BOOST_CONSTEXPR boost::cstdfloat::detail::float_internal128_t imag(const complex&); + #endif + template<> boost::cstdfloat::detail::float_internal128_t abs (const complex&); + template<> boost::cstdfloat::detail::float_internal128_t arg (const complex&); + template<> boost::cstdfloat::detail::float_internal128_t norm(const complex&); + + template<> complex sqrt (const complex&); + template<> complex sin (const complex&); + template<> complex cos (const complex&); + template<> complex tan (const complex&); + #if !defined(BOOST_NO_CXX11_CONSTEXPR) + template<> complex asin (const complex&); + template<> complex acos (const complex&); + template<> complex atan (const complex&); + #endif + template<> complex exp (const complex&); + template<> complex log (const complex&); + template<> complex log10(const complex&); + template<> complex pow (const complex&, const boost::cstdfloat::detail::float_internal128_t&); + template<> complex pow (const complex&, const complex&); + template<> complex pow (const boost::cstdfloat::detail::float_internal128_t&, const complex&); + template<> complex sinh (const complex&); + template<> complex cosh (const complex&); + template<> complex tanh (const complex&); + #if !defined(BOOST_NO_CXX11_CONSTEXPR) + template<> complex asinh(const complex&); + template<> complex acosh(const complex&); + template<> complex atanh(const complex&); + #endif + + template<> + class complex + { + public: + typedef boost::cstdfloat::detail::float_internal128_t value_type; + + #if defined(BOOST_NO_CXX11_CONSTEXPR) + + complex(const complex& z) : re(z.re), + im(z.im) { } + + complex(value_type Re = BOOST_FLOAT128_C(0.0), + value_type Im = BOOST_FLOAT128_C(0.0)) : re(Re), + im(Im) { } + + explicit complex(const complex&); + explicit complex(const complex&); + explicit complex(const complex&); + + const value_type& real() const { return re; } + const value_type& imag() const { return im; } + value_type& real() { return re; } + value_type& imag() { return im; } + + #else + + BOOST_CONSTEXPR complex(const complex& z) : re(z.re), + im(z.im) { } + + BOOST_CONSTEXPR complex(value_type Re = BOOST_FLOAT128_C(0.0), + value_type Im = BOOST_FLOAT128_C(0.0)) : re(Re), + im(Im) { } + + explicit BOOST_CONSTEXPR complex(const complex&); + explicit BOOST_CONSTEXPR complex(const complex&); + explicit BOOST_CONSTEXPR complex(const complex&); + + BOOST_CONSTEXPR value_type real() { return re; } + BOOST_CONSTEXPR value_type imag() { return im; } + + #endif + + void real(value_type v) { re = v; } + void imag(value_type v) { im = v; } + + complex& operator+=(value_type v) + { + re += v; + return *this; + } + + complex& operator-=(value_type v) + { + re -= v; + return *this; + } + + complex& operator*=(value_type v) + { + re *= v; + im *= v; + return *this; + } + + complex& operator/=(value_type v) + { + re /= v; + im /= v; + return *this; + } + + template complex& operator+=(const complex& x) + { + re += value_type(x.re); + im += value_type(x.im); + return *this; + } + + template complex& operator-=(const complex& x) + { + re -= value_type(x.re); + im -= value_type(x.im); + return *this; + } + + template complex& operator*=(const complex& x) + { + const value_type re_x(x.re); + const value_type im_x(x.im); + const value_type tmp_re((re * re_x) - (im * im_x)); + const value_type tmp_im((re * im_x) + (im * re_x)); + re = tmp_re; + im = tmp_im; + return *this; + } + + template complex& operator/=(const complex& x) + { + const value_type re_x(x.re); + const value_type im_x(x.im); + const value_type one_over_denom = value_type(1) / std::sqrt((re_x * re_x) + (im_x * im_x)); + const value_type tmp_re = ((re * re_x) + (im * im_x)) * one_over_denom; + const value_type tmp_im = ((im * re_x) - (re * im_x)) * one_over_denom; + re = tmp_re; + im = tmp_im; + return *this; + } + + template + complex& operator=(const complex& z) + { + re = z.real(); + im = z.imag(); + return *this; + } + + complex& operator=(const value_type& v) { re = v; im = value_type(0); return *this; } + + private: + value_type re; + value_type im; + }; + + #if defined(BOOST_NO_CXX11_CONSTEXPR) + complex::complex(const complex& f) : re(boost::cstdfloat::detail::float_internal128_t( f.real())), im(boost::cstdfloat::detail::float_internal128_t( f.imag())) { } + complex::complex(const complex& d) : re(boost::cstdfloat::detail::float_internal128_t( d.real())), im(boost::cstdfloat::detail::float_internal128_t( d.imag())) { } + complex::complex(const complex& ld) : re(boost::cstdfloat::detail::float_internal128_t(ld.real())), im(boost::cstdfloat::detail::float_internal128_t(ld.imag())) { } + #else + BOOST_CONSTEXPR complex::complex(const complex& f) : re(boost::cstdfloat::detail::float_internal128_t( f.real())), im(boost::cstdfloat::detail::float_internal128_t( f.imag())) { } + BOOST_CONSTEXPR complex::complex(const complex& d) : re(boost::cstdfloat::detail::float_internal128_t( d.real())), im(boost::cstdfloat::detail::float_internal128_t( d.imag())) { } + BOOST_CONSTEXPR complex::complex(const complex& ld) : re(boost::cstdfloat::detail::float_internal128_t(ld.real())), im(boost::cstdfloat::detail::float_internal128_t(ld.imag())) { } + #endif + } // namespace std + + namespace boost { namespace cstdfloat { namespace detail { + template std::complex iz_helper__x(const std::complex& x) + { + const T tmp_r = x.real(); + return std::complex(-x.imag(), tmp_r); + } + } } } // boost::cstdfloat::detail + + namespace std + { + // 26.4.7, specific values. + #if defined(BOOST_NO_CXX11_CONSTEXPR) + template<> boost::cstdfloat::detail::float_internal128_t& real(complex& x) { return x.real(); } + template<> boost::cstdfloat::detail::float_internal128_t& imag(complex& x) { return x.imag(); } + #else + template<> BOOST_CONSTEXPR boost::cstdfloat::detail::float_internal128_t real(const complex& x) { return x.real(); } + template<> BOOST_CONSTEXPR boost::cstdfloat::detail::float_internal128_t imag(const complex& x) { return x.imag(); } + #endif + template<> boost::cstdfloat::detail::float_internal128_t abs (const complex& x) { return std::sqrt((real(x) * real(x)) + (imag(x) * imag(x))); } + template<> boost::cstdfloat::detail::float_internal128_t arg (const complex& x); + template<> boost::cstdfloat::detail::float_internal128_t norm(const complex& x) { return (real(x) * real(x)) + (imag(x) * imag(x)); } + template<> complex conj (const complex& x); + #if !defined(BOOST_NO_CXX11_CONSTEXPR) + template<> complex proj (const complex& x); + #endif + template<> complex polar(const boost::cstdfloat::detail::float_internal128_t& x, const boost::cstdfloat::detail::float_internal128_t&); + + // Global add, sub, mul, div. + template<> complex operator+(const complex& u, const complex& v) { return complex(u.real() + v.real(), u.imag() + v.imag()); } + template<> complex operator-(const complex& u, const complex& v) { return complex(u.real() - v.real(), u.imag() - v.imag()); } + template<> complex operator*(const complex& u, const complex& v) + { + const typename complex::value_type ur = u.real(); + const typename complex::value_type ui = u.imag(); + const typename complex::value_type vr = v.real(); + const typename complex::value_type vi = v.imag(); + return complex((ur * vr) - (ui * vi), (ur * vi) + (ui * vr)); + } + template<> inline complex operator/(const complex& u, const complex& v) + { + const boost::cstdfloat::detail::float_internal128_t one_over_denom = BOOST_FLOAT128_C(1.0) / std::norm(v); + const boost::cstdfloat::detail::float_internal128_t tmp_re = ((u.real() * v.real()) + (u.imag() * v.imag())) * one_over_denom; + const boost::cstdfloat::detail::float_internal128_t tmp_im = ((u.imag() * v.real()) - (u.real() * v.imag())) * one_over_denom; + return complex(tmp_re, tmp_im); + } + + template<> complex operator+(const complex& u, const boost::cstdfloat::detail::float_internal128_t& v) { return complex(u.real() + v, u.imag()); } + template<> complex operator-(const complex& u, const boost::cstdfloat::detail::float_internal128_t& v) { return complex(u.real() - v, u.imag()); } + template<> complex operator*(const complex& u, const boost::cstdfloat::detail::float_internal128_t& v) { return complex(u.real() * v, u.imag() * v); } + template<> complex operator/(const complex& u, const boost::cstdfloat::detail::float_internal128_t& v) { return complex(u.real() / v, u.imag() / v); } + + template<> complex operator+(const boost::cstdfloat::detail::float_internal128_t& u, const complex& v) { return complex(u + v.real(), v.imag()); } + template<> complex operator-(const boost::cstdfloat::detail::float_internal128_t& u, const complex& v) { return complex(u - v.real(), -v.imag()); } + template<> complex operator*(const boost::cstdfloat::detail::float_internal128_t& u, const complex& v) { return complex(u * v.real(), u * v.imag()); } + template<> complex operator/(const boost::cstdfloat::detail::float_internal128_t& u, const complex& v) { const boost::cstdfloat::detail::float_internal128_t v_norm = norm(v); return complex((u * v.real()) / v_norm, (-u * v.imag()) / v_norm); } + + // Unary plus / minus. + template<> complex operator+(const complex& u) { return u; } + template<> complex operator-(const complex& u) { return complex(-u.real(), -u.imag()); } + + // Equality and inequality. + template<> BOOST_CONSTEXPR bool operator==(const complex& u, const complex& v) { return ((u.real() == v.real()) && (u.imag() == v.imag())); } + template<> BOOST_CONSTEXPR bool operator!=(const complex& u, const complex& v) { return ((u.real() != v.real()) || (u.imag() != v.imag())); } + + template<> BOOST_CONSTEXPR bool operator==(const complex& u, const boost::cstdfloat::detail::float_internal128_t& v) { return ((u.real() == v) && (u.imag() == BOOST_FLOAT128_C(0.0))); } + template<> BOOST_CONSTEXPR bool operator!=(const complex& u, const boost::cstdfloat::detail::float_internal128_t& v) { return ((u.real() != v) || (u.imag() != BOOST_FLOAT128_C(0.0))); } + + template<> BOOST_CONSTEXPR bool operator==(const boost::cstdfloat::detail::float_internal128_t& u, const complex& v) { return ((u == v.real()) && (v.imag() == BOOST_FLOAT128_C(0.0))); } + template<> BOOST_CONSTEXPR bool operator!=(const boost::cstdfloat::detail::float_internal128_t& u, const complex& v) { return ((u != v.real()) || (v.imag() != BOOST_FLOAT128_C(0.0))); } + + // 26.4.8, transcendentals + template<> complex sqrt(const complex& x) + { + // sqrt(*this) = (s, I / 2s) for R >= 0 + // (|I| / 2s, +-s) for R < 0 + // where s = sqrt{ [ |R| + sqrt(R^2 + I^2) ] / 2 }, + // and the +- sign is the same as the sign of I. + + const boost::cstdfloat::detail::float_internal128_t zr = x.real(); + const boost::cstdfloat::detail::float_internal128_t s = std::sqrt((std::fabs(zr) + std::abs(x)) / 2); + const boost::cstdfloat::detail::float_internal128_t zi = x.imag(); + + if(zr >= 0) + { + return complex(s, (zi / s) / 2); + } + else + { + const bool imag_is_pos = (zi >= 0); + + return complex((std::fabs(zi) / s) / 2, (imag_is_pos ? s : -s)); + } + } + + template<> complex sin(const complex& x) + { + const boost::cstdfloat::detail::float_internal128_t sin_x = std::sin (x.real()); + const boost::cstdfloat::detail::float_internal128_t cos_x = std::cos (x.real()); + const boost::cstdfloat::detail::float_internal128_t sinh_y = std::sinh(x.imag()); + const boost::cstdfloat::detail::float_internal128_t cosh_y = std::cosh(x.imag()); + return complex(sin_x * cosh_y, cos_x * sinh_y); + } + + template<> complex cos(const complex& x) + { + const boost::cstdfloat::detail::float_internal128_t sin_x = std::sin (x.real()); + const boost::cstdfloat::detail::float_internal128_t cos_x = std::cos (x.real()); + const boost::cstdfloat::detail::float_internal128_t sinh_y = std::sinh(x.imag()); + const boost::cstdfloat::detail::float_internal128_t cosh_y = std::cosh(x.imag()); + return complex(cos_x * cosh_y, -(sin_x * sinh_y)); + } + + template<> complex tan(const complex& x) + { + return std::sin(x) / std::cos(x); + } + + template + inline std::basic_ostream& operator<< (std::basic_ostream& os, const std::complex&); + + template + inline std::basic_ostream& operator<< (std::basic_ostream& os, const std::complex& x) + { + std::basic_ostringstream ostr; + ostr.flags(os.flags()); + ostr.imbue(os.getloc()); + ostr.precision(os.precision()); + + ostr << '(' << x.real() << ',' << x.imag() << ')'; + + return (os << ostr.str()); + } + } // namespace std + + #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_COMPLEX (i.e., has libquadmath complex support) + + #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT (i.e., has libquadmath support) // This is the end of the preamble. Now we use the results // of the queries that have been obtained in the preamble. @@ -706,6 +1627,8 @@ BOOST_STATIC_ASSERT(std::numeric_limits::radix == 2); BOOST_STATIC_ASSERT(std::numeric_limits::digits == 11); BOOST_STATIC_ASSERT(std::numeric_limits::max_exponent == 16); + + #undef BOOST_CSTDFLOAT_HAS_FLOAT16_NATIVE_TYPE #endif #if(BOOST_CSTDFLOAT_HAS_FLOAT32_NATIVE_TYPE == 1) @@ -717,6 +1640,8 @@ BOOST_STATIC_ASSERT(std::numeric_limits::radix == 2); BOOST_STATIC_ASSERT(std::numeric_limits::digits == 24); BOOST_STATIC_ASSERT(std::numeric_limits::max_exponent == 128); + + #undef BOOST_CSTDFLOAT_HAS_FLOAT32_NATIVE_TYPE #endif #if(BOOST_CSTDFLOAT_HAS_FLOAT64_NATIVE_TYPE == 1) @@ -728,6 +1653,8 @@ BOOST_STATIC_ASSERT(std::numeric_limits::radix == 2); BOOST_STATIC_ASSERT(std::numeric_limits::digits == 53); BOOST_STATIC_ASSERT(std::numeric_limits::max_exponent == 1024); + + #undef BOOST_CSTDFLOAT_HAS_FLOAT64_NATIVE_TYPE #endif #if(BOOST_CSTDFLOAT_HAS_FLOAT80_NATIVE_TYPE == 1) @@ -739,6 +1666,8 @@ BOOST_STATIC_ASSERT(std::numeric_limits::radix == 2); BOOST_STATIC_ASSERT(std::numeric_limits::digits == 64); BOOST_STATIC_ASSERT(std::numeric_limits::max_exponent == 16384); + + #undef BOOST_CSTDFLOAT_HAS_FLOAT80_NATIVE_TYPE #endif #if(BOOST_CSTDFLOAT_HAS_FLOAT128_NATIVE_TYPE == 1) @@ -750,6 +1679,8 @@ BOOST_STATIC_ASSERT(std::numeric_limits::radix == 2); BOOST_STATIC_ASSERT(std::numeric_limits::digits == 113); BOOST_STATIC_ASSERT(std::numeric_limits::max_exponent == 16384); + + #undef BOOST_CSTDFLOAT_HAS_FLOAT128_NATIVE_TYPE #endif #if (BOOST_CSTDFLOAT_MAXIMUM_AVAILABLE_WIDTH == 16) @@ -765,6 +1696,8 @@ #else #error The maximum available floating-point width for is undefined. #endif + + #undef BOOST_CSTDFLOAT_MAXIMUM_AVAILABLE_WIDTH } // namespace boost