From 5ce1cfc2525b06c0a9e38531813781de0281c96d Mon Sep 17 00:00:00 2001 From: DongHun Kwak Date: Thu, 5 Dec 2019 15:24:45 +0900 Subject: Imported Upstream version 1.71.0 --- boost/math/cstdfloat/cstdfloat_cmath.hpp | 1560 ++++++++++----- boost/math/differentiation/autodiff.hpp | 2053 ++++++++++++++++++++ boost/math/differentiation/autodiff_cpp11.hpp | 383 ++++ .../interpolators/cardinal_quadratic_b_spline.hpp | 57 + boost/math/interpolators/catmull_rom.hpp | 64 +- .../detail/cardinal_quadratic_b_spline_detail.hpp | 206 ++ .../detail/vector_barycentric_rational_detail.hpp | 193 ++ .../detail/whittaker_shannon_detail.hpp | 126 ++ .../interpolators/vector_barycentric_rational.hpp | 82 + boost/math/interpolators/whittaker_shannon.hpp | 47 + .../detail/ooura_fourier_integrals_detail.hpp | 651 +++++++ boost/math/quadrature/ooura_fourier_integrals.hpp | 68 + boost/math/special_functions/beta.hpp | 7 +- boost/math/special_functions/cos_pi.hpp | 7 +- .../detail/unchecked_factorial.hpp | 42 +- boost/math/special_functions/ellint_1.hpp | 14 +- boost/math/special_functions/ellint_2.hpp | 13 +- boost/math/special_functions/ellint_3.hpp | 19 +- boost/math/special_functions/ellint_d.hpp | 2 +- boost/math/special_functions/lambert_w.hpp | 28 +- boost/math/special_functions/prime.hpp | 2 +- boost/math/special_functions/sin_pi.hpp | 13 +- boost/math/special_functions/trunc.hpp | 51 + boost/math/tools/roots.hpp | 175 +- boost/math/tools/univariate_statistics.hpp | 43 +- 25 files changed, 5277 insertions(+), 629 deletions(-) create mode 100644 boost/math/differentiation/autodiff.hpp create mode 100644 boost/math/differentiation/autodiff_cpp11.hpp create mode 100644 boost/math/interpolators/cardinal_quadratic_b_spline.hpp create mode 100644 boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp create mode 100644 boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp create mode 100644 boost/math/interpolators/detail/whittaker_shannon_detail.hpp create mode 100644 boost/math/interpolators/vector_barycentric_rational.hpp create mode 100644 boost/math/interpolators/whittaker_shannon.hpp create mode 100644 boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp create mode 100644 boost/math/quadrature/ooura_fourier_integrals.hpp (limited to 'boost/math') diff --git a/boost/math/cstdfloat/cstdfloat_cmath.hpp b/boost/math/cstdfloat/cstdfloat_cmath.hpp index 83fb480dc6..ba77b2a7ad 100644 --- a/boost/math/cstdfloat/cstdfloat_cmath.hpp +++ b/boost/math/cstdfloat/cstdfloat_cmath.hpp @@ -10,595 +10,1085 @@ // Implement quadruple-precision support. #ifndef _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_ - #define _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_ - - #include - #include - - #if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT) - - #include - #include - #include - #include - #include - - #if defined(_WIN32) && defined(__GNUC__) - // Several versions of Mingw and probably cygwin too have broken - // libquadmath implementations that segfault as soon as you call - // expq or any function that depends on it. - #define BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS - #endif - - // Here is a helper function used for raising the value of a given - // floating-point type to the power of n, where n has integral type. - namespace boost { namespace math { namespace cstdfloat { namespace detail { - - template - inline float_type pown(const float_type& x, const integer_type p) - { - const bool isneg = (x < 0); - const bool isnan = (x != x); - const bool isinf = ((!isneg) ? bool(+x > (std::numeric_limits::max)()) - : bool(-x > (std::numeric_limits::max)())); - - if(isnan) { return x; } - - if(isinf) { return std::numeric_limits::quiet_NaN(); } - - const bool x_is_neg = (x < 0); - const float_type abs_x = (x_is_neg ? -x : x); - - if(p < static_cast(0)) - { - if(abs_x < (std::numeric_limits::min)()) - { - return (x_is_neg ? -std::numeric_limits::infinity() - : +std::numeric_limits::infinity()); - } - else - { - return float_type(1) / pown(x, static_cast(-p)); - } - } - - if(p == static_cast(0)) - { - return float_type(1); - } - else - { - if(p == static_cast(1)) { return x; } +#define _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_ + +#include +#include + +#if defined(BOOST_CSTDFLOAT_HAS_INTERNAL_FLOAT128_T) && defined(BOOST_MATH_USE_FLOAT128) && !defined(BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#if defined(_WIN32) && defined(__GNUC__) + // Several versions of Mingw and probably cygwin too have broken + // libquadmath implementations that segfault as soon as you call + // expq or any function that depends on it. +#define BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS +#endif - if(abs_x > (std::numeric_limits::max)()) - { - return (x_is_neg ? -std::numeric_limits::infinity() - : +std::numeric_limits::infinity()); +// Here is a helper function used for raising the value of a given +// floating-point type to the power of n, where n has integral type. +namespace boost { + namespace math { + namespace cstdfloat { + namespace detail { + + template + inline float_type pown(const float_type& x, const integer_type p) + { + const bool isneg = (x < 0); + const bool isnan = (x != x); + const bool isinf = ((!isneg) ? bool(+x > (std::numeric_limits::max)()) + : bool(-x > (std::numeric_limits::max)())); + + if (isnan) { return x; } + + if (isinf) { return std::numeric_limits::quiet_NaN(); } + + const bool x_is_neg = (x < 0); + const float_type abs_x = (x_is_neg ? -x : x); + + if (p < static_cast(0)) + { + if (abs_x < (std::numeric_limits::min)()) + { + return (x_is_neg ? -std::numeric_limits::infinity() + : +std::numeric_limits::infinity()); + } + else + { + return float_type(1) / pown(x, static_cast(-p)); + } + } + + if (p == static_cast(0)) + { + return float_type(1); + } + else + { + if (p == static_cast(1)) { return x; } + + if (abs_x > (std::numeric_limits::max)()) + { + return (x_is_neg ? -std::numeric_limits::infinity() + : +std::numeric_limits::infinity()); + } + + if (p == static_cast(2)) { return (x * x); } + else if (p == static_cast(3)) { return ((x * x) * x); } + else if (p == static_cast(4)) { const float_type x2 = (x * x); return (x2 * x2); } + else + { + // The variable xn stores the binary powers of x. + float_type result(((p % integer_type(2)) != integer_type(0)) ? x : float_type(1)); + float_type xn(x); + + integer_type p2 = p; + + while (integer_type(p2 /= 2) != integer_type(0)) + { + // Square xn for each binary power. + xn *= xn; + + const bool has_binary_power = (integer_type(p2 % integer_type(2)) != integer_type(0)); + + if (has_binary_power) + { + // Multiply the result with each binary power contained in the exponent. + result *= xn; + } + } + + return result; + } + } + } + + } } + } +} // boost::math::cstdfloat::detail + +// We will now define preprocessor symbols representing quadruple-precision functions. +#if defined(BOOST_INTEL) +#define BOOST_CSTDFLOAT_FLOAT128_LDEXP __ldexpq +#define BOOST_CSTDFLOAT_FLOAT128_FREXP __frexpq +#define BOOST_CSTDFLOAT_FLOAT128_FABS __fabsq +#define BOOST_CSTDFLOAT_FLOAT128_FLOOR __floorq +#define BOOST_CSTDFLOAT_FLOAT128_CEIL __ceilq +#if !defined(BOOST_CSTDFLOAT_FLOAT128_SQRT) +#define BOOST_CSTDFLOAT_FLOAT128_SQRT __sqrtq +#endif +#define BOOST_CSTDFLOAT_FLOAT128_TRUNC __truncq +#define BOOST_CSTDFLOAT_FLOAT128_EXP __expq +#define BOOST_CSTDFLOAT_FLOAT128_EXPM1 __expm1q +#define BOOST_CSTDFLOAT_FLOAT128_POW __powq +#define BOOST_CSTDFLOAT_FLOAT128_LOG __logq +#define BOOST_CSTDFLOAT_FLOAT128_LOG10 __log10q +#define BOOST_CSTDFLOAT_FLOAT128_SIN __sinq +#define BOOST_CSTDFLOAT_FLOAT128_COS __cosq +#define BOOST_CSTDFLOAT_FLOAT128_TAN __tanq +#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_ASINH __asinhq +#define BOOST_CSTDFLOAT_FLOAT128_ACOSH __acoshq +#define BOOST_CSTDFLOAT_FLOAT128_ATANH __atanhq +#define BOOST_CSTDFLOAT_FLOAT128_FMOD __fmodq +#define BOOST_CSTDFLOAT_FLOAT128_ATAN2 __atan2q +#define BOOST_CSTDFLOAT_FLOAT128_LGAMMA __lgammaq +#define BOOST_CSTDFLOAT_FLOAT128_TGAMMA __tgammaq +// begin more functions +#define BOOST_CSTDFLOAT_FLOAT128_REMAINDER __remainderq +#define BOOST_CSTDFLOAT_FLOAT128_REMQUO __remquoq +#define BOOST_CSTDFLOAT_FLOAT128_FMA __fmaq +#define BOOST_CSTDFLOAT_FLOAT128_FMAX __fmaxq +#define BOOST_CSTDFLOAT_FLOAT128_FMIN __fminq +#define BOOST_CSTDFLOAT_FLOAT128_FDIM __fdimq +#define BOOST_CSTDFLOAT_FLOAT128_NAN __nanq +//#define BOOST_CSTDFLOAT_FLOAT128_EXP2 __exp2q +#define BOOST_CSTDFLOAT_FLOAT128_LOG2 __log2q +#define BOOST_CSTDFLOAT_FLOAT128_LOG1P __log1pq +#define BOOST_CSTDFLOAT_FLOAT128_CBRT __cbrtq +#define BOOST_CSTDFLOAT_FLOAT128_HYPOT __hypotq +#define BOOST_CSTDFLOAT_FLOAT128_ERF __erfq +#define BOOST_CSTDFLOAT_FLOAT128_ERFC __erfcq +#define BOOST_CSTDFLOAT_FLOAT128_LLROUND __llroundq +#define BOOST_CSTDFLOAT_FLOAT128_LROUND __lroundq +#define BOOST_CSTDFLOAT_FLOAT128_ROUND __roundq +#define BOOST_CSTDFLOAT_FLOAT128_NEARBYINT __nearbyintq +#define BOOST_CSTDFLOAT_FLOAT128_LLRINT __llrintq +#define BOOST_CSTDFLOAT_FLOAT128_LRINT __lrintq +#define BOOST_CSTDFLOAT_FLOAT128_RINT __rintq +#define BOOST_CSTDFLOAT_FLOAT128_MODF __modfq +#define BOOST_CSTDFLOAT_FLOAT128_SCALBLN __scalblnq +#define BOOST_CSTDFLOAT_FLOAT128_SCALBN __scalbnq +#define BOOST_CSTDFLOAT_FLOAT128_ILOGB __ilogbq +#define BOOST_CSTDFLOAT_FLOAT128_LOGB __logbq +#define BOOST_CSTDFLOAT_FLOAT128_NEXTAFTER __nextafterq +//#define BOOST_CSTDFLOAT_FLOAT128_NEXTTOWARD __nexttowardq +#define BOOST_CSTDFLOAT_FLOAT128_COPYSIGN __copysignq +#define BOOST_CSTDFLOAT_FLOAT128_SIGNBIT __signbitq +//#define BOOST_CSTDFLOAT_FLOAT128_FPCLASSIFY __fpclassifyq +//#define BOOST_CSTDFLOAT_FLOAT128_ISFINITE __isfiniteq +#define BOOST_CSTDFLOAT_FLOAT128_ISINF __isinfq +#define BOOST_CSTDFLOAT_FLOAT128_ISNAN __isnanq +//#define BOOST_CSTDFLOAT_FLOAT128_ISNORMAL __isnormalq +//#define BOOST_CSTDFLOAT_FLOAT128_ISGREATER __isgreaterq +//#define BOOST_CSTDFLOAT_FLOAT128_ISGREATEREQUAL __isgreaterequalq +//#define BOOST_CSTDFLOAT_FLOAT128_ISLESS __islessq +//#define BOOST_CSTDFLOAT_FLOAT128_ISLESSEQUAL __islessequalq +//#define BOOST_CSTDFLOAT_FLOAT128_ISLESSGREATER __islessgreaterq +//#define BOOST_CSTDFLOAT_FLOAT128_ISUNORDERED __isunorderedq +// end more functions +#elif defined(__GNUC__) +#define BOOST_CSTDFLOAT_FLOAT128_LDEXP ldexpq +#define BOOST_CSTDFLOAT_FLOAT128_FREXP frexpq +#define BOOST_CSTDFLOAT_FLOAT128_FABS fabsq +#define BOOST_CSTDFLOAT_FLOAT128_FLOOR floorq +#define BOOST_CSTDFLOAT_FLOAT128_CEIL ceilq +#if !defined(BOOST_CSTDFLOAT_FLOAT128_SQRT) +#define BOOST_CSTDFLOAT_FLOAT128_SQRT sqrtq +#endif +#define BOOST_CSTDFLOAT_FLOAT128_TRUNC truncq +#define BOOST_CSTDFLOAT_FLOAT128_POW powq +#define BOOST_CSTDFLOAT_FLOAT128_LOG logq +#define BOOST_CSTDFLOAT_FLOAT128_LOG10 log10q +#define BOOST_CSTDFLOAT_FLOAT128_SIN sinq +#define BOOST_CSTDFLOAT_FLOAT128_COS cosq +#define BOOST_CSTDFLOAT_FLOAT128_TAN tanq +#define BOOST_CSTDFLOAT_FLOAT128_ASIN asinq +#define BOOST_CSTDFLOAT_FLOAT128_ACOS acosq +#define BOOST_CSTDFLOAT_FLOAT128_ATAN atanq +#define BOOST_CSTDFLOAT_FLOAT128_FMOD fmodq +#define BOOST_CSTDFLOAT_FLOAT128_ATAN2 atan2q +#define BOOST_CSTDFLOAT_FLOAT128_LGAMMA lgammaq +#if !defined(BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS) +#define BOOST_CSTDFLOAT_FLOAT128_EXP expq +#define BOOST_CSTDFLOAT_FLOAT128_EXPM1 expm1q +#define BOOST_CSTDFLOAT_FLOAT128_SINH sinhq +#define BOOST_CSTDFLOAT_FLOAT128_COSH coshq +#define BOOST_CSTDFLOAT_FLOAT128_TANH tanhq +#define BOOST_CSTDFLOAT_FLOAT128_ASINH asinhq +#define BOOST_CSTDFLOAT_FLOAT128_ACOSH acoshq +#define BOOST_CSTDFLOAT_FLOAT128_ATANH atanhq +#define BOOST_CSTDFLOAT_FLOAT128_TGAMMA tgammaq +#else // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS +#define BOOST_CSTDFLOAT_FLOAT128_EXP expq_patch +#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_ASINH asinhq_patch +#define BOOST_CSTDFLOAT_FLOAT128_ACOSH acoshq_patch +#define BOOST_CSTDFLOAT_FLOAT128_ATANH atanhq_patch +#define BOOST_CSTDFLOAT_FLOAT128_TGAMMA tgammaq_patch +#endif // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS +// begin more functions +#define BOOST_CSTDFLOAT_FLOAT128_REMAINDER remainderq +#define BOOST_CSTDFLOAT_FLOAT128_REMQUO remquoq +#define BOOST_CSTDFLOAT_FLOAT128_FMA fmaq +#define BOOST_CSTDFLOAT_FLOAT128_FMAX fmaxq +#define BOOST_CSTDFLOAT_FLOAT128_FMIN fminq +#define BOOST_CSTDFLOAT_FLOAT128_FDIM fdimq +#define BOOST_CSTDFLOAT_FLOAT128_NAN nanq +//#define BOOST_CSTDFLOAT_FLOAT128_EXP2 exp2q +#define BOOST_CSTDFLOAT_FLOAT128_LOG2 log2q +#define BOOST_CSTDFLOAT_FLOAT128_LOG1P log1pq +#define BOOST_CSTDFLOAT_FLOAT128_CBRT cbrtq +#define BOOST_CSTDFLOAT_FLOAT128_HYPOT hypotq +#define BOOST_CSTDFLOAT_FLOAT128_ERF erfq +#define BOOST_CSTDFLOAT_FLOAT128_ERFC erfcq +#define BOOST_CSTDFLOAT_FLOAT128_LLROUND llroundq +#define BOOST_CSTDFLOAT_FLOAT128_LROUND lroundq +#define BOOST_CSTDFLOAT_FLOAT128_ROUND roundq +#define BOOST_CSTDFLOAT_FLOAT128_NEARBYINT nearbyintq +#define BOOST_CSTDFLOAT_FLOAT128_LLRINT llrintq +#define BOOST_CSTDFLOAT_FLOAT128_LRINT lrintq +#define BOOST_CSTDFLOAT_FLOAT128_RINT rintq +#define BOOST_CSTDFLOAT_FLOAT128_MODF modfq +#define BOOST_CSTDFLOAT_FLOAT128_SCALBLN scalblnq +#define BOOST_CSTDFLOAT_FLOAT128_SCALBN scalbnq +#define BOOST_CSTDFLOAT_FLOAT128_ILOGB ilogbq +#define BOOST_CSTDFLOAT_FLOAT128_LOGB logbq +#define BOOST_CSTDFLOAT_FLOAT128_NEXTAFTER nextafterq +//#define BOOST_CSTDFLOAT_FLOAT128_NEXTTOWARD nexttowardq +#define BOOST_CSTDFLOAT_FLOAT128_COPYSIGN copysignq +#define BOOST_CSTDFLOAT_FLOAT128_SIGNBIT signbitq +//#define BOOST_CSTDFLOAT_FLOAT128_FPCLASSIFY fpclassifyq +//#define BOOST_CSTDFLOAT_FLOAT128_ISFINITE isfiniteq +#define BOOST_CSTDFLOAT_FLOAT128_ISINF isinfq +#define BOOST_CSTDFLOAT_FLOAT128_ISNAN isnanq +//#define BOOST_CSTDFLOAT_FLOAT128_ISNORMAL isnormalq +//#define BOOST_CSTDFLOAT_FLOAT128_ISGREATER isgreaterq +//#define BOOST_CSTDFLOAT_FLOAT128_ISGREATEREQUAL isgreaterequalq +//#define BOOST_CSTDFLOAT_FLOAT128_ISLESS islessq +//#define BOOST_CSTDFLOAT_FLOAT128_ISLESSEQUAL islessequalq +//#define BOOST_CSTDFLOAT_FLOAT128_ISLESSGREATER islessgreaterq +//#define BOOST_CSTDFLOAT_FLOAT128_ISUNORDERED isunorderedq +// end more functions +#endif - if (p == static_cast(2)) { return (x * x); } - else if(p == static_cast(3)) { return ((x * x) * x); } - else if(p == static_cast(4)) { const float_type x2 = (x * x); return (x2 * x2); } - else - { - // The variable xn stores the binary powers of x. - float_type result(((p % integer_type(2)) != integer_type(0)) ? x : float_type(1)); - float_type xn (x); - - integer_type p2 = p; - - while(integer_type(p2 /= 2) != integer_type(0)) - { - // Square xn for each binary power. - xn *= xn; - - const bool has_binary_power = (integer_type(p2 % integer_type(2)) != integer_type(0)); - - if(has_binary_power) - { - // Multiply the result with each binary power contained in the exponent. - result *= xn; - } - } - - return result; - } - } - } - - } } } } // boost::math::cstdfloat::detail - - // We will now define preprocessor symbols representing quadruple-precision functions. - #if defined(BOOST_INTEL) - #define BOOST_CSTDFLOAT_FLOAT128_LDEXP __ldexpq - #define BOOST_CSTDFLOAT_FLOAT128_FREXP __frexpq - #define BOOST_CSTDFLOAT_FLOAT128_FABS __fabsq - #define BOOST_CSTDFLOAT_FLOAT128_FLOOR __floorq - #define BOOST_CSTDFLOAT_FLOAT128_CEIL __ceilq - #if !defined(BOOST_CSTDFLOAT_FLOAT128_SQRT) - #define BOOST_CSTDFLOAT_FLOAT128_SQRT __sqrtq - #endif - #define BOOST_CSTDFLOAT_FLOAT128_TRUNC __truncq - #define BOOST_CSTDFLOAT_FLOAT128_EXP __expq - #define BOOST_CSTDFLOAT_FLOAT128_EXPM1 __expm1q - #define BOOST_CSTDFLOAT_FLOAT128_POW __powq - #define BOOST_CSTDFLOAT_FLOAT128_LOG __logq - #define BOOST_CSTDFLOAT_FLOAT128_LOG10 __log10q - #define BOOST_CSTDFLOAT_FLOAT128_SIN __sinq - #define BOOST_CSTDFLOAT_FLOAT128_COS __cosq - #define BOOST_CSTDFLOAT_FLOAT128_TAN __tanq - #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_ASINH __asinhq - #define BOOST_CSTDFLOAT_FLOAT128_ACOSH __acoshq - #define BOOST_CSTDFLOAT_FLOAT128_ATANH __atanhq - #define BOOST_CSTDFLOAT_FLOAT128_FMOD __fmodq - #define BOOST_CSTDFLOAT_FLOAT128_ATAN2 __atan2q - #define BOOST_CSTDFLOAT_FLOAT128_LGAMMA __lgammaq - #define BOOST_CSTDFLOAT_FLOAT128_TGAMMA __tgammaq - #elif defined(__GNUC__) - #define BOOST_CSTDFLOAT_FLOAT128_LDEXP ldexpq - #define BOOST_CSTDFLOAT_FLOAT128_FREXP frexpq - #define BOOST_CSTDFLOAT_FLOAT128_FABS fabsq - #define BOOST_CSTDFLOAT_FLOAT128_FLOOR floorq - #define BOOST_CSTDFLOAT_FLOAT128_CEIL ceilq - #if !defined(BOOST_CSTDFLOAT_FLOAT128_SQRT) - #define BOOST_CSTDFLOAT_FLOAT128_SQRT sqrtq - #endif - #define BOOST_CSTDFLOAT_FLOAT128_TRUNC truncq - #define BOOST_CSTDFLOAT_FLOAT128_POW powq - #define BOOST_CSTDFLOAT_FLOAT128_LOG logq - #define BOOST_CSTDFLOAT_FLOAT128_LOG10 log10q - #define BOOST_CSTDFLOAT_FLOAT128_SIN sinq - #define BOOST_CSTDFLOAT_FLOAT128_COS cosq - #define BOOST_CSTDFLOAT_FLOAT128_TAN tanq - #define BOOST_CSTDFLOAT_FLOAT128_ASIN asinq - #define BOOST_CSTDFLOAT_FLOAT128_ACOS acosq - #define BOOST_CSTDFLOAT_FLOAT128_ATAN atanq - #define BOOST_CSTDFLOAT_FLOAT128_FMOD fmodq - #define BOOST_CSTDFLOAT_FLOAT128_ATAN2 atan2q - #define BOOST_CSTDFLOAT_FLOAT128_LGAMMA lgammaq - #if !defined(BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS) - #define BOOST_CSTDFLOAT_FLOAT128_EXP expq - #define BOOST_CSTDFLOAT_FLOAT128_EXPM1 expm1q_internal - #define BOOST_CSTDFLOAT_FLOAT128_SINH sinhq - #define BOOST_CSTDFLOAT_FLOAT128_COSH coshq - #define BOOST_CSTDFLOAT_FLOAT128_TANH tanhq - #define BOOST_CSTDFLOAT_FLOAT128_ASINH asinhq - #define BOOST_CSTDFLOAT_FLOAT128_ACOSH acoshq - #define BOOST_CSTDFLOAT_FLOAT128_ATANH atanhq - #define BOOST_CSTDFLOAT_FLOAT128_TGAMMA tgammaq - #else // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS - #define BOOST_CSTDFLOAT_FLOAT128_EXP expq_patch - #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_ASINH asinhq_patch - #define BOOST_CSTDFLOAT_FLOAT128_ACOSH acoshq_patch - #define BOOST_CSTDFLOAT_FLOAT128_ATANH atanhq_patch - #define BOOST_CSTDFLOAT_FLOAT128_TGAMMA tgammaq_patch - #endif // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS - #endif - - // Implement quadruple-precision functions in the namespace - // boost::math::cstdfloat::detail. Subsequently inject these into the - // std namespace via *using* directive. - - // Begin with some forward function declarations. Also implement patches - // for compilers that have broken float128 exponential functions. - - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LDEXP (boost::math::cstdfloat::detail::float_internal128_t, int) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FREXP (boost::math::cstdfloat::detail::float_internal128_t, int*) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FABS (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FLOOR (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_CEIL (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SQRT (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TRUNC (boost::math::cstdfloat::detail::float_internal128_t) throw(); - 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(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG10 (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SIN (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COS (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TAN (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASIN (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOS (boost::math::cstdfloat::detail::float_internal128_t) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN (boost::math::cstdfloat::detail::float_internal128_t) throw(); - 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(); - 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(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LGAMMA(boost::math::cstdfloat::detail::float_internal128_t) throw(); - - #if !defined(BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS) - - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::math::cstdfloat::detail::float_internal128_t x) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH (boost::math::cstdfloat::detail::float_internal128_t x) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH (boost::math::cstdfloat::detail::float_internal128_t x) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH (boost::math::cstdfloat::detail::float_internal128_t x) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASINH (boost::math::cstdfloat::detail::float_internal128_t x) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOSH (boost::math::cstdfloat::detail::float_internal128_t x) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATANH (boost::math::cstdfloat::detail::float_internal128_t x) throw(); - extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::math::cstdfloat::detail::float_internal128_t x) throw(); - - #else // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS - - // Forward declaration of the patched exponent function, exp(x). - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::math::cstdfloat::detail::float_internal128_t x); - - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXPM1 (boost::math::cstdfloat::detail::float_internal128_t x) - { - // Compute exp(x) - 1 for x small. - - // Use an order-36 polynomial approximation of the exponential function - // in the range of (-ln2 < x < ln2). Scale the argument to this range - // and subsequently multiply the result by 2^n accordingly. - - // Derive the polynomial coefficients with Mathematica(R) by generating - // a table of high-precision values of exp(x) in the range (-ln2 < x < ln2) - // and subsequently applying the built-in *Fit* function. - - // 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] - - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - - float_type sum; - - if(x > BOOST_FLOAT128_C(0.693147180559945309417232121458176568075500134360255)) - { +// Implement quadruple-precision functions in the namespace +// boost::math::cstdfloat::detail. Subsequently inject these into the +// std namespace via *using* directive. + +// Begin with some forward function declarations. Also implement patches +// for compilers that have broken float128 exponential functions. + +extern "C" int quadmath_snprintf(char*, std::size_t, const char*, ...) throw(); + +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LDEXP(boost::math::cstdfloat::detail::float_internal128_t, int) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FREXP(boost::math::cstdfloat::detail::float_internal128_t, int*) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FABS(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FLOOR(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_CEIL(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SQRT(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TRUNC(boost::math::cstdfloat::detail::float_internal128_t) throw(); +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(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG10(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SIN(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COS(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TAN(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASIN(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOS(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATAN(boost::math::cstdfloat::detail::float_internal128_t) throw(); +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(); +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(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LGAMMA(boost::math::cstdfloat::detail::float_internal128_t) throw(); + +// begin more functions +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_REMAINDER(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_REMQUO(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t, int*) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FMA(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FMAX(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FMIN(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_FDIM(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_NAN(const char*) throw(); +//extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP2 (boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG2(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOG1P(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_CBRT(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_HYPOT(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ERF(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ERFC(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" long long int BOOST_CSTDFLOAT_FLOAT128_LLROUND(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" long int BOOST_CSTDFLOAT_FLOAT128_LROUND(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ROUND(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_NEARBYINT(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" long long int BOOST_CSTDFLOAT_FLOAT128_LLRINT(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" long int BOOST_CSTDFLOAT_FLOAT128_LRINT(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_RINT(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_MODF(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t*) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SCALBLN(boost::math::cstdfloat::detail::float_internal128_t, long int) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SCALBN(boost::math::cstdfloat::detail::float_internal128_t, int) throw(); +extern "C" int BOOST_CSTDFLOAT_FLOAT128_ILOGB(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_LOGB(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_NEXTAFTER(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_NEXTTOWARD (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COPYSIGN(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" int BOOST_CSTDFLOAT_FLOAT128_SIGNBIT(boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_FPCLASSIFY (boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISFINITE (boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISINF(boost::math::cstdfloat::detail::float_internal128_t) throw(); +extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISNAN(boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ISNORMAL (boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISGREATER (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISGREATEREQUAL(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISLESS (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISLESSEQUAL (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISLESSGREATER(boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); +//extern "C" int BOOST_CSTDFLOAT_FLOAT128_ISUNORDERED (boost::math::cstdfloat::detail::float_internal128_t, boost::math::cstdfloat::detail::float_internal128_t) throw(); + // end more functions + +#if !defined(BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS) + +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXPM1(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASINH(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOSH(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATANH(boost::math::cstdfloat::detail::float_internal128_t x) throw(); +extern "C" boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::math::cstdfloat::detail::float_internal128_t x) throw(); + +#else // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS + +// Forward declaration of the patched exponent function, exp(x). +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP(boost::math::cstdfloat::detail::float_internal128_t x); + +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXPM1(boost::math::cstdfloat::detail::float_internal128_t x) +{ + // Compute exp(x) - 1 for x small. + + // Use an order-36 polynomial approximation of the exponential function + // in the range of (-ln2 < x < ln2). Scale the argument to this range + // and subsequently multiply the result by 2^n accordingly. + + // Derive the polynomial coefficients with Mathematica(R) by generating + // a table of high-precision values of exp(x) in the range (-ln2 < x < ln2) + // and subsequently applying the built-in *Fit* function. + + // 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] + + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + + float_type sum; + + if (x > BOOST_FLOAT128_C(0.693147180559945309417232121458176568075500134360255)) + { sum = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x) - float_type(1); - } - else - { + } + else + { // Compute the polynomial approximation of exp(alpha). - sum = (((((((((((((((((((((((((((((((((((( float_type(BOOST_FLOAT128_C(2.69291698127774166063293705964720493864630783729857438187365E-42)) * x - + float_type(BOOST_FLOAT128_C(9.70937085471487654794114679403710456028986572118859594614033E-41))) * x - + float_type(BOOST_FLOAT128_C(3.38715585158055097155585505318085512156885389014410753080500E-39))) * x - + float_type(BOOST_FLOAT128_C(1.15162718532861050809222658798662695267019717760563645440433E-37))) * x - + float_type(BOOST_FLOAT128_C(3.80039074689434663295873584133017767349635602413675471702393E-36))) * x - + float_type(BOOST_FLOAT128_C(1.21612504934087520075905434734158045947460467096773246215239E-34))) * x - + float_type(BOOST_FLOAT128_C(3.76998762883139753126119821241037824830069851253295480396224E-33))) * x - + float_type(BOOST_FLOAT128_C(1.13099628863830344684998293828608215735777107850991029729440E-31))) * x - + float_type(BOOST_FLOAT128_C(3.27988923706982293204067897468714277771890104022419696770352E-30))) * x - + float_type(BOOST_FLOAT128_C(9.18368986379558482800593745627556950089950023355628325088207E-29))) * x - + float_type(BOOST_FLOAT128_C(2.47959626322479746949155352659617642905315302382639380521497E-27))) * x - + float_type(BOOST_FLOAT128_C(6.44695028438447337900255966737803112935639344283098705091949E-26))) * x - + float_type(BOOST_FLOAT128_C(1.61173757109611834904452725462599961406036904573072897122957E-24))) * x - + float_type(BOOST_FLOAT128_C(3.86817017063068403772269360016918092488847584660382953555804E-23))) * x - + float_type(BOOST_FLOAT128_C(8.89679139245057328674891109315654704307721758924206107351744E-22))) * x - + float_type(BOOST_FLOAT128_C(1.95729410633912612308475595397946731738088422488032228717097E-20))) * x - + float_type(BOOST_FLOAT128_C(4.11031762331216485847799061511674191805055663711439605760231E-19))) * x - + float_type(BOOST_FLOAT128_C(8.22063524662432971695598123977873600603370758794431071426640E-18))) * x - + float_type(BOOST_FLOAT128_C(1.56192069685862264622163643500633782667263448653185159383285E-16))) * x - + float_type(BOOST_FLOAT128_C(2.81145725434552076319894558300988749849555291507956994126835E-15))) * x - + float_type(BOOST_FLOAT128_C(4.77947733238738529743820749111754320727153728139716409114011E-14))) * x - + float_type(BOOST_FLOAT128_C(7.64716373181981647590113198578807092707697416852226691068627E-13))) * x - + float_type(BOOST_FLOAT128_C(1.14707455977297247138516979786821056670509688396295740818677E-11))) * x - + float_type(BOOST_FLOAT128_C(1.60590438368216145993923771701549479323291461578567184216302E-10))) * x - + float_type(BOOST_FLOAT128_C(2.08767569878680989792100903212014323125428376052986408239620E-09))) * x - + float_type(BOOST_FLOAT128_C(2.50521083854417187750521083854417187750523408006206780016659E-08))) * x - + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573195144226062684604E-07))) * x - + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573191310049321957902E-06))) * x - + float_type(BOOST_FLOAT128_C(0.00002480158730158730158730158730158730158730158730149317774))) * x - + float_type(BOOST_FLOAT128_C(0.00019841269841269841269841269841269841269841269841293575920))) * x - + float_type(BOOST_FLOAT128_C(0.00138888888888888888888888888888888888888888888888889071045))) * x - + float_type(BOOST_FLOAT128_C(0.00833333333333333333333333333333333333333333333333332986595))) * x - + float_type(BOOST_FLOAT128_C(0.04166666666666666666666666666666666666666666666666666664876))) * x - + float_type(BOOST_FLOAT128_C(0.16666666666666666666666666666666666666666666666666666669048))) * x - + float_type(BOOST_FLOAT128_C(0.50000000000000000000000000000000000000000000000000000000006))) * x - + float_type(BOOST_FLOAT128_C(0.99999999999999999999999999999999999999999999999999999999995))) * x); - } - - return sum; - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP (boost::math::cstdfloat::detail::float_internal128_t x) - { - // Patch the expq() function for a subset of broken GCC compilers - // like GCC 4.7, 4.8 on MinGW. - - // Use an order-36 polynomial approximation of the exponential function - // in the range of (-ln2 < x < ln2). Scale the argument to this range - // and subsequently multiply the result by 2^n accordingly. - - // Derive the polynomial coefficients with Mathematica(R) by generating - // a table of high-precision values of exp(x) in the range (-ln2 < x < ln2) - // and subsequently applying the built-in *Fit* function. - - // 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] - - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - - // Scale the argument x to the range (-ln2 < x < ln2). - 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 != x) - { + sum = ((((((((((((((((((((((((((((((((((((float_type(BOOST_FLOAT128_C(2.69291698127774166063293705964720493864630783729857438187365E-42)) * x + + float_type(BOOST_FLOAT128_C(9.70937085471487654794114679403710456028986572118859594614033E-41))) * x + + float_type(BOOST_FLOAT128_C(3.38715585158055097155585505318085512156885389014410753080500E-39))) * x + + float_type(BOOST_FLOAT128_C(1.15162718532861050809222658798662695267019717760563645440433E-37))) * x + + float_type(BOOST_FLOAT128_C(3.80039074689434663295873584133017767349635602413675471702393E-36))) * x + + float_type(BOOST_FLOAT128_C(1.21612504934087520075905434734158045947460467096773246215239E-34))) * x + + float_type(BOOST_FLOAT128_C(3.76998762883139753126119821241037824830069851253295480396224E-33))) * x + + float_type(BOOST_FLOAT128_C(1.13099628863830344684998293828608215735777107850991029729440E-31))) * x + + float_type(BOOST_FLOAT128_C(3.27988923706982293204067897468714277771890104022419696770352E-30))) * x + + float_type(BOOST_FLOAT128_C(9.18368986379558482800593745627556950089950023355628325088207E-29))) * x + + float_type(BOOST_FLOAT128_C(2.47959626322479746949155352659617642905315302382639380521497E-27))) * x + + float_type(BOOST_FLOAT128_C(6.44695028438447337900255966737803112935639344283098705091949E-26))) * x + + float_type(BOOST_FLOAT128_C(1.61173757109611834904452725462599961406036904573072897122957E-24))) * x + + float_type(BOOST_FLOAT128_C(3.86817017063068403772269360016918092488847584660382953555804E-23))) * x + + float_type(BOOST_FLOAT128_C(8.89679139245057328674891109315654704307721758924206107351744E-22))) * x + + float_type(BOOST_FLOAT128_C(1.95729410633912612308475595397946731738088422488032228717097E-20))) * x + + float_type(BOOST_FLOAT128_C(4.11031762331216485847799061511674191805055663711439605760231E-19))) * x + + float_type(BOOST_FLOAT128_C(8.22063524662432971695598123977873600603370758794431071426640E-18))) * x + + float_type(BOOST_FLOAT128_C(1.56192069685862264622163643500633782667263448653185159383285E-16))) * x + + float_type(BOOST_FLOAT128_C(2.81145725434552076319894558300988749849555291507956994126835E-15))) * x + + float_type(BOOST_FLOAT128_C(4.77947733238738529743820749111754320727153728139716409114011E-14))) * x + + float_type(BOOST_FLOAT128_C(7.64716373181981647590113198578807092707697416852226691068627E-13))) * x + + float_type(BOOST_FLOAT128_C(1.14707455977297247138516979786821056670509688396295740818677E-11))) * x + + float_type(BOOST_FLOAT128_C(1.60590438368216145993923771701549479323291461578567184216302E-10))) * x + + float_type(BOOST_FLOAT128_C(2.08767569878680989792100903212014323125428376052986408239620E-09))) * x + + float_type(BOOST_FLOAT128_C(2.50521083854417187750521083854417187750523408006206780016659E-08))) * x + + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573195144226062684604E-07))) * x + + float_type(BOOST_FLOAT128_C(2.75573192239858906525573192239858906525573191310049321957902E-06))) * x + + float_type(BOOST_FLOAT128_C(0.00002480158730158730158730158730158730158730158730149317774))) * x + + float_type(BOOST_FLOAT128_C(0.00019841269841269841269841269841269841269841269841293575920))) * x + + float_type(BOOST_FLOAT128_C(0.00138888888888888888888888888888888888888888888888889071045))) * x + + float_type(BOOST_FLOAT128_C(0.00833333333333333333333333333333333333333333333333332986595))) * x + + float_type(BOOST_FLOAT128_C(0.04166666666666666666666666666666666666666666666666666664876))) * x + + float_type(BOOST_FLOAT128_C(0.16666666666666666666666666666666666666666666666666666669048))) * x + + float_type(BOOST_FLOAT128_C(0.50000000000000000000000000000000000000000000000000000000006))) * x + + float_type(BOOST_FLOAT128_C(0.99999999999999999999999999999999999999999999999999999999995))) * x); + } + + return sum; +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_EXP(boost::math::cstdfloat::detail::float_internal128_t x) +{ + // Patch the expq() function for a subset of broken GCC compilers + // like GCC 4.7, 4.8 on MinGW. + + // Use an order-36 polynomial approximation of the exponential function + // in the range of (-ln2 < x < ln2). Scale the argument to this range + // and subsequently multiply the result by 2^n accordingly. + + // Derive the polynomial coefficients with Mathematica(R) by generating + // a table of high-precision values of exp(x) in the range (-ln2 < x < ln2) + // and subsequently applying the built-in *Fit* function. + + // 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] + + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + + // Scale the argument x to the range (-ln2 < x < ln2). + 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 != x) + { // The argument is NaN. return std::numeric_limits::quiet_NaN(); - } - else if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x) > BOOST_FLOAT128_C(+0.693147180559945309417232121458176568075500134360255)) - { + } + else if (::BOOST_CSTDFLOAT_FLOAT128_FABS(x) > BOOST_FLOAT128_C(+0.693147180559945309417232121458176568075500134360255)) + { // The absolute value of the argument exceeds ln2. n = static_cast(::BOOST_CSTDFLOAT_FLOAT128_FLOOR(x_over_ln2)); - } - else if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x) < BOOST_FLOAT128_C(+0.693147180559945309417232121458176568075500134360255)) - { + } + else if (::BOOST_CSTDFLOAT_FLOAT128_FABS(x) < BOOST_FLOAT128_C(+0.693147180559945309417232121458176568075500134360255)) + { // The absolute value of the argument is less than ln2. n = static_cast(0); - } - else - { + } + else + { // The absolute value of the argument is exactly equal to ln2 (in the sense of floating-point equality). return float_type(2); - } + } - // Check if the argument is very near an integer. - const float_type floor_of_x = ::BOOST_CSTDFLOAT_FLOAT128_FLOOR(x); + // Check if the argument is very near an integer. + const float_type floor_of_x = ::BOOST_CSTDFLOAT_FLOAT128_FLOOR(x); - if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x - floor_of_x) < float_type(BOOST_CSTDFLOAT_FLOAT128_EPS)) - { + if (::BOOST_CSTDFLOAT_FLOAT128_FABS(x - floor_of_x) < float_type(BOOST_CSTDFLOAT_FLOAT128_EPS)) + { // Return e^n for arguments very near an integer. return boost::math::cstdfloat::detail::pown(BOOST_FLOAT128_C(2.71828182845904523536028747135266249775724709369996), static_cast(floor_of_x)); - } - - // Compute the scaled argument alpha. - const float_type alpha = x - (n * BOOST_FLOAT128_C(0.693147180559945309417232121458176568075500134360255)); - - // Compute the polynomial approximation of expm1(alpha) and add to it - // in order to obtain the scaled result. - const float_type scaled_result = ::BOOST_CSTDFLOAT_FLOAT128_EXPM1(alpha) + float_type(1); - - // Rescale the result and return it. - return scaled_result * boost::math::cstdfloat::detail::pown(float_type(2), n); - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH (boost::math::cstdfloat::detail::float_internal128_t x) - { - // Patch the sinhq() function for a subset of broken GCC compilers - // like GCC 4.7, 4.8 on MinGW. - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - - // Here, we use the following: - // Set: ex = exp(x) - // Set: em1 = expm1(x) - // Then - // sinh(x) = (ex - 1/ex) / 2 ; for |x| >= 1 - // sinh(x) = (2em1 + em1^2) / (2ex) ; for |x| < 1 - - const float_type ex = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x); - - if(::BOOST_CSTDFLOAT_FLOAT128_FABS(x) < float_type(+1)) - { + } + + // Compute the scaled argument alpha. + const float_type alpha = x - (n * BOOST_FLOAT128_C(0.693147180559945309417232121458176568075500134360255)); + + // Compute the polynomial approximation of expm1(alpha) and add to it + // in order to obtain the scaled result. + const float_type scaled_result = ::BOOST_CSTDFLOAT_FLOAT128_EXPM1(alpha) + float_type(1); + + // Rescale the result and return it. + return scaled_result * boost::math::cstdfloat::detail::pown(float_type(2), n); +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_SINH(boost::math::cstdfloat::detail::float_internal128_t x) +{ + // Patch the sinhq() function for a subset of broken GCC compilers + // like GCC 4.7, 4.8 on MinGW. + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + + // Here, we use the following: + // Set: ex = exp(x) + // Set: em1 = expm1(x) + // Then + // sinh(x) = (ex - 1/ex) / 2 ; for |x| >= 1 + // sinh(x) = (2em1 + em1^2) / (2ex) ; for |x| < 1 + + const float_type ex = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x); + + if (::BOOST_CSTDFLOAT_FLOAT128_FABS(x) < float_type(+1)) + { const float_type em1 = ::BOOST_CSTDFLOAT_FLOAT128_EXPM1(x); return ((em1 * 2) + (em1 * em1)) / (ex * 2); - } - else - { + } + else + { return (ex - (float_type(1) / ex)) / 2; - } - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH (boost::math::cstdfloat::detail::float_internal128_t x) - { - // Patch the coshq() function for a subset of broken GCC compilers - // like GCC 4.7, 4.8 on MinGW. - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - const float_type ex = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x); - return (ex + (float_type(1) / ex)) / 2; - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH (boost::math::cstdfloat::detail::float_internal128_t x) - { - // Patch the tanhq() function for a subset of broken GCC compilers - // like GCC 4.7, 4.8 on MinGW. - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - const float_type ex_plus = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x); - const float_type ex_minus = (float_type(1) / ex_plus); - return (ex_plus - ex_minus) / (ex_plus + ex_minus); - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASINH(boost::math::cstdfloat::detail::float_internal128_t x) throw() - { - // Patch the asinh() function since quadmath does not have it. - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - return ::BOOST_CSTDFLOAT_FLOAT128_LOG(x + ::BOOST_CSTDFLOAT_FLOAT128_SQRT((x * x) + float_type(1))); - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOSH(boost::math::cstdfloat::detail::float_internal128_t x) throw() - { - // Patch the acosh() function since quadmath does not have it. - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - const float_type zp(x + float_type(1)); - const float_type zm(x - float_type(1)); - - return ::BOOST_CSTDFLOAT_FLOAT128_LOG(x + (zp * ::BOOST_CSTDFLOAT_FLOAT128_SQRT(zm / zp))); - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATANH(boost::math::cstdfloat::detail::float_internal128_t x) throw() - { - // Patch the atanh() function since quadmath does not have it. - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - return ( ::BOOST_CSTDFLOAT_FLOAT128_LOG(float_type(1) + x) - - ::BOOST_CSTDFLOAT_FLOAT128_LOG(float_type(1) - x)) / 2; - } - inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::math::cstdfloat::detail::float_internal128_t x) throw() - { - // Patch the tgammaq() function for a subset of broken GCC compilers - // like GCC 4.7, 4.8 on MinGW. - typedef boost::math::cstdfloat::detail::float_internal128_t float_type; - - if(x > float_type(0)) - { + } +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_COSH(boost::math::cstdfloat::detail::float_internal128_t x) +{ + // Patch the coshq() function for a subset of broken GCC compilers + // like GCC 4.7, 4.8 on MinGW. + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + const float_type ex = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x); + return (ex + (float_type(1) / ex)) / 2; +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TANH(boost::math::cstdfloat::detail::float_internal128_t x) +{ + // Patch the tanhq() function for a subset of broken GCC compilers + // like GCC 4.7, 4.8 on MinGW. + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + const float_type ex_plus = ::BOOST_CSTDFLOAT_FLOAT128_EXP(x); + const float_type ex_minus = (float_type(1) / ex_plus); + return (ex_plus - ex_minus) / (ex_plus + ex_minus); +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ASINH(boost::math::cstdfloat::detail::float_internal128_t x) throw() +{ + // Patch the asinh() function since quadmath does not have it. + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + return ::BOOST_CSTDFLOAT_FLOAT128_LOG(x + ::BOOST_CSTDFLOAT_FLOAT128_SQRT((x * x) + float_type(1))); +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ACOSH(boost::math::cstdfloat::detail::float_internal128_t x) throw() +{ + // Patch the acosh() function since quadmath does not have it. + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + const float_type zp(x + float_type(1)); + const float_type zm(x - float_type(1)); + + return ::BOOST_CSTDFLOAT_FLOAT128_LOG(x + (zp * ::BOOST_CSTDFLOAT_FLOAT128_SQRT(zm / zp))); +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_ATANH(boost::math::cstdfloat::detail::float_internal128_t x) throw() +{ + // Patch the atanh() function since quadmath does not have it. + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + return (::BOOST_CSTDFLOAT_FLOAT128_LOG(float_type(1) + x) + - ::BOOST_CSTDFLOAT_FLOAT128_LOG(float_type(1) - x)) / 2; +} +inline boost::math::cstdfloat::detail::float_internal128_t BOOST_CSTDFLOAT_FLOAT128_TGAMMA(boost::math::cstdfloat::detail::float_internal128_t x) throw() +{ + // Patch the tgammaq() function for a subset of broken GCC compilers + // like GCC 4.7, 4.8 on MinGW. + typedef boost::math::cstdfloat::detail::float_internal128_t float_type; + + if (x > float_type(0)) + { return ::BOOST_CSTDFLOAT_FLOAT128_EXP(::BOOST_CSTDFLOAT_FLOAT128_LGAMMA(x)); - } - else if(x < float_type(0)) - { + } + else if (x < float_type(0)) + { // For x < 0, compute tgamma(-x) and use the reflection formula. - const float_type positive_x = -x; - float_type gamma_value = ::BOOST_CSTDFLOAT_FLOAT128_TGAMMA(positive_x); - const float_type floor_of_positive_x = ::BOOST_CSTDFLOAT_FLOAT128_FLOOR (positive_x); + const float_type positive_x = -x; + float_type gamma_value = ::BOOST_CSTDFLOAT_FLOAT128_TGAMMA(positive_x); + const float_type floor_of_positive_x = ::BOOST_CSTDFLOAT_FLOAT128_FLOOR(positive_x); // Take the reflection checks (slightly adapted) from . const bool floor_of_z_is_equal_to_z = (positive_x == ::BOOST_CSTDFLOAT_FLOAT128_FLOOR(positive_x)); BOOST_CONSTEXPR_OR_CONST float_type my_pi = BOOST_FLOAT128_C(3.14159265358979323846264338327950288419716939937511); - if(floor_of_z_is_equal_to_z) + if (floor_of_z_is_equal_to_z) { - const bool is_odd = ((boost::int32_t(floor_of_positive_x) % boost::int32_t(2)) != boost::int32_t(0)); + const bool is_odd = ((boost::int32_t(floor_of_positive_x) % boost::int32_t(2)) != boost::int32_t(0)); - return (is_odd ? -std::numeric_limits::infinity() - : +std::numeric_limits::infinity()); + return (is_odd ? -std::numeric_limits::infinity() + : +std::numeric_limits::infinity()); } const float_type sinpx_value = x * ::BOOST_CSTDFLOAT_FLOAT128_SIN(my_pi * x); gamma_value *= sinpx_value; - const bool result_is_too_large_to_represent = ( (::BOOST_CSTDFLOAT_FLOAT128_FABS(gamma_value) < float_type(1)) - && (((std::numeric_limits::max)() * ::BOOST_CSTDFLOAT_FLOAT128_FABS(gamma_value)) < my_pi)); + const bool result_is_too_large_to_represent = ((::BOOST_CSTDFLOAT_FLOAT128_FABS(gamma_value) < float_type(1)) + && (((std::numeric_limits::max)() * ::BOOST_CSTDFLOAT_FLOAT128_FABS(gamma_value)) < my_pi)); - if(result_is_too_large_to_represent) + if (result_is_too_large_to_represent) { - const bool is_odd = ((boost::int32_t(floor_of_positive_x) % boost::int32_t(2)) != boost::int32_t(0)); + const bool is_odd = ((boost::int32_t(floor_of_positive_x) % boost::int32_t(2)) != boost::int32_t(0)); - return (is_odd ? -std::numeric_limits::infinity() - : +std::numeric_limits::infinity()); + return (is_odd ? -std::numeric_limits::infinity() + : +std::numeric_limits::infinity()); } gamma_value = -my_pi / gamma_value; - if((gamma_value > float_type(0)) || (gamma_value < float_type(0))) + if ((gamma_value > float_type(0)) || (gamma_value < float_type(0))) { - return gamma_value; + return gamma_value; } else { - // The value of gamma is too small to represent. Return 0.0 here. - return float_type(0); + // The value of gamma is too small to represent. Return 0.0 here. + return float_type(0); } - } - else - { + } + else + { // Gamma of zero is complex infinity. Return NaN here. return std::numeric_limits::quiet_NaN(); - } - } - #endif // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS - - // Define the quadruple-precision functions in the namespace boost::math::cstdfloat::detail. - - namespace boost { namespace math { namespace cstdfloat { namespace detail { - 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); } - 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); } - inline boost::math::cstdfloat::detail::float_internal128_t fabs (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FABS (x); } - inline boost::math::cstdfloat::detail::float_internal128_t abs (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FABS (x); } - inline boost::math::cstdfloat::detail::float_internal128_t floor (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FLOOR (x); } - inline boost::math::cstdfloat::detail::float_internal128_t ceil (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_CEIL (x); } - inline boost::math::cstdfloat::detail::float_internal128_t sqrt (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SQRT (x); } - inline boost::math::cstdfloat::detail::float_internal128_t trunc (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TRUNC (x); } - inline boost::math::cstdfloat::detail::float_internal128_t exp (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_EXP (x); } - 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); } - 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)); } - inline boost::math::cstdfloat::detail::float_internal128_t log (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG (x); } - inline boost::math::cstdfloat::detail::float_internal128_t log10 (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG10 (x); } - inline boost::math::cstdfloat::detail::float_internal128_t sin (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SIN (x); } - inline boost::math::cstdfloat::detail::float_internal128_t cos (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_COS (x); } - inline boost::math::cstdfloat::detail::float_internal128_t tan (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TAN (x); } - inline boost::math::cstdfloat::detail::float_internal128_t asin (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ASIN (x); } - inline boost::math::cstdfloat::detail::float_internal128_t acos (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ACOS (x); } - inline boost::math::cstdfloat::detail::float_internal128_t atan (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ATAN (x); } - inline boost::math::cstdfloat::detail::float_internal128_t sinh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SINH (x); } - inline boost::math::cstdfloat::detail::float_internal128_t cosh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_COSH (x); } - inline boost::math::cstdfloat::detail::float_internal128_t tanh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TANH (x); } - inline boost::math::cstdfloat::detail::float_internal128_t asinh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ASINH (x); } - inline boost::math::cstdfloat::detail::float_internal128_t acosh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ACOSH (x); } - inline boost::math::cstdfloat::detail::float_internal128_t atanh (boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ATANH (x); } - 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); } - 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); } - inline boost::math::cstdfloat::detail::float_internal128_t lgamma(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LGAMMA(x); } - inline boost::math::cstdfloat::detail::float_internal128_t tgamma(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TGAMMA(x); } - } } } } // boost::math::cstdfloat::detail - - // We will now inject the quadruple-precision functions - // into the std namespace. This is done via *using* directive. - namespace std - { - using boost::math::cstdfloat::detail::ldexp; - using boost::math::cstdfloat::detail::frexp; - using boost::math::cstdfloat::detail::fabs; + } +} +#endif // BOOST_CSTDFLOAT_BROKEN_FLOAT128_MATH_FUNCTIONS + +// Define the quadruple-precision functions in the namespace boost::math::cstdfloat::detail. + +namespace boost { + namespace math { + namespace cstdfloat { + namespace detail { + 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); } + 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); } + inline boost::math::cstdfloat::detail::float_internal128_t fabs(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FABS(x); } + inline boost::math::cstdfloat::detail::float_internal128_t abs(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FABS(x); } + inline boost::math::cstdfloat::detail::float_internal128_t floor(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_FLOOR(x); } + inline boost::math::cstdfloat::detail::float_internal128_t ceil(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_CEIL(x); } + inline boost::math::cstdfloat::detail::float_internal128_t sqrt(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SQRT(x); } + inline boost::math::cstdfloat::detail::float_internal128_t trunc(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TRUNC(x); } + inline boost::math::cstdfloat::detail::float_internal128_t exp(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_EXP(x); } + inline boost::math::cstdfloat::detail::float_internal128_t expm1(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_EXPM1(x); } + 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); } + 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)); } + inline boost::math::cstdfloat::detail::float_internal128_t log(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG(x); } + inline boost::math::cstdfloat::detail::float_internal128_t log10(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG10(x); } + inline boost::math::cstdfloat::detail::float_internal128_t sin(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SIN(x); } + inline boost::math::cstdfloat::detail::float_internal128_t cos(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_COS(x); } + inline boost::math::cstdfloat::detail::float_internal128_t tan(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TAN(x); } + inline boost::math::cstdfloat::detail::float_internal128_t asin(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ASIN(x); } + inline boost::math::cstdfloat::detail::float_internal128_t acos(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ACOS(x); } + inline boost::math::cstdfloat::detail::float_internal128_t atan(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ATAN(x); } + inline boost::math::cstdfloat::detail::float_internal128_t sinh(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SINH(x); } + inline boost::math::cstdfloat::detail::float_internal128_t cosh(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_COSH(x); } + inline boost::math::cstdfloat::detail::float_internal128_t tanh(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TANH(x); } + inline boost::math::cstdfloat::detail::float_internal128_t asinh(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ASINH(x); } + inline boost::math::cstdfloat::detail::float_internal128_t acosh(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ACOSH(x); } + inline boost::math::cstdfloat::detail::float_internal128_t atanh(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ATANH(x); } + 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); } + 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); } + inline boost::math::cstdfloat::detail::float_internal128_t lgamma(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LGAMMA(x); } + inline boost::math::cstdfloat::detail::float_internal128_t tgamma(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_TGAMMA(x); } + // begin more functions + inline boost::math::cstdfloat::detail::float_internal128_t remainder(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_REMAINDER(x, y); } + inline boost::math::cstdfloat::detail::float_internal128_t remquo(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y, int* z) { return ::BOOST_CSTDFLOAT_FLOAT128_REMQUO(x, y, z); } + inline boost::math::cstdfloat::detail::float_internal128_t fma(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y, boost::math::cstdfloat::detail::float_internal128_t z) { return BOOST_CSTDFLOAT_FLOAT128_FMA(x, y, z); } + + inline boost::math::cstdfloat::detail::float_internal128_t fmax(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_FMAX(x, y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + fmax(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return ::BOOST_CSTDFLOAT_FLOAT128_FMAX(x, y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + fmax(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_FMAX(x, y); } + inline boost::math::cstdfloat::detail::float_internal128_t fmin(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_FMIN(x, y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + fmin(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return ::BOOST_CSTDFLOAT_FLOAT128_FMIN(x, y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + fmin(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_FMIN(x, y); } + + inline boost::math::cstdfloat::detail::float_internal128_t fdim(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_FDIM(x, y); } + inline boost::math::cstdfloat::detail::float_internal128_t nanq(const char* x) { return ::BOOST_CSTDFLOAT_FLOAT128_NAN(x); } + inline boost::math::cstdfloat::detail::float_internal128_t exp2(boost::math::cstdfloat::detail::float_internal128_t x) + { + return ::BOOST_CSTDFLOAT_FLOAT128_POW(boost::math::cstdfloat::detail::float_internal128_t(2), x); + } + inline boost::math::cstdfloat::detail::float_internal128_t log2(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG2(x); } + inline boost::math::cstdfloat::detail::float_internal128_t log1p(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOG1P(x); } + inline boost::math::cstdfloat::detail::float_internal128_t cbrt(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_CBRT(x); } + inline boost::math::cstdfloat::detail::float_internal128_t hypot(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y, boost::math::cstdfloat::detail::float_internal128_t z) { return ::BOOST_CSTDFLOAT_FLOAT128_SQRT(x*x + y * y + z * z); } + inline boost::math::cstdfloat::detail::float_internal128_t hypot(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_HYPOT(x, y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + hypot(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return ::BOOST_CSTDFLOAT_FLOAT128_HYPOT(x, y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + hypot(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_HYPOT(x, y); } + + + inline boost::math::cstdfloat::detail::float_internal128_t erf(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ERF(x); } + inline boost::math::cstdfloat::detail::float_internal128_t erfc(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ERFC(x); } + inline long long int llround(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LLROUND(x); } + inline long int lround(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LROUND(x); } + inline boost::math::cstdfloat::detail::float_internal128_t round(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ROUND(x); } + inline boost::math::cstdfloat::detail::float_internal128_t nearbyint(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_NEARBYINT(x); } + inline long long int llrint(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LLRINT(x); } + inline long int lrint(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LRINT(x); } + inline boost::math::cstdfloat::detail::float_internal128_t rint(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_RINT(x); } + inline boost::math::cstdfloat::detail::float_internal128_t modf(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t* y) { return ::BOOST_CSTDFLOAT_FLOAT128_MODF(x, y); } + inline boost::math::cstdfloat::detail::float_internal128_t scalbln(boost::math::cstdfloat::detail::float_internal128_t x, long int y) { return ::BOOST_CSTDFLOAT_FLOAT128_SCALBLN(x, y); } + inline boost::math::cstdfloat::detail::float_internal128_t scalbn(boost::math::cstdfloat::detail::float_internal128_t x, int y) { return ::BOOST_CSTDFLOAT_FLOAT128_SCALBN(x, y); } + inline int ilogb(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ILOGB(x); } + inline boost::math::cstdfloat::detail::float_internal128_t logb(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_LOGB(x); } + inline boost::math::cstdfloat::detail::float_internal128_t nextafter(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_NEXTAFTER(x, y); } + inline boost::math::cstdfloat::detail::float_internal128_t nexttoward(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return -(::BOOST_CSTDFLOAT_FLOAT128_NEXTAFTER(-x, -y)); } + inline boost::math::cstdfloat::detail::float_internal128_t copysign BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_COPYSIGN(x, y); } + inline bool signbit BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_SIGNBIT(x); } + inline int fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x) + { + if (::BOOST_CSTDFLOAT_FLOAT128_ISNAN(x)) + return FP_NAN; + else if (::BOOST_CSTDFLOAT_FLOAT128_ISINF(x)) + return FP_INFINITE; + else if (x == BOOST_FLOAT128_C(0.0)) + return FP_ZERO; + + if (::BOOST_CSTDFLOAT_FLOAT128_FABS(x) < BOOST_CSTDFLOAT_FLOAT128_MIN) + return FP_SUBNORMAL; + else + return FP_NORMAL; + } + inline bool isfinite BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x) + { + return !::BOOST_CSTDFLOAT_FLOAT128_ISNAN(x) && !::BOOST_CSTDFLOAT_FLOAT128_ISINF(x); + } + inline bool isinf BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ISINF(x); } + inline bool isnan BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x) { return ::BOOST_CSTDFLOAT_FLOAT128_ISNAN(x); } + inline bool isnormal BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x) { return boost::math::cstdfloat::detail::fpclassify BOOST_PREVENT_MACRO_SUBSTITUTION(x) == FP_NORMAL; } + inline bool isgreater BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) + { + if (isnan BOOST_PREVENT_MACRO_SUBSTITUTION(x) || isnan BOOST_PREVENT_MACRO_SUBSTITUTION(y)) + return false; + return x > y; + } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isgreater BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return isgreater BOOST_PREVENT_MACRO_SUBSTITUTION(x, (boost::math::cstdfloat::detail::float_internal128_t)y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isgreater BOOST_PREVENT_MACRO_SUBSTITUTION(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return isgreater BOOST_PREVENT_MACRO_SUBSTITUTION((boost::math::cstdfloat::detail::float_internal128_t)x, y); } + + inline bool isgreaterequal BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) + { + if (isnan BOOST_PREVENT_MACRO_SUBSTITUTION(x) || isnan BOOST_PREVENT_MACRO_SUBSTITUTION(y)) + return false; + return x >= y; + } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isgreaterequal BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return isgreaterequal BOOST_PREVENT_MACRO_SUBSTITUTION(x, (boost::math::cstdfloat::detail::float_internal128_t)y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isgreaterequal BOOST_PREVENT_MACRO_SUBSTITUTION(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return isgreaterequal BOOST_PREVENT_MACRO_SUBSTITUTION((boost::math::cstdfloat::detail::float_internal128_t)x, y); } + + inline bool isless BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) + { + if (isnan BOOST_PREVENT_MACRO_SUBSTITUTION(x) || isnan BOOST_PREVENT_MACRO_SUBSTITUTION(y)) + return false; + return x < y; + } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isless BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return isless BOOST_PREVENT_MACRO_SUBSTITUTION(x, (boost::math::cstdfloat::detail::float_internal128_t)y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isless BOOST_PREVENT_MACRO_SUBSTITUTION(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return isless BOOST_PREVENT_MACRO_SUBSTITUTION((boost::math::cstdfloat::detail::float_internal128_t)x, y); } + + + inline bool islessequal BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) + { + if (isnan BOOST_PREVENT_MACRO_SUBSTITUTION(x) || isnan BOOST_PREVENT_MACRO_SUBSTITUTION(y)) + return false; + return x <= y; + } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + islessequal BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return islessequal BOOST_PREVENT_MACRO_SUBSTITUTION(x, (boost::math::cstdfloat::detail::float_internal128_t)y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + islessequal BOOST_PREVENT_MACRO_SUBSTITUTION(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return islessequal BOOST_PREVENT_MACRO_SUBSTITUTION((boost::math::cstdfloat::detail::float_internal128_t)x, y); } + + + inline bool islessgreater BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) + { + if (isnan BOOST_PREVENT_MACRO_SUBSTITUTION(x) || isnan BOOST_PREVENT_MACRO_SUBSTITUTION(y)) + return false; + return (x < y) || (x > y); + } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + islessgreater BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return islessgreater BOOST_PREVENT_MACRO_SUBSTITUTION(x, (boost::math::cstdfloat::detail::float_internal128_t)y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + islessgreater BOOST_PREVENT_MACRO_SUBSTITUTION(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return islessgreater BOOST_PREVENT_MACRO_SUBSTITUTION((boost::math::cstdfloat::detail::float_internal128_t)x, y); } + + + inline bool isunordered BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, boost::math::cstdfloat::detail::float_internal128_t y) { return ::BOOST_CSTDFLOAT_FLOAT128_ISNAN(x) || ::BOOST_CSTDFLOAT_FLOAT128_ISNAN(y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isunordered BOOST_PREVENT_MACRO_SUBSTITUTION(boost::math::cstdfloat::detail::float_internal128_t x, T y) { return isunordered BOOST_PREVENT_MACRO_SUBSTITUTION(x, (boost::math::cstdfloat::detail::float_internal128_t)y); } + template + inline typename boost::enable_if_c< + boost::is_convertible::value + && !boost::is_same::value, boost::math::cstdfloat::detail::float_internal128_t>::type + isunordered BOOST_PREVENT_MACRO_SUBSTITUTION(T x, boost::math::cstdfloat::detail::float_internal128_t y) { return isunordered BOOST_PREVENT_MACRO_SUBSTITUTION((boost::math::cstdfloat::detail::float_internal128_t)x, y); } + + + // end more functions + } + } + } +} // boost::math::cstdfloat::detail + +// We will now inject the quadruple-precision functions +// into the std namespace. This is done via *using* directive. +namespace std +{ + using boost::math::cstdfloat::detail::ldexp; + using boost::math::cstdfloat::detail::frexp; + using boost::math::cstdfloat::detail::fabs; #if !(defined(_GLIBCXX_USE_FLOAT128) && defined(__GNUC__) && (__GNUC__ >= 7)) - using boost::math::cstdfloat::detail::abs; + using boost::math::cstdfloat::detail::abs; #endif - using boost::math::cstdfloat::detail::floor; - using boost::math::cstdfloat::detail::ceil; - using boost::math::cstdfloat::detail::sqrt; - using boost::math::cstdfloat::detail::trunc; - using boost::math::cstdfloat::detail::exp; - using boost::math::cstdfloat::detail::pow; - using boost::math::cstdfloat::detail::log; - using boost::math::cstdfloat::detail::log10; - using boost::math::cstdfloat::detail::sin; - using boost::math::cstdfloat::detail::cos; - using boost::math::cstdfloat::detail::tan; - using boost::math::cstdfloat::detail::asin; - using boost::math::cstdfloat::detail::acos; - using boost::math::cstdfloat::detail::atan; - using boost::math::cstdfloat::detail::sinh; - using boost::math::cstdfloat::detail::cosh; - using boost::math::cstdfloat::detail::tanh; - using boost::math::cstdfloat::detail::asinh; - using boost::math::cstdfloat::detail::acosh; - using boost::math::cstdfloat::detail::atanh; - using boost::math::cstdfloat::detail::fmod; - using boost::math::cstdfloat::detail::atan2; - using boost::math::cstdfloat::detail::lgamma; - using boost::math::cstdfloat::detail::tgamma; - } // namespace std - - // We will now remove the preprocessor symbols representing quadruple-precision - // functions from the preprocessor. - - #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_TRUNC - #undef BOOST_CSTDFLOAT_FLOAT128_EXP - #undef BOOST_CSTDFLOAT_FLOAT128_EXPM1 - #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_ASINH - #undef BOOST_CSTDFLOAT_FLOAT128_ACOSH - #undef BOOST_CSTDFLOAT_FLOAT128_ATANH - #undef BOOST_CSTDFLOAT_FLOAT128_FMOD - #undef BOOST_CSTDFLOAT_FLOAT128_ATAN2 - #undef BOOST_CSTDFLOAT_FLOAT128_LGAMMA - #undef BOOST_CSTDFLOAT_FLOAT128_TGAMMA - - #endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT (i.e., the user would like to have libquadmath support) + using boost::math::cstdfloat::detail::floor; + using boost::math::cstdfloat::detail::ceil; + using boost::math::cstdfloat::detail::sqrt; + using boost::math::cstdfloat::detail::trunc; + using boost::math::cstdfloat::detail::exp; + using boost::math::cstdfloat::detail::expm1; + using boost::math::cstdfloat::detail::pow; + using boost::math::cstdfloat::detail::log; + using boost::math::cstdfloat::detail::log10; + using boost::math::cstdfloat::detail::sin; + using boost::math::cstdfloat::detail::cos; + using boost::math::cstdfloat::detail::tan; + using boost::math::cstdfloat::detail::asin; + using boost::math::cstdfloat::detail::acos; + using boost::math::cstdfloat::detail::atan; + using boost::math::cstdfloat::detail::sinh; + using boost::math::cstdfloat::detail::cosh; + using boost::math::cstdfloat::detail::tanh; + using boost::math::cstdfloat::detail::asinh; + using boost::math::cstdfloat::detail::acosh; + using boost::math::cstdfloat::detail::atanh; + using boost::math::cstdfloat::detail::fmod; + using boost::math::cstdfloat::detail::atan2; + using boost::math::cstdfloat::detail::lgamma; + using boost::math::cstdfloat::detail::tgamma; + + // begin more functions + using boost::math::cstdfloat::detail::remainder; + using boost::math::cstdfloat::detail::remquo; + using boost::math::cstdfloat::detail::fma; + using boost::math::cstdfloat::detail::fmax; + using boost::math::cstdfloat::detail::fmin; + using boost::math::cstdfloat::detail::fdim; + using boost::math::cstdfloat::detail::nanq; + using boost::math::cstdfloat::detail::exp2; + using boost::math::cstdfloat::detail::log2; + using boost::math::cstdfloat::detail::log1p; + using boost::math::cstdfloat::detail::cbrt; + using boost::math::cstdfloat::detail::hypot; + using boost::math::cstdfloat::detail::erf; + using boost::math::cstdfloat::detail::erfc; + using boost::math::cstdfloat::detail::llround; + using boost::math::cstdfloat::detail::lround; + using boost::math::cstdfloat::detail::round; + using boost::math::cstdfloat::detail::nearbyint; + using boost::math::cstdfloat::detail::llrint; + using boost::math::cstdfloat::detail::lrint; + using boost::math::cstdfloat::detail::rint; + using boost::math::cstdfloat::detail::modf; + using boost::math::cstdfloat::detail::scalbln; + using boost::math::cstdfloat::detail::scalbn; + using boost::math::cstdfloat::detail::ilogb; + using boost::math::cstdfloat::detail::logb; + using boost::math::cstdfloat::detail::nextafter; + using boost::math::cstdfloat::detail::nexttoward; + using boost::math::cstdfloat::detail::copysign; + using boost::math::cstdfloat::detail::signbit; + using boost::math::cstdfloat::detail::fpclassify; + using boost::math::cstdfloat::detail::isfinite; + using boost::math::cstdfloat::detail::isinf; + using boost::math::cstdfloat::detail::isnan; + using boost::math::cstdfloat::detail::isnormal; + using boost::math::cstdfloat::detail::isgreater; + using boost::math::cstdfloat::detail::isgreaterequal; + using boost::math::cstdfloat::detail::isless; + using boost::math::cstdfloat::detail::islessequal; + using boost::math::cstdfloat::detail::islessgreater; + using boost::math::cstdfloat::detail::isunordered; + // end more functions + + // + // Very basic iostream operator: + // + inline std::ostream& operator << (std::ostream& os, __float128 m_value) + { + std::streamsize digits = os.precision(); + std::ios_base::fmtflags f = os.flags(); + std::string s; + + char buf[100]; + boost::scoped_array buf2; + std::string format = "%"; + if (f & std::ios_base::showpos) + format += "+"; + if (f & std::ios_base::showpoint) + format += "#"; + format += ".*"; + if (digits == 0) + digits = 36; + format += "Q"; + if (f & std::ios_base::scientific) + format += "e"; + else if (f & std::ios_base::fixed) + format += "f"; + else + format += "g"; + + int v = quadmath_snprintf(buf, 100, format.c_str(), digits, m_value); + + if ((v < 0) || (v >= 99)) + { + int v_max = v; + buf2.reset(new char[v + 3]); + v = quadmath_snprintf(&buf2[0], v_max + 3, format.c_str(), digits, m_value); + if (v >= v_max + 3) + { + BOOST_THROW_EXCEPTION(std::runtime_error("Formatting of float128_type failed.")); + } + s = &buf2[0]; + } + else + s = buf; + std::streamsize ss = os.width(); + if (ss > static_cast(s.size())) + { + char fill = os.fill(); + if ((os.flags() & std::ios_base::left) == std::ios_base::left) + s.append(static_cast(ss - s.size()), fill); + else + s.insert(static_cast(0), static_cast(ss - s.size()), fill); + } + + return os << s; + } + + +} // namespace std + +// We will now remove the preprocessor symbols representing quadruple-precision +// functions from the preprocessor. + +#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_TRUNC +#undef BOOST_CSTDFLOAT_FLOAT128_EXP +#undef BOOST_CSTDFLOAT_FLOAT128_EXPM1 +#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_ASINH +#undef BOOST_CSTDFLOAT_FLOAT128_ACOSH +#undef BOOST_CSTDFLOAT_FLOAT128_ATANH +#undef BOOST_CSTDFLOAT_FLOAT128_FMOD +#undef BOOST_CSTDFLOAT_FLOAT128_ATAN2 +#undef BOOST_CSTDFLOAT_FLOAT128_LGAMMA +#undef BOOST_CSTDFLOAT_FLOAT128_TGAMMA + +// begin more functions +#undef BOOST_CSTDFLOAT_FLOAT128_REMAINDER +#undef BOOST_CSTDFLOAT_FLOAT128_REMQUO +#undef BOOST_CSTDFLOAT_FLOAT128_FMA +#undef BOOST_CSTDFLOAT_FLOAT128_FMAX +#undef BOOST_CSTDFLOAT_FLOAT128_FMIN +#undef BOOST_CSTDFLOAT_FLOAT128_FDIM +#undef BOOST_CSTDFLOAT_FLOAT128_NAN +#undef BOOST_CSTDFLOAT_FLOAT128_EXP2 +#undef BOOST_CSTDFLOAT_FLOAT128_LOG2 +#undef BOOST_CSTDFLOAT_FLOAT128_LOG1P +#undef BOOST_CSTDFLOAT_FLOAT128_CBRT +#undef BOOST_CSTDFLOAT_FLOAT128_HYPOT +#undef BOOST_CSTDFLOAT_FLOAT128_ERF +#undef BOOST_CSTDFLOAT_FLOAT128_ERFC +#undef BOOST_CSTDFLOAT_FLOAT128_LLROUND +#undef BOOST_CSTDFLOAT_FLOAT128_LROUND +#undef BOOST_CSTDFLOAT_FLOAT128_ROUND +#undef BOOST_CSTDFLOAT_FLOAT128_NEARBYINT +#undef BOOST_CSTDFLOAT_FLOAT128_LLRINT +#undef BOOST_CSTDFLOAT_FLOAT128_LRINT +#undef BOOST_CSTDFLOAT_FLOAT128_RINT +#undef BOOST_CSTDFLOAT_FLOAT128_MODF +#undef BOOST_CSTDFLOAT_FLOAT128_SCALBLN +#undef BOOST_CSTDFLOAT_FLOAT128_SCALBN +#undef BOOST_CSTDFLOAT_FLOAT128_ILOGB +#undef BOOST_CSTDFLOAT_FLOAT128_LOGB +#undef BOOST_CSTDFLOAT_FLOAT128_NEXTAFTER +#undef BOOST_CSTDFLOAT_FLOAT128_NEXTTOWARD +#undef BOOST_CSTDFLOAT_FLOAT128_COPYSIGN +#undef BOOST_CSTDFLOAT_FLOAT128_SIGNBIT +#undef BOOST_CSTDFLOAT_FLOAT128_FPCLASSIFY +#undef BOOST_CSTDFLOAT_FLOAT128_ISFINITE +#undef BOOST_CSTDFLOAT_FLOAT128_ISINF +#undef BOOST_CSTDFLOAT_FLOAT128_ISNAN +#undef BOOST_CSTDFLOAT_FLOAT128_ISNORMAL +#undef BOOST_CSTDFLOAT_FLOAT128_ISGREATER +#undef BOOST_CSTDFLOAT_FLOAT128_ISGREATEREQUAL +#undef BOOST_CSTDFLOAT_FLOAT128_ISLESS +#undef BOOST_CSTDFLOAT_FLOAT128_ISLESSEQUAL +#undef BOOST_CSTDFLOAT_FLOAT128_ISLESSGREATER +#undef BOOST_CSTDFLOAT_FLOAT128_ISUNORDERED +// end more functions + +#endif // Not BOOST_CSTDFLOAT_NO_LIBQUADMATH_SUPPORT (i.e., the user would like to have libquadmath support) #endif // _BOOST_CSTDFLOAT_CMATH_2014_02_15_HPP_ + diff --git a/boost/math/differentiation/autodiff.hpp b/boost/math/differentiation/autodiff.hpp new file mode 100644 index 0000000000..e98eecab6d --- /dev/null +++ b/boost/math/differentiation/autodiff.hpp @@ -0,0 +1,2053 @@ +// Copyright Matthew Pulver 2018 - 2019. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP +#define BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace math { +namespace differentiation { +// Automatic Differentiation v1 +inline namespace autodiff_v1 { +namespace detail { + +template +struct promote_args_n { + using type = typename tools::promote_args_2::type>::type; +}; + +template +struct promote_args_n { + using type = typename tools::promote_arg::type; +}; + +} // namespace detail + +template +using promote = typename detail::promote_args_n::type; + +namespace detail { + +template +class fvar; + +template +struct is_fvar_impl : std::false_type {}; + +template +struct is_fvar_impl> : std::true_type {}; + +template +using is_fvar = is_fvar_impl>; + +template +struct nest_fvar { + using type = fvar::type, Order>; +}; + +template +struct nest_fvar { + using type = fvar; +}; + +template +struct get_depth_impl : std::integral_constant {}; + +template +struct get_depth_impl> + : std::integral_constant::value + 1> {}; + +template +using get_depth = get_depth_impl>; + +template +struct get_order_sum_t : std::integral_constant {}; + +template +struct get_order_sum_t> + : std::integral_constant::value + Order> {}; + +template +using get_order_sum = get_order_sum_t>; + +template +struct get_root_type { + using type = RealType; +}; + +template +struct get_root_type> { + using type = typename get_root_type::type; +}; + +template +struct type_at { + using type = RealType; +}; + +template +struct type_at, Depth> { + using type = typename conditional, + typename type_at::type>::type; +}; + +template +using get_type_at = typename type_at::type; + +// Satisfies Boost's Conceptual Requirements for Real Number Types. +// https://www.boost.org/libs/math/doc/html/math_toolkit/real_concepts.html +template +class fvar { + std::array v; + + public: + using root_type = typename get_root_type::type; // RealType in the root fvar. + + fvar() = default; + + // Initialize a variable or constant. + fvar(root_type const&, bool const is_variable); + + // RealType(cr) | RealType | RealType is copy constructible. + fvar(fvar const&) = default; + + // Be aware of implicit casting from one fvar<> type to another by this copy constructor. + template + fvar(fvar const&); + + // RealType(ca) | RealType | RealType is copy constructible from the arithmetic types. + explicit fvar(root_type const&); // Initialize a constant. (No epsilon terms.) + + template + fvar(RealType2 const& ca); // Supports any RealType2 for which static_cast(ca) compiles. + + // r = cr | RealType& | Assignment operator. + fvar& operator=(fvar const&) = default; + + // r = ca | RealType& | Assignment operator from the arithmetic types. + // Handled by constructor that takes a single parameter of generic type. + // fvar& operator=(root_type const&); // Set a constant. + + // r += cr | RealType& | Adds cr to r. + template + fvar& operator+=(fvar const&); + + // r += ca | RealType& | Adds ar to r. + fvar& operator+=(root_type const&); + + // r -= cr | RealType& | Subtracts cr from r. + template + fvar& operator-=(fvar const&); + + // r -= ca | RealType& | Subtracts ca from r. + fvar& operator-=(root_type const&); + + // r *= cr | RealType& | Multiplies r by cr. + template + fvar& operator*=(fvar const&); + + // r *= ca | RealType& | Multiplies r by ca. + fvar& operator*=(root_type const&); + + // r /= cr | RealType& | Divides r by cr. + template + fvar& operator/=(fvar const&); + + // r /= ca | RealType& | Divides r by ca. + fvar& operator/=(root_type const&); + + // -r | RealType | Unary Negation. + fvar operator-() const; + + // +r | RealType& | Identity Operation. + fvar const& operator+() const; + + // cr + cr2 | RealType | Binary Addition + template + promote> operator+(fvar const&) const; + + // cr + ca | RealType | Binary Addition + fvar operator+(root_type const&) const; + + // ca + cr | RealType | Binary Addition + template + friend fvar operator+(typename fvar::root_type const&, + fvar const&); + + // cr - cr2 | RealType | Binary Subtraction + template + promote> operator-(fvar const&) const; + + // cr - ca | RealType | Binary Subtraction + fvar operator-(root_type const&) const; + + // ca - cr | RealType | Binary Subtraction + template + friend fvar operator-(typename fvar::root_type const&, + fvar const&); + + // cr * cr2 | RealType | Binary Multiplication + template + promote> operator*(fvar const&)const; + + // cr * ca | RealType | Binary Multiplication + fvar operator*(root_type const&)const; + + // ca * cr | RealType | Binary Multiplication + template + friend fvar operator*(typename fvar::root_type const&, + fvar const&); + + // cr / cr2 | RealType | Binary Subtraction + template + promote> operator/(fvar const&) const; + + // cr / ca | RealType | Binary Subtraction + fvar operator/(root_type const&) const; + + // ca / cr | RealType | Binary Subtraction + template + friend fvar operator/(typename fvar::root_type const&, + fvar const&); + + // For all comparison overloads, only the root term is compared. + + // cr == cr2 | bool | Equality Comparison + template + bool operator==(fvar const&) const; + + // cr == ca | bool | Equality Comparison + bool operator==(root_type const&) const; + + // ca == cr | bool | Equality Comparison + template + friend bool operator==(typename fvar::root_type const&, fvar const&); + + // cr != cr2 | bool | Inequality Comparison + template + bool operator!=(fvar const&) const; + + // cr != ca | bool | Inequality Comparison + bool operator!=(root_type const&) const; + + // ca != cr | bool | Inequality Comparison + template + friend bool operator!=(typename fvar::root_type const&, fvar const&); + + // cr <= cr2 | bool | Less than equal to. + template + bool operator<=(fvar const&) const; + + // cr <= ca | bool | Less than equal to. + bool operator<=(root_type const&) const; + + // ca <= cr | bool | Less than equal to. + template + friend bool operator<=(typename fvar::root_type const&, fvar const&); + + // cr >= cr2 | bool | Greater than equal to. + template + bool operator>=(fvar const&) const; + + // cr >= ca | bool | Greater than equal to. + bool operator>=(root_type const&) const; + + // ca >= cr | bool | Greater than equal to. + template + friend bool operator>=(typename fvar::root_type const&, fvar const&); + + // cr < cr2 | bool | Less than comparison. + template + bool operator<(fvar const&) const; + + // cr < ca | bool | Less than comparison. + bool operator<(root_type const&) const; + + // ca < cr | bool | Less than comparison. + template + friend bool operator<(typename fvar::root_type const&, fvar const&); + + // cr > cr2 | bool | Greater than comparison. + template + bool operator>(fvar const&) const; + + // cr > ca | bool | Greater than comparison. + bool operator>(root_type const&) const; + + // ca > cr | bool | Greater than comparison. + template + friend bool operator>(typename fvar::root_type const&, fvar const&); + + // Will throw std::out_of_range if Order < order. + template + get_type_at at(size_t order, Orders... orders) const; + + template + get_type_at derivative(Orders... orders) const; + + const RealType& operator[](size_t) const; + + fvar inverse() const; // Multiplicative inverse. + + fvar& negate(); // Negate and return reference to *this. + + static constexpr size_t depth = get_depth::value; // Number of nested std::array. + + static constexpr size_t order_sum = get_order_sum::value; + + explicit operator root_type() const; // Must be explicit, otherwise overloaded operators are ambiguous. + + template >>::type> + explicit operator T() const; // Must be explicit; multiprecision has trouble without the std::enable_if + + fvar& set_root(root_type const&); + + // Apply coefficients using horner method. + template + promote, Fvar, Fvars...> apply_coefficients(size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const; + + template + fvar apply_coefficients(size_t const order, Func const& f) const; + + // Use when function returns derivative(i)/factorial(i) and may have some infinite derivatives. + template + promote, Fvar, Fvars...> apply_coefficients_nonhorner(size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const; + + template + fvar apply_coefficients_nonhorner(size_t const order, Func const& f) const; + + // Apply derivatives using horner method. + template + promote, Fvar, Fvars...> apply_derivatives(size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const; + + template + fvar apply_derivatives(size_t const order, Func const& f) const; + + // Use when function returns derivative(i) and may have some infinite derivatives. + template + promote, Fvar, Fvars...> apply_derivatives_nonhorner(size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const; + + template + fvar apply_derivatives_nonhorner(size_t const order, Func const& f) const; + + private: + RealType epsilon_inner_product(size_t z0, + size_t isum0, + size_t m0, + fvar const& cr, + size_t z1, + size_t isum1, + size_t m1, + size_t j) const; + + fvar epsilon_multiply(size_t z0, size_t isum0, fvar const& cr, size_t z1, size_t isum1) const; + + fvar epsilon_multiply(size_t z0, size_t isum0, root_type const& ca) const; + + fvar inverse_apply() const; + + fvar& multiply_assign_by_root_type(bool is_root, root_type const&); + + template + friend class fvar; + + template + friend std::ostream& operator<<(std::ostream&, fvar const&); + + // C++11 Compatibility +#ifdef BOOST_NO_CXX17_IF_CONSTEXPR + template + void fvar_cpp11(std::true_type, RootType const& ca, bool const is_variable); + + template + void fvar_cpp11(std::false_type, RootType const& ca, bool const is_variable); + + template + get_type_at at_cpp11(std::true_type, size_t order, Orders... orders) const; + + template + get_type_at at_cpp11(std::false_type, size_t order, Orders... orders) const; + + template + fvar epsilon_multiply_cpp11(std::true_type, + SizeType z0, + size_t isum0, + fvar const& cr, + size_t z1, + size_t isum1) const; + + template + fvar epsilon_multiply_cpp11(std::false_type, + SizeType z0, + size_t isum0, + fvar const& cr, + size_t z1, + size_t isum1) const; + + template + fvar epsilon_multiply_cpp11(std::true_type, SizeType z0, size_t isum0, root_type const& ca) const; + + template + fvar epsilon_multiply_cpp11(std::false_type, SizeType z0, size_t isum0, root_type const& ca) const; + + template + fvar& multiply_assign_by_root_type_cpp11(std::true_type, bool is_root, RootType const& ca); + + template + fvar& multiply_assign_by_root_type_cpp11(std::false_type, bool is_root, RootType const& ca); + + template + fvar& negate_cpp11(std::true_type, RootType const&); + + template + fvar& negate_cpp11(std::false_type, RootType const&); + + template + fvar& set_root_cpp11(std::true_type, RootType const& root); + + template + fvar& set_root_cpp11(std::false_type, RootType const& root); +#endif +}; + +// C++11 compatibility +#ifdef BOOST_NO_CXX17_IF_CONSTEXPR +#define BOOST_AUTODIFF_IF_CONSTEXPR +#else +#define BOOST_AUTODIFF_IF_CONSTEXPR constexpr +#endif + +// Standard Library Support Requirements + +// fabs(cr1) | RealType +template +fvar fabs(fvar const&); + +// abs(cr1) | RealType +template +fvar abs(fvar const&); + +// ceil(cr1) | RealType +template +fvar ceil(fvar const&); + +// floor(cr1) | RealType +template +fvar floor(fvar const&); + +// exp(cr1) | RealType +template +fvar exp(fvar const&); + +// pow(cr, ca) | RealType +template +fvar pow(fvar const&, typename fvar::root_type const&); + +// pow(ca, cr) | RealType +template +fvar pow(typename fvar::root_type const&, fvar const&); + +// pow(cr1, cr2) | RealType +template +promote, fvar> pow(fvar const&, + fvar const&); + +// sqrt(cr1) | RealType +template +fvar sqrt(fvar const&); + +// log(cr1) | RealType +template +fvar log(fvar const&); + +// frexp(cr1, &i) | RealType +template +fvar frexp(fvar const&, int*); + +// ldexp(cr1, i) | RealType +template +fvar ldexp(fvar const&, int); + +// cos(cr1) | RealType +template +fvar cos(fvar const&); + +// sin(cr1) | RealType +template +fvar sin(fvar const&); + +// asin(cr1) | RealType +template +fvar asin(fvar const&); + +// tan(cr1) | RealType +template +fvar tan(fvar const&); + +// atan(cr1) | RealType +template +fvar atan(fvar const&); + +// atan2(cr, ca) | RealType +template +fvar atan2(fvar const&, typename fvar::root_type const&); + +// atan2(ca, cr) | RealType +template +fvar atan2(typename fvar::root_type const&, fvar const&); + +// atan2(cr1, cr2) | RealType +template +promote, fvar> atan2(fvar const&, + fvar const&); + +// fmod(cr1,cr2) | RealType +template +promote, fvar> fmod(fvar const&, + fvar const&); + +// round(cr1) | RealType +template +fvar round(fvar const&); + +// iround(cr1) | int +template +int iround(fvar const&); + +template +long lround(fvar const&); + +template +long long llround(fvar const&); + +// trunc(cr1) | RealType +template +fvar trunc(fvar const&); + +template +long double truncl(fvar const&); + +// itrunc(cr1) | int +template +int itrunc(fvar const&); + +template +long long lltrunc(fvar const&); + +// Additional functions +template +fvar acos(fvar const&); + +template +fvar acosh(fvar const&); + +template +fvar asinh(fvar const&); + +template +fvar atanh(fvar const&); + +template +fvar cosh(fvar const&); + +template +fvar digamma(fvar const&); + +template +fvar erf(fvar const&); + +template +fvar erfc(fvar const&); + +template +fvar lambert_w0(fvar const&); + +template +fvar lgamma(fvar const&); + +template +fvar sinc(fvar const&); + +template +fvar sinh(fvar const&); + +template +fvar tanh(fvar const&); + +template +fvar tgamma(fvar const&); + +template +struct zero : std::integral_constant {}; + +} // namespace detail + +template +using autodiff_fvar = typename detail::nest_fvar::type; + +template +autodiff_fvar make_fvar(RealType const& ca) { + return autodiff_fvar(ca, true); +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +namespace detail { + +template +auto make_fvar_for_tuple(std::index_sequence, RealType const& ca) { + return make_fvar::value..., Order>(ca); +} + +template +auto make_ftuple_impl(std::index_sequence, RealTypes const&... ca) { + return std::make_tuple(make_fvar_for_tuple(std::make_index_sequence{}, ca)...); +} + +} // namespace detail + +template +auto make_ftuple(RealTypes const&... ca) { + static_assert(sizeof...(Orders) == sizeof...(RealTypes), + "Number of Orders must match number of function parameters."); + return detail::make_ftuple_impl(std::index_sequence_for{}, ca...); +} +#endif + +namespace detail { + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +template +fvar::fvar(root_type const& ca, bool const is_variable) { + if constexpr (is_fvar::value) { + v.front() = RealType(ca, is_variable); + if constexpr (0 < Order) + std::fill(v.begin() + 1, v.end(), static_cast(0)); + } else { + v.front() = ca; + if constexpr (0 < Order) + v[1] = static_cast(static_cast(is_variable)); + if constexpr (1 < Order) + std::fill(v.begin() + 2, v.end(), static_cast(0)); + } +} +#endif + +template +template +fvar::fvar(fvar const& cr) { + for (size_t i = 0; i <= (std::min)(Order, Order2); ++i) + v[i] = static_cast(cr.v[i]); + if BOOST_AUTODIFF_IF_CONSTEXPR (Order2 < Order) + std::fill(v.begin() + (Order2 + 1), v.end(), static_cast(0)); +} + +template +fvar::fvar(root_type const& ca) : v{{static_cast(ca)}} {} + +// Can cause compiler error if RealType2 cannot be cast to root_type. +template +template +fvar::fvar(RealType2 const& ca) : v{{static_cast(ca)}} {} + +/* +template +fvar& fvar::operator=(root_type const& ca) +{ + v.front() = static_cast(ca); + if constexpr (0 < Order) + std::fill(v.begin()+1, v.end(), static_cast(0)); + return *this; +} +*/ + +template +template +fvar& fvar::operator+=(fvar const& cr) { + for (size_t i = 0; i <= (std::min)(Order, Order2); ++i) + v[i] += cr.v[i]; + return *this; +} + +template +fvar& fvar::operator+=(root_type const& ca) { + v.front() += ca; + return *this; +} + +template +template +fvar& fvar::operator-=(fvar const& cr) { + for (size_t i = 0; i <= Order; ++i) + v[i] -= cr.v[i]; + return *this; +} + +template +fvar& fvar::operator-=(root_type const& ca) { + v.front() -= ca; + return *this; +} + +template +template +fvar& fvar::operator*=(fvar const& cr) { + using diff_t = typename std::array::difference_type; + promote const zero(0); + if BOOST_AUTODIFF_IF_CONSTEXPR (Order <= Order2) + for (size_t i = 0, j = Order; i <= Order; ++i, --j) + v[j] = std::inner_product(v.cbegin(), v.cend() - diff_t(i), cr.v.crbegin() + diff_t(i), zero); + else { + for (size_t i = 0, j = Order; i <= Order - Order2; ++i, --j) + v[j] = std::inner_product(cr.v.cbegin(), cr.v.cend(), v.crbegin() + diff_t(i), zero); + for (size_t i = Order - Order2 + 1, j = Order2 - 1; i <= Order; ++i, --j) + v[j] = std::inner_product(cr.v.cbegin(), cr.v.cbegin() + diff_t(j + 1), v.crbegin() + diff_t(i), zero); + } + return *this; +} + +template +fvar& fvar::operator*=(root_type const& ca) { + return multiply_assign_by_root_type(true, ca); +} + +template +template +fvar& fvar::operator/=(fvar const& cr) { + using diff_t = typename std::array::difference_type; + RealType const zero(0); + v.front() /= cr.v.front(); + if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2) + for (size_t i = 1, j = Order2 - 1, k = Order; i <= Order; ++i, --j, --k) + (v[i] -= std::inner_product( + cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), v.crbegin() + diff_t(k), zero)) /= cr.v.front(); + else if BOOST_AUTODIFF_IF_CONSTEXPR (0 < Order2) + for (size_t i = 1, j = Order2 - 1, k = Order; i <= Order; ++i, j && --j, --k) + (v[i] -= std::inner_product( + cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), v.crbegin() + diff_t(k), zero)) /= cr.v.front(); + else + for (size_t i = 1; i <= Order; ++i) + v[i] /= cr.v.front(); + return *this; +} + +template +fvar& fvar::operator/=(root_type const& ca) { + std::for_each(v.begin(), v.end(), [&ca](RealType& x) { x /= ca; }); + return *this; +} + +template +fvar fvar::operator-() const { + fvar retval(*this); + retval.negate(); + return retval; +} + +template +fvar const& fvar::operator+() const { + return *this; +} + +template +template +promote, fvar> fvar::operator+( + fvar const& cr) const { + promote, fvar> retval; + for (size_t i = 0; i <= (std::min)(Order, Order2); ++i) + retval.v[i] = v[i] + cr.v[i]; + if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2) + for (size_t i = Order + 1; i <= Order2; ++i) + retval.v[i] = cr.v[i]; + else if BOOST_AUTODIFF_IF_CONSTEXPR (Order2 < Order) + for (size_t i = Order2 + 1; i <= Order; ++i) + retval.v[i] = v[i]; + return retval; +} + +template +fvar fvar::operator+(root_type const& ca) const { + fvar retval(*this); + retval.v.front() += ca; + return retval; +} + +template +fvar operator+(typename fvar::root_type const& ca, + fvar const& cr) { + return cr + ca; +} + +template +template +promote, fvar> fvar::operator-( + fvar const& cr) const { + promote, fvar> retval; + for (size_t i = 0; i <= (std::min)(Order, Order2); ++i) + retval.v[i] = v[i] - cr.v[i]; + if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2) + for (auto i = Order + 1; i <= Order2; ++i) + retval.v[i] = -cr.v[i]; + else if BOOST_AUTODIFF_IF_CONSTEXPR (Order2 < Order) + for (auto i = Order2 + 1; i <= Order; ++i) + retval.v[i] = v[i]; + return retval; +} + +template +fvar fvar::operator-(root_type const& ca) const { + fvar retval(*this); + retval.v.front() -= ca; + return retval; +} + +template +fvar operator-(typename fvar::root_type const& ca, + fvar const& cr) { + fvar mcr = -cr; // Has same address as retval in operator-() due to NRVO. + mcr += ca; + return mcr; // <-- This allows for NRVO. The following does not. --> return mcr += ca; +} + +template +template +promote, fvar> fvar::operator*( + fvar const& cr) const { + using diff_t = typename std::array::difference_type; + promote const zero(0); + promote, fvar> retval; + if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2) + for (size_t i = 0, j = Order, k = Order2; i <= Order2; ++i, j && --j, --k) + retval.v[i] = std::inner_product(v.cbegin(), v.cend() - diff_t(j), cr.v.crbegin() + diff_t(k), zero); + else + for (size_t i = 0, j = Order2, k = Order; i <= Order; ++i, j && --j, --k) + retval.v[i] = std::inner_product(cr.v.cbegin(), cr.v.cend() - diff_t(j), v.crbegin() + diff_t(k), zero); + return retval; +} + +template +fvar fvar::operator*(root_type const& ca) const { + fvar retval(*this); + retval *= ca; + return retval; +} + +template +fvar operator*(typename fvar::root_type const& ca, + fvar const& cr) { + return cr * ca; +} + +template +template +promote, fvar> fvar::operator/( + fvar const& cr) const { + using diff_t = typename std::array::difference_type; + promote const zero(0); + promote, fvar> retval; + retval.v.front() = v.front() / cr.v.front(); + if BOOST_AUTODIFF_IF_CONSTEXPR (Order < Order2) { + for (size_t i = 1, j = Order2 - 1; i <= Order; ++i, --j) + retval.v[i] = + (v[i] - std::inner_product( + cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero)) / + cr.v.front(); + for (size_t i = Order + 1, j = Order2 - Order - 1; i <= Order2; ++i, --j) + retval.v[i] = + -std::inner_product( + cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero) / + cr.v.front(); + } else if BOOST_AUTODIFF_IF_CONSTEXPR (0 < Order2) + for (size_t i = 1, j = Order2 - 1, k = Order; i <= Order; ++i, j && --j, --k) + retval.v[i] = + (v[i] - std::inner_product( + cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(k), zero)) / + cr.v.front(); + else + for (size_t i = 1; i <= Order; ++i) + retval.v[i] = v[i] / cr.v.front(); + return retval; +} + +template +fvar fvar::operator/(root_type const& ca) const { + fvar retval(*this); + retval /= ca; + return retval; +} + +template +fvar operator/(typename fvar::root_type const& ca, + fvar const& cr) { + using diff_t = typename std::array::difference_type; + fvar retval; + retval.v.front() = ca / cr.v.front(); + if BOOST_AUTODIFF_IF_CONSTEXPR (0 < Order) { + RealType const zero(0); + for (size_t i = 1, j = Order - 1; i <= Order; ++i, --j) + retval.v[i] = + -std::inner_product( + cr.v.cbegin() + 1, cr.v.cend() - diff_t(j), retval.v.crbegin() + diff_t(j + 1), zero) / + cr.v.front(); + } + return retval; +} + +template +template +bool fvar::operator==(fvar const& cr) const { + return v.front() == cr.v.front(); +} + +template +bool fvar::operator==(root_type const& ca) const { + return v.front() == ca; +} + +template +bool operator==(typename fvar::root_type const& ca, fvar const& cr) { + return ca == cr.v.front(); +} + +template +template +bool fvar::operator!=(fvar const& cr) const { + return v.front() != cr.v.front(); +} + +template +bool fvar::operator!=(root_type const& ca) const { + return v.front() != ca; +} + +template +bool operator!=(typename fvar::root_type const& ca, fvar const& cr) { + return ca != cr.v.front(); +} + +template +template +bool fvar::operator<=(fvar const& cr) const { + return v.front() <= cr.v.front(); +} + +template +bool fvar::operator<=(root_type const& ca) const { + return v.front() <= ca; +} + +template +bool operator<=(typename fvar::root_type const& ca, fvar const& cr) { + return ca <= cr.v.front(); +} + +template +template +bool fvar::operator>=(fvar const& cr) const { + return v.front() >= cr.v.front(); +} + +template +bool fvar::operator>=(root_type const& ca) const { + return v.front() >= ca; +} + +template +bool operator>=(typename fvar::root_type const& ca, fvar const& cr) { + return ca >= cr.v.front(); +} + +template +template +bool fvar::operator<(fvar const& cr) const { + return v.front() < cr.v.front(); +} + +template +bool fvar::operator<(root_type const& ca) const { + return v.front() < ca; +} + +template +bool operator<(typename fvar::root_type const& ca, fvar const& cr) { + return ca < cr.v.front(); +} + +template +template +bool fvar::operator>(fvar const& cr) const { + return v.front() > cr.v.front(); +} + +template +bool fvar::operator>(root_type const& ca) const { + return v.front() > ca; +} + +template +bool operator>(typename fvar::root_type const& ca, fvar const& cr) { + return ca > cr.v.front(); +} + + /*** Other methods and functions ***/ + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +// f : order -> derivative(order)/factorial(order) +// Use this when you have the polynomial coefficients, rather than just the derivatives. E.g. See atan2(). +template +template +promote, Fvar, Fvars...> fvar::apply_coefficients( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + size_t i = (std::min)(order, order_sum); + promote, Fvar, Fvars...> accumulator = cr.apply_coefficients( + order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward(fvars)...); + while (i--) + (accumulator *= epsilon) += cr.apply_coefficients( + order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward(fvars)...); + return accumulator; +} +#endif + +// f : order -> derivative(order)/factorial(order) +// Use this when you have the polynomial coefficients, rather than just the derivatives. E.g. See atan(). +template +template +fvar fvar::apply_coefficients(size_t const order, Func const& f) const { + fvar const epsilon = fvar(*this).set_root(0); +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR + size_t i = (std::min)(order, order_sum); +#else // ODR-use of static constexpr + size_t i = order < order_sum ? order : order_sum; +#endif + fvar accumulator = f(i); + while (i--) + (accumulator *= epsilon) += f(i); + return accumulator; +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +// f : order -> derivative(order) +template +template +promote, Fvar, Fvars...> fvar::apply_coefficients_nonhorner( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + fvar epsilon_i = fvar(1); // epsilon to the power of i + promote, Fvar, Fvars...> accumulator = cr.apply_coefficients_nonhorner( + order, + [&f](auto... indices) { return f(0, static_cast(indices)...); }, + std::forward(fvars)...); + size_t const i_max = (std::min)(order, order_sum); + for (size_t i = 1; i <= i_max; ++i) { + epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); + accumulator += epsilon_i.epsilon_multiply( + i, + 0, + cr.apply_coefficients_nonhorner( + order - i, + [&f, i](auto... indices) { return f(i, static_cast(indices)...); }, + std::forward(fvars)...), + 0, + 0); + } + return accumulator; +} +#endif + +// f : order -> coefficient(order) +template +template +fvar fvar::apply_coefficients_nonhorner(size_t const order, + Func const& f) const { + fvar const epsilon = fvar(*this).set_root(0); + fvar epsilon_i = fvar(1); // epsilon to the power of i + fvar accumulator = fvar(f(0u)); +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR + size_t const i_max = (std::min)(order, order_sum); +#else // ODR-use of static constexpr + size_t const i_max = order < order_sum ? order : order_sum; +#endif + for (size_t i = 1; i <= i_max; ++i) { + epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); + accumulator += epsilon_i.epsilon_multiply(i, 0, f(i)); + } + return accumulator; +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +// f : order -> derivative(order) +template +template +promote, Fvar, Fvars...> fvar::apply_derivatives( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + size_t i = (std::min)(order, order_sum); + promote, Fvar, Fvars...> accumulator = + cr.apply_derivatives( + order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward(fvars)...) / + factorial(static_cast(i)); + while (i--) + (accumulator *= epsilon) += + cr.apply_derivatives( + order - i, [&f, i](auto... indices) { return f(i, indices...); }, std::forward(fvars)...) / + factorial(static_cast(i)); + return accumulator; +} +#endif + +// f : order -> derivative(order) +template +template +fvar fvar::apply_derivatives(size_t const order, Func const& f) const { + fvar const epsilon = fvar(*this).set_root(0); +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR + size_t i = (std::min)(order, order_sum); +#else // ODR-use of static constexpr + size_t i = order < order_sum ? order : order_sum; +#endif + fvar accumulator = f(i) / factorial(static_cast(i)); + while (i--) + (accumulator *= epsilon) += f(i) / factorial(static_cast(i)); + return accumulator; +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +// f : order -> derivative(order) +template +template +promote, Fvar, Fvars...> fvar::apply_derivatives_nonhorner( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + fvar epsilon_i = fvar(1); // epsilon to the power of i + promote, Fvar, Fvars...> accumulator = cr.apply_derivatives_nonhorner( + order, + [&f](auto... indices) { return f(0, static_cast(indices)...); }, + std::forward(fvars)...); + size_t const i_max = (std::min)(order, order_sum); + for (size_t i = 1; i <= i_max; ++i) { + epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); + accumulator += epsilon_i.epsilon_multiply( + i, + 0, + cr.apply_derivatives_nonhorner( + order - i, + [&f, i](auto... indices) { return f(i, static_cast(indices)...); }, + std::forward(fvars)...) / + factorial(static_cast(i)), + 0, + 0); + } + return accumulator; +} +#endif + +// f : order -> derivative(order) +template +template +fvar fvar::apply_derivatives_nonhorner(size_t const order, + Func const& f) const { + fvar const epsilon = fvar(*this).set_root(0); + fvar epsilon_i = fvar(1); // epsilon to the power of i + fvar accumulator = fvar(f(0u)); +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR + size_t const i_max = (std::min)(order, order_sum); +#else // ODR-use of static constexpr + size_t const i_max = order < order_sum ? order : order_sum; +#endif + for (size_t i = 1; i <= i_max; ++i) { + epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); + accumulator += epsilon_i.epsilon_multiply(i, 0, f(i) / factorial(static_cast(i))); + } + return accumulator; +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +// Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)" +template +template +get_type_at fvar::at(size_t order, Orders... orders) const { + if constexpr (0 < sizeof...(Orders)) + return v.at(order).at(static_cast(orders)...); + else + return v.at(order); +} +#endif + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +// Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)" +template +template +get_type_at, sizeof...(Orders)> fvar::derivative( + Orders... orders) const { + static_assert(sizeof...(Orders) <= depth, + "Number of parameters to derivative(...) cannot exceed fvar::depth."); + return at(static_cast(orders)...) * + (... * factorial(static_cast(orders))); +} +#endif + +template +const RealType& fvar::operator[](size_t i) const { + return v[i]; +} + +template +RealType fvar::epsilon_inner_product(size_t z0, + size_t const isum0, + size_t const m0, + fvar const& cr, + size_t z1, + size_t const isum1, + size_t const m1, + size_t const j) const { + static_assert(is_fvar::value, "epsilon_inner_product() must have 1 < depth."); + RealType accumulator = RealType(); + auto const i0_max = m1 < j ? j - m1 : 0; + for (auto i0 = m0, i1 = j - m0; i0 <= i0_max; ++i0, --i1) + accumulator += v[i0].epsilon_multiply(z0, isum0 + i0, cr.v[i1], z1, isum1 + i1); + return accumulator; +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +template +fvar fvar::epsilon_multiply(size_t z0, + size_t isum0, + fvar const& cr, + size_t z1, + size_t isum1) const { + using diff_t = typename std::array::difference_type; + RealType const zero(0); + size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; + size_t const m1 = order_sum + isum1 < Order + z1 ? Order + z1 - (order_sum + isum1) : 0; + size_t const i_max = m0 + m1 < Order ? Order - (m0 + m1) : 0; + fvar retval = fvar(); + if constexpr (is_fvar::value) + for (size_t i = 0, j = Order; i <= i_max; ++i, --j) + retval.v[j] = epsilon_inner_product(z0, isum0, m0, cr, z1, isum1, m1, j); + else + for (size_t i = 0, j = Order; i <= i_max; ++i, --j) + retval.v[j] = std::inner_product( + v.cbegin() + diff_t(m0), v.cend() - diff_t(i + m1), cr.v.crbegin() + diff_t(i + m0), zero); + return retval; +} +#endif + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +// When called from outside this method, z0 should be non-zero. Otherwise if z0=0 then it will give an +// incorrect result of 0 when the root value is 0 and ca=inf, when instead the correct product is nan. +// If z0=0 then use the regular multiply operator*() instead. +template +fvar fvar::epsilon_multiply(size_t z0, + size_t isum0, + root_type const& ca) const { + fvar retval(*this); + size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; + if constexpr (is_fvar::value) + for (size_t i = m0; i <= Order; ++i) + retval.v[i] = retval.v[i].epsilon_multiply(z0, isum0 + i, ca); + else + for (size_t i = m0; i <= Order; ++i) + if (retval.v[i] != static_cast(0)) + retval.v[i] *= ca; + return retval; +} +#endif + +template +fvar fvar::inverse() const { + return static_cast(*this) == 0 ? inverse_apply() : 1 / *this; +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +template +fvar& fvar::negate() { + if constexpr (is_fvar::value) + std::for_each(v.begin(), v.end(), [](RealType& r) { r.negate(); }); + else + std::for_each(v.begin(), v.end(), [](RealType& a) { a = -a; }); + return *this; +} +#endif + +// This gives log(0.0) = depth(1)(-inf,inf,-inf,inf,-inf,inf) +// 1 / *this: log(0.0) = depth(1)(-inf,inf,-inf,-nan,-nan,-nan) +template +fvar fvar::inverse_apply() const { + root_type derivatives[order_sum + 1]; // LCOV_EXCL_LINE This causes a false negative on lcov coverage test. + root_type const x0 = static_cast(*this); + *derivatives = 1 / x0; + for (size_t i = 1; i <= order_sum; ++i) + derivatives[i] = -derivatives[i - 1] * i / x0; + return apply_derivatives_nonhorner(order_sum, [&derivatives](size_t j) { return derivatives[j]; }); +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +template +fvar& fvar::multiply_assign_by_root_type(bool is_root, + root_type const& ca) { + auto itr = v.begin(); + if constexpr (is_fvar::value) { + itr->multiply_assign_by_root_type(is_root, ca); + for (++itr; itr != v.end(); ++itr) + itr->multiply_assign_by_root_type(false, ca); + } else { + if (is_root || *itr != 0) + *itr *= ca; // Skip multiplication of 0 by ca=inf to avoid nan, except when is_root. + for (++itr; itr != v.end(); ++itr) + if (*itr != 0) + *itr *= ca; + } + return *this; +} +#endif + +template +fvar::operator root_type() const { + return static_cast(v.front()); +} + +template +template +fvar::operator T() const { + return static_cast(static_cast(v.front())); +} + +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR +template +fvar& fvar::set_root(root_type const& root) { + if constexpr (is_fvar::value) + v.front().set_root(root); + else + v.front() = root; + return *this; +} +#endif + +// Standard Library Support Requirements + +template +fvar fabs(fvar const& cr) { + typename fvar::root_type const zero(0); + return cr < zero ? -cr + : cr == zero ? fvar() // Canonical fabs'(0) = 0. + : cr; // Propagate NaN. +} + +template +fvar abs(fvar const& cr) { + return fabs(cr); +} + +template +fvar ceil(fvar const& cr) { + using std::ceil; + return fvar(ceil(static_cast::root_type>(cr))); +} + +template +fvar floor(fvar const& cr) { + using std::floor; + return fvar(floor(static_cast::root_type>(cr))); +} + +template +fvar exp(fvar const& cr) { + using std::exp; + constexpr size_t order = fvar::order_sum; + using root_type = typename fvar::root_type; + root_type const d0 = exp(static_cast(cr)); + return cr.apply_derivatives(order, [&d0](size_t) { return d0; }); +} + +template +fvar pow(fvar const& x, + typename fvar::root_type const& y) { + using std::pow; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const x0 = static_cast(x); + root_type derivatives[order + 1]{pow(x0, y)}; + for (size_t i = 0; i < order && y - i != 0; ++i) + derivatives[i + 1] = (y - i) * derivatives[i] / x0; + return x.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i]; }); +} + +template +fvar pow(typename fvar::root_type const& x, + fvar const& y) { + BOOST_MATH_STD_USING + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const y0 = static_cast(y); + root_type derivatives[order + 1]; + *derivatives = pow(x, y0); + root_type const logx = log(x); + for (size_t i = 0; i < order; ++i) + derivatives[i + 1] = derivatives[i] * logx; + return y.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i]; }); +} + +template +promote, fvar> pow(fvar const& x, + fvar const& y) { + BOOST_MATH_STD_USING + using return_type = promote, fvar>; + using root_type = typename return_type::root_type; + constexpr size_t order = return_type::order_sum; + root_type const x0 = static_cast(x); + root_type const y0 = static_cast(y); + root_type dxydx[order + 1]{pow(x0, y0)}; + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return return_type(*dxydx); + else { + for (size_t i = 0; i < order && y0 - i != 0; ++i) + dxydx[i + 1] = (y0 - i) * dxydx[i] / x0; + std::array, order + 1> lognx; + lognx.front() = fvar(1); +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR + lognx[1] = log(make_fvar(x0)); +#else // for compilers that compile this branch when order=0. + lognx[(std::min)(size_t(1), order)] = log(make_fvar(x0)); +#endif + for (size_t i = 1; i < order; ++i) + lognx[i + 1] = lognx[i] * lognx[1]; + auto const f = [&dxydx, &lognx](size_t i, size_t j) { + size_t binomial = 1; + root_type sum = dxydx[i] * static_cast(lognx[j]); + for (size_t k = 1; k <= i; ++k) { + (binomial *= (i - k + 1)) /= k; // binomial_coefficient(i,k) + sum += binomial * dxydx[i - k] * lognx[j].derivative(k); + } + return sum; + }; + if (fabs(x0) < std::numeric_limits::epsilon()) + return x.apply_derivatives_nonhorner(order, f, y); + return x.apply_derivatives(order, f, y); + } +} + +template +fvar sqrt(fvar const& cr) { + using std::sqrt; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type derivatives[order + 1]; + root_type const x = static_cast(cr); + *derivatives = sqrt(x); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(*derivatives); + else { + root_type numerator = 0.5; + root_type powers = 1; +#ifndef BOOST_NO_CXX17_IF_CONSTEXPR + derivatives[1] = numerator / *derivatives; +#else // for compilers that compile this branch when order=0. + derivatives[(std::min)(size_t(1), order)] = numerator / *derivatives; +#endif + using diff_t = typename std::array::difference_type; + for (size_t i = 2; i <= order; ++i) { + numerator *= static_cast(-0.5) * ((static_cast(i) << 1) - 3); + powers *= x; + derivatives[i] = numerator / (powers * *derivatives); + } + auto const f = [&derivatives](size_t i) { return derivatives[i]; }; + if (cr < std::numeric_limits::epsilon()) + return cr.apply_derivatives_nonhorner(order, f); + return cr.apply_derivatives(order, f); + } +} + +// Natural logarithm. If cr==0 then derivative(i) may have nans due to nans from inverse(). +template +fvar log(fvar const& cr) { + using std::log; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = log(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto const d1 = make_fvar(static_cast(cr)).inverse(); // log'(x) = 1 / x + return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar frexp(fvar const& cr, int* exp) { + using multiprecision::exp2; + using std::exp2; + using std::frexp; + using root_type = typename fvar::root_type; + frexp(static_cast(cr), exp); + return cr * static_cast(exp2(-*exp)); +} + +template +fvar ldexp(fvar const& cr, int exp) { + // argument to std::exp2 must be casted to root_type, otherwise std::exp2 returns double (always) + using multiprecision::exp2; + using std::exp2; + return cr * exp2(static_cast::root_type>(exp)); +} + +template +fvar cos(fvar const& cr) { + BOOST_MATH_STD_USING + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = cos(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + root_type const d1 = -sin(static_cast(cr)); + root_type const derivatives[4]{d0, d1, -d0, -d1}; + return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 3]; }); + } +} + +template +fvar sin(fvar const& cr) { + BOOST_MATH_STD_USING + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = sin(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + root_type const d1 = cos(static_cast(cr)); + root_type const derivatives[4]{d0, d1, -d0, -d1}; + return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 3]; }); + } +} + +template +fvar asin(fvar const& cr) { + using std::asin; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = asin(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); + auto const d1 = sqrt((x *= x).negate() += 1).inverse(); // asin'(x) = 1 / sqrt(1-x*x). + return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar tan(fvar const& cr) { + using std::tan; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = tan(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto c = cos(make_fvar(static_cast(cr))); + auto const d1 = (c *= c).inverse(); // tan'(x) = 1 / cos(x)^2 + return cr.apply_coefficients_nonhorner(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar atan(fvar const& cr) { + using std::atan; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = atan(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); + auto const d1 = ((x *= x) += 1).inverse(); // atan'(x) = 1 / (x*x+1). + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar atan2(fvar const& cr, + typename fvar::root_type const& ca) { + using std::atan2; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = atan2(static_cast(cr), ca); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto y = make_fvar(static_cast(cr)); + auto const d1 = ca / ((y *= y) += (ca * ca)); // (d/dy)atan2(y,x) = x / (y*y+x*x) + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar atan2(typename fvar::root_type const& ca, + fvar const& cr) { + using std::atan2; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = atan2(ca, static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); + auto const d1 = -ca / ((x *= x) += (ca * ca)); // (d/dx)atan2(y,x) = -y / (x*x+y*y) + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +promote, fvar> atan2(fvar const& cr1, + fvar const& cr2) { + using std::atan2; + using return_type = promote, fvar>; + using root_type = typename return_type::root_type; + constexpr size_t order = return_type::order_sum; + root_type const y = static_cast(cr1); + root_type const x = static_cast(cr2); + root_type const d00 = atan2(y, x); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return return_type(d00); + else { + constexpr size_t order1 = fvar::order_sum; + constexpr size_t order2 = fvar::order_sum; + auto x01 = make_fvar::root_type, order2 - 1>(x); + auto const d01 = -y / ((x01 *= x01) += (y * y)); + auto y10 = make_fvar::root_type, order1 - 1>(y); + auto x10 = make_fvar::root_type, 0, order2>(x); + auto const d10 = x10 / ((x10 * x10) + (y10 *= y10)); + auto const f = [&d00, &d01, &d10](size_t i, size_t j) { + return i ? d10[i - 1][j] / i : j ? d01[j - 1] / j : d00; + }; + return cr1.apply_coefficients(order, f, cr2); + } +} + +template +promote, fvar> fmod(fvar const& cr1, + fvar const& cr2) { + using boost::math::trunc; + auto const numer = static_cast::root_type>(cr1); + auto const denom = static_cast::root_type>(cr2); + return cr1 - cr2 * trunc(numer / denom); +} + +template +fvar round(fvar const& cr) { + using boost::math::round; + return fvar(round(static_cast::root_type>(cr))); +} + +template +int iround(fvar const& cr) { + using boost::math::iround; + return iround(static_cast::root_type>(cr)); +} + +template +long lround(fvar const& cr) { + using boost::math::lround; + return lround(static_cast::root_type>(cr)); +} + +template +long long llround(fvar const& cr) { + using boost::math::llround; + return llround(static_cast::root_type>(cr)); +} + +template +fvar trunc(fvar const& cr) { + using boost::math::trunc; + return fvar(trunc(static_cast::root_type>(cr))); +} + +template +long double truncl(fvar const& cr) { + using std::truncl; + return truncl(static_cast::root_type>(cr)); +} + +template +int itrunc(fvar const& cr) { + using boost::math::itrunc; + return itrunc(static_cast::root_type>(cr)); +} + +template +long long lltrunc(fvar const& cr) { + using boost::math::lltrunc; + return lltrunc(static_cast::root_type>(cr)); +} + +template +std::ostream& operator<<(std::ostream& out, fvar const& cr) { + out << "depth(" << cr.depth << ")(" << cr.v.front(); + for (size_t i = 1; i <= Order; ++i) + out << ',' << cr.v[i]; + return out << ')'; +} + +// Additional functions + +template +fvar acos(fvar const& cr) { + using std::acos; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = acos(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); + auto const d1 = sqrt((x *= x).negate() += 1).inverse().negate(); // acos'(x) = -1 / sqrt(1-x*x). + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar acosh(fvar const& cr) { + using boost::math::acosh; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = acosh(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); + auto const d1 = sqrt((x *= x) -= 1).inverse(); // acosh'(x) = 1 / sqrt(x*x-1). + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar asinh(fvar const& cr) { + using boost::math::asinh; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = asinh(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); + auto const d1 = sqrt((x *= x) += 1).inverse(); // asinh'(x) = 1 / sqrt(x*x+1). + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar atanh(fvar const& cr) { + using boost::math::atanh; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = atanh(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); + auto const d1 = ((x *= x).negate() += 1).inverse(); // atanh'(x) = 1 / (1-x*x) + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar cosh(fvar const& cr) { + BOOST_MATH_STD_USING + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = cosh(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + root_type const derivatives[2]{d0, sinh(static_cast(cr))}; + return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 1]; }); + } +} + +template +fvar digamma(fvar const& cr) { + using boost::math::digamma; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const x = static_cast(cr); + root_type const d0 = digamma(x); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + static_assert(order <= static_cast(std::numeric_limits::max()), + "order exceeds maximum derivative for boost::math::polygamma()."); + return cr.apply_derivatives( + order, [&x, &d0](size_t i) { return i ? boost::math::polygamma(static_cast(i), x) : d0; }); + } +} + +template +fvar erf(fvar const& cr) { + using boost::math::erf; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = erf(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); // d1 = 2/sqrt(pi)*exp(-x*x) + auto const d1 = 2 * constants::one_div_root_pi() * exp((x *= x).negate()); + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar erfc(fvar const& cr) { + using boost::math::erfc; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = erfc(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + auto x = make_fvar(static_cast(cr)); // erfc'(x) = -erf'(x) + auto const d1 = -2 * constants::one_div_root_pi() * exp((x *= x).negate()); + return cr.apply_coefficients(order, [&d0, &d1](size_t i) { return i ? d1[i - 1] / i : d0; }); + } +} + +template +fvar lambert_w0(fvar const& cr) { + using std::exp; + using boost::math::lambert_w0; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type derivatives[order + 1]; + *derivatives = lambert_w0(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(*derivatives); + else { + root_type const expw = exp(*derivatives); + derivatives[1] = 1 / (static_cast(cr) + expw); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 1) + return cr.apply_derivatives_nonhorner(order, [&derivatives](size_t i) { return derivatives[i]; }); + else { + using diff_t = typename std::array::difference_type; + root_type d1powers = derivatives[1] * derivatives[1]; + root_type const x = derivatives[1] * expw; + derivatives[2] = d1powers * (-1 - x); + std::array coef{{-1, -1}}; // as in derivatives[2]. + for (size_t n = 3; n <= order; ++n) { + coef[n - 1] = coef[n - 2] * -static_cast(2 * n - 3); + for (size_t j = n - 2; j != 0; --j) + (coef[j] *= -static_cast(n - 1)) -= (n + j - 2) * coef[j - 1]; + coef[0] *= -static_cast(n - 1); + d1powers *= derivatives[1]; + derivatives[n] = + d1powers * std::accumulate(coef.crend() - diff_t(n - 1), + coef.crend(), + coef[n - 1], + [&x](root_type const& a, root_type const& b) { return a * x + b; }); + } + return cr.apply_derivatives_nonhorner(order, [&derivatives](size_t i) { return derivatives[i]; }); + } + } +} + +template +fvar lgamma(fvar const& cr) { + using std::lgamma; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const x = static_cast(cr); + root_type const d0 = lgamma(x); + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(d0); + else { + static_assert(order <= static_cast(std::numeric_limits::max()) + 1, + "order exceeds maximum derivative for boost::math::polygamma()."); + return cr.apply_derivatives( + order, [&x, &d0](size_t i) { return i ? boost::math::polygamma(static_cast(i - 1), x) : d0; }); + } +} + +template +fvar sinc(fvar const& cr) { + if (cr != 0) + return sin(cr) / cr; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type taylor[order + 1]{1}; // sinc(0) = 1 + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(*taylor); + else { + for (size_t n = 2; n <= order; n += 2) + taylor[n] = (1 - static_cast(n & 2)) / factorial(static_cast(n + 1)); + return cr.apply_coefficients_nonhorner(order, [&taylor](size_t i) { return taylor[i]; }); + } +} + +template +fvar sinh(fvar const& cr) { + BOOST_MATH_STD_USING + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + root_type const d0 = sinh(static_cast(cr)); + if BOOST_AUTODIFF_IF_CONSTEXPR (fvar::order_sum == 0) + return fvar(d0); + else { + root_type const derivatives[2]{d0, cosh(static_cast(cr))}; + return cr.apply_derivatives(order, [&derivatives](size_t i) { return derivatives[i & 1]; }); + } +} + +template +fvar tanh(fvar const& cr) { + fvar retval = exp(cr * 2); + fvar const denom = retval + 1; + (retval -= 1) /= denom; + return retval; +} + +template +fvar tgamma(fvar const& cr) { + using std::tgamma; + using root_type = typename fvar::root_type; + constexpr size_t order = fvar::order_sum; + if BOOST_AUTODIFF_IF_CONSTEXPR (order == 0) + return fvar(tgamma(static_cast(cr))); + else { + if (cr < 0) + return constants::pi() / (sin(constants::pi() * cr) * tgamma(1 - cr)); + return exp(lgamma(cr)).set_root(tgamma(static_cast(cr))); + } +} + +} // namespace detail +} // namespace autodiff_v1 +} // namespace differentiation +} // namespace math +} // namespace boost + +namespace std { + +// boost::math::tools::digits() is handled by this std::numeric_limits<> specialization, +// and similarly for max_value, min_value, log_max_value, log_min_value, and epsilon. +template +class numeric_limits> + : public numeric_limits::root_type> { +}; + +} // namespace std + +namespace boost { +namespace math { +namespace tools { +namespace detail { + +template +using autodiff_fvar_type = differentiation::detail::fvar; + +template +using autodiff_root_type = typename autodiff_fvar_type::root_type; +} // namespace detail + +// See boost/math/tools/promotion.hpp +template +struct promote_args_2, + detail::autodiff_fvar_type> { + using type = detail::autodiff_fvar_type::type, +#ifndef BOOST_NO_CXX14_CONSTEXPR + (std::max)(Order0, Order1)>; +#else + Order0; +#endif +}; + +template +struct promote_args> { + using type = detail::autodiff_fvar_type::type, Order>; +}; + +template +struct promote_args_2, RealType1> { + using type = detail::autodiff_fvar_type::type, Order0>; +}; + +template +struct promote_args_2> { + using type = detail::autodiff_fvar_type::type, Order1>; +}; + +template +inline BOOST_MATH_CONSTEXPR destination_t real_cast(detail::autodiff_fvar_type const& from_v) + BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(destination_t) && BOOST_MATH_IS_FLOAT(RealType)) { + return real_cast(static_cast>(from_v)); +} + +} // namespace tools + +namespace policies { + +template +using fvar_t = differentiation::detail::fvar; +template +struct evaluation, Policy> { + using type = fvar_t::type, Order>; +}; + +template +struct evaluation, Policy> { + using type = + fvar_t::type, Order>; +}; + +} // namespace policies +} // namespace math +} // namespace boost + +#ifdef BOOST_NO_CXX17_IF_CONSTEXPR +#include "autodiff_cpp11.hpp" +#endif + +#endif // BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP diff --git a/boost/math/differentiation/autodiff_cpp11.hpp b/boost/math/differentiation/autodiff_cpp11.hpp new file mode 100644 index 0000000000..42140a4824 --- /dev/null +++ b/boost/math/differentiation/autodiff_cpp11.hpp @@ -0,0 +1,383 @@ +// Copyright Matthew Pulver 2018 - 2019. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +// Contributors: +// * Kedar R. Bhat - C++11 compatibility. + +// Notes: +// * Any changes to this file should always be downstream from autodiff.cpp. +// C++17 is a higher-level language and is easier to maintain. For example, a number of functions which are +// lucidly read in autodiff.cpp are forced to be split into multiple structs/functions in this file for +// C++11. +// * Use of typename RootType and SizeType is a hack to prevent Visual Studio 2015 from compiling functions +// that are never called, that would otherwise produce compiler errors. Also forces functions to be inline. + +#ifndef BOOST_MATH_DIFFERENTIATION_AUTODIFF_HPP +#error \ + "Do not #include this file directly. This should only be #included by autodiff.hpp for C++11 compatibility." +#endif + +#include + +namespace boost { +namespace math { +namespace differentiation { +inline namespace autodiff_v1 { +namespace detail { + +template +fvar::fvar(root_type const& ca, bool const is_variable) { + fvar_cpp11(is_fvar{}, ca, is_variable); +} + +template +template +void fvar::fvar_cpp11(std::true_type, RootType const& ca, bool const is_variable) { + v.front() = RealType(ca, is_variable); + if (0 < Order) + std::fill(v.begin() + 1, v.end(), static_cast(0)); +} + +template +template +void fvar::fvar_cpp11(std::false_type, RootType const& ca, bool const is_variable) { + v.front() = ca; + if (0 < Order) { + v[1] = static_cast(static_cast(is_variable)); + if (1 < Order) + std::fill(v.begin() + 2, v.end(), static_cast(0)); + } +} + +template +template +get_type_at fvar::at_cpp11(std::true_type, + size_t order, + Orders...) const { + return v.at(order); +} + +template +template +get_type_at fvar::at_cpp11(std::false_type, + size_t order, + Orders... orders) const { + return v.at(order).at(orders...); +} + +// Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)" +template +template +get_type_at fvar::at(size_t order, Orders... orders) const { + return at_cpp11(std::integral_constant{}, order, orders...); +} + +template +constexpr T product(Ts...) { + return static_cast(1); +} + +template +constexpr T product(T factor, Ts... factors) { + return factor * product(factors...); +} + +// Can throw "std::out_of_range: array::at: __n (which is 7) >= _Nm (which is 7)" +template +template +get_type_at, sizeof...(Orders)> fvar::derivative( + Orders... orders) const { + static_assert(sizeof...(Orders) <= depth, + "Number of parameters to derivative(...) cannot exceed fvar::depth."); + return at(static_cast(orders)...) * + product(boost::math::factorial(static_cast(orders))...); +} + +template +class Curry { + Func const& f_; + size_t const i_; + + public: + template // typename SizeType to force inline constructor. + Curry(Func const& f, SizeType i) : f_(f), i_(static_cast(i)) {} + template + RootType operator()(Indices... indices) const { + using unsigned_t = typename std::make_unsigned::type...>::type; + return f_(i_, static_cast(indices)...); + } +}; + +template +template +promote, Fvar, Fvars...> fvar::apply_coefficients( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + size_t i = order < order_sum ? order : order_sum; + using return_type = promote, Fvar, Fvars...>; + return_type accumulator = cr.apply_coefficients( + order - i, Curry(f, i), std::forward(fvars)...); + while (i--) + (accumulator *= epsilon) += cr.apply_coefficients( + order - i, Curry(f, i), std::forward(fvars)...); + return accumulator; +} + +template +template +promote, Fvar, Fvars...> fvar::apply_coefficients_nonhorner( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + fvar epsilon_i = fvar(1); // epsilon to the power of i + using return_type = promote, Fvar, Fvars...>; + return_type accumulator = cr.apply_coefficients_nonhorner( + order, Curry(f, 0), std::forward(fvars)...); + size_t const i_max = order < order_sum ? order : order_sum; + for (size_t i = 1; i <= i_max; ++i) { + epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); + accumulator += epsilon_i.epsilon_multiply( + i, + 0, + cr.apply_coefficients_nonhorner( + order - i, Curry(f, i), std::forward(fvars)...), + 0, + 0); + } + return accumulator; +} + +template +template +promote, Fvar, Fvars...> fvar::apply_derivatives( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + size_t i = order < order_sum ? order : order_sum; + using return_type = promote, Fvar, Fvars...>; + return_type accumulator = + cr.apply_derivatives( + order - i, Curry(f, i), std::forward(fvars)...) / + factorial(static_cast(i)); + while (i--) + (accumulator *= epsilon) += + cr.apply_derivatives( + order - i, Curry(f, i), std::forward(fvars)...) / + factorial(static_cast(i)); + return accumulator; +} + +template +template +promote, Fvar, Fvars...> fvar::apply_derivatives_nonhorner( + size_t const order, + Func const& f, + Fvar const& cr, + Fvars&&... fvars) const { + fvar const epsilon = fvar(*this).set_root(0); + fvar epsilon_i = fvar(1); // epsilon to the power of i + using return_type = promote, Fvar, Fvars...>; + return_type accumulator = cr.apply_derivatives_nonhorner( + order, Curry(f, 0), std::forward(fvars)...); + size_t const i_max = order < order_sum ? order : order_sum; + for (size_t i = 1; i <= i_max; ++i) { + epsilon_i = epsilon_i.epsilon_multiply(i - 1, 0, epsilon, 1, 0); + accumulator += epsilon_i.epsilon_multiply( + i, + 0, + cr.apply_derivatives_nonhorner( + order - i, Curry(f, i), std::forward(fvars)...) / + factorial(static_cast(i)), + 0, + 0); + } + return accumulator; +} + +template +template +fvar fvar::epsilon_multiply_cpp11(std::true_type, + SizeType z0, + size_t isum0, + fvar const& cr, + size_t z1, + size_t isum1) const { + size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; + size_t const m1 = order_sum + isum1 < Order + z1 ? Order + z1 - (order_sum + isum1) : 0; + size_t const i_max = m0 + m1 < Order ? Order - (m0 + m1) : 0; + fvar retval = fvar(); + for (size_t i = 0, j = Order; i <= i_max; ++i, --j) + retval.v[j] = epsilon_inner_product(z0, isum0, m0, cr, z1, isum1, m1, j); + return retval; +} + +template +template +fvar fvar::epsilon_multiply_cpp11(std::false_type, + SizeType z0, + size_t isum0, + fvar const& cr, + size_t z1, + size_t isum1) const { + using ssize_t = typename std::make_signed::type; + RealType const zero(0); + size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; + size_t const m1 = order_sum + isum1 < Order + z1 ? Order + z1 - (order_sum + isum1) : 0; + size_t const i_max = m0 + m1 < Order ? Order - (m0 + m1) : 0; + fvar retval = fvar(); + for (size_t i = 0, j = Order; i <= i_max; ++i, --j) + retval.v[j] = std::inner_product( + v.cbegin() + ssize_t(m0), v.cend() - ssize_t(i + m1), cr.v.crbegin() + ssize_t(i + m0), zero); + return retval; +} + +template +fvar fvar::epsilon_multiply(size_t z0, + size_t isum0, + fvar const& cr, + size_t z1, + size_t isum1) const { + return epsilon_multiply_cpp11(is_fvar{}, z0, isum0, cr, z1, isum1); +} + +template +template +fvar fvar::epsilon_multiply_cpp11(std::true_type, + SizeType z0, + size_t isum0, + root_type const& ca) const { + fvar retval(*this); + size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; + for (size_t i = m0; i <= Order; ++i) + retval.v[i] = retval.v[i].epsilon_multiply(z0, isum0 + i, ca); + return retval; +} + +template +template +fvar fvar::epsilon_multiply_cpp11(std::false_type, + SizeType z0, + size_t isum0, + root_type const& ca) const { + fvar retval(*this); + size_t const m0 = order_sum + isum0 < Order + z0 ? Order + z0 - (order_sum + isum0) : 0; + for (size_t i = m0; i <= Order; ++i) + if (retval.v[i] != static_cast(0)) + retval.v[i] *= ca; + return retval; +} + +template +fvar fvar::epsilon_multiply(size_t z0, + size_t isum0, + root_type const& ca) const { + return epsilon_multiply_cpp11(is_fvar{}, z0, isum0, ca); +} + +template +template +fvar& fvar::multiply_assign_by_root_type_cpp11(std::true_type, + bool is_root, + RootType const& ca) { + auto itr = v.begin(); + itr->multiply_assign_by_root_type(is_root, ca); + for (++itr; itr != v.end(); ++itr) + itr->multiply_assign_by_root_type(false, ca); + return *this; +} + +template +template +fvar& fvar::multiply_assign_by_root_type_cpp11(std::false_type, + bool is_root, + RootType const& ca) { + auto itr = v.begin(); + if (is_root || *itr != 0) + *itr *= ca; // Skip multiplication of 0 by ca=inf to avoid nan, except when is_root. + for (++itr; itr != v.end(); ++itr) + if (*itr != 0) + *itr *= ca; + return *this; +} + +template +fvar& fvar::multiply_assign_by_root_type(bool is_root, + root_type const& ca) { + return multiply_assign_by_root_type_cpp11(is_fvar{}, is_root, ca); +} + +template +template +fvar& fvar::negate_cpp11(std::true_type, RootType const&) { + std::for_each(v.begin(), v.end(), [](RealType& r) { r.negate(); }); + return *this; +} + +template +template +fvar& fvar::negate_cpp11(std::false_type, RootType const&) { + std::for_each(v.begin(), v.end(), [](RealType& a) { a = -a; }); + return *this; +} + +template +fvar& fvar::negate() { + return negate_cpp11(is_fvar{}, static_cast(*this)); +} + +template +template +fvar& fvar::set_root_cpp11(std::true_type, RootType const& root) { + v.front().set_root(root); + return *this; +} + +template +template +fvar& fvar::set_root_cpp11(std::false_type, RootType const& root) { + v.front() = root; + return *this; +} + +template +fvar& fvar::set_root(root_type const& root) { + return set_root_cpp11(is_fvar{}, root); +} + +template +auto make_fvar_for_tuple(mp11::index_sequence, RealType const& ca) + -> decltype(make_fvar::value..., Order>(ca)) { + return make_fvar::value..., Order>(ca); +} + +template +auto make_ftuple_impl(mp11::index_sequence, RealTypes const&... ca) + -> decltype(std::make_tuple(make_fvar_for_tuple(mp11::make_index_sequence{}, + ca)...)) { + return std::make_tuple(make_fvar_for_tuple(mp11::make_index_sequence{}, ca)...); +} + +} // namespace detail + +template +auto make_ftuple(RealTypes const&... ca) + -> decltype(detail::make_ftuple_impl(mp11::index_sequence_for{}, + ca...)) { + static_assert(sizeof...(Orders) == sizeof...(RealTypes), + "Number of Orders must match number of function parameters."); + return detail::make_ftuple_impl(mp11::index_sequence_for{}, ca...); +} + +} // namespace autodiff_v1 +} // namespace differentiation +} // namespace math +} // namespace boost diff --git a/boost/math/interpolators/cardinal_quadratic_b_spline.hpp b/boost/math/interpolators/cardinal_quadratic_b_spline.hpp new file mode 100644 index 0000000000..a5b150f2f1 --- /dev/null +++ b/boost/math/interpolators/cardinal_quadratic_b_spline.hpp @@ -0,0 +1,57 @@ +// Copyright Nick Thompson, 2019 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_HPP +#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_HPP +#include +#include + + +namespace boost{ namespace math{ namespace interpolators { + +template +class cardinal_quadratic_b_spline +{ +public: + // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them. + // y[0] = y(a), y[n - 1] = y(b), step_size = (b - a)/(n -1). + cardinal_quadratic_b_spline(const Real* const y, + size_t n, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), + Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()) + : impl_(std::make_shared>(y, n, t0, h, left_endpoint_derivative, right_endpoint_derivative)) + {} + + // Oh the bizarre error messages if we template this on a RandomAccessContainer: + cardinal_quadratic_b_spline(std::vector const & y, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), + Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()) + : impl_(std::make_shared>(y.data(), y.size(), t0, h, left_endpoint_derivative, right_endpoint_derivative)) + {} + + + Real operator()(Real t) const { + return impl_->operator()(t); + } + + Real prime(Real t) const { + return impl_->prime(t); + } + + Real t_max() const { + return impl_->t_max(); + } + +private: + std::shared_ptr> impl_; +}; + +}}} +#endif diff --git a/boost/math/interpolators/catmull_rom.hpp b/boost/math/interpolators/catmull_rom.hpp index 05ad666de2..20207871db 100644 --- a/boost/math/interpolators/catmull_rom.hpp +++ b/boost/math/interpolators/catmull_rom.hpp @@ -13,6 +13,25 @@ #include #include #include +#include + +namespace std_workaround { + +#if defined(__cpp_lib_nonmember_container_access) || (defined(BOOST_MSVC) && (BOOST_MSVC >= 1900)) + using std::size; +#else + template + inline BOOST_CONSTEXPR std::size_t size(const C& c) + { + return c.size(); + } + template + inline BOOST_CONSTEXPR std::size_t size(const T(&array)[N]) BOOST_NOEXCEPT + { + return N; + } +#endif +} namespace boost{ namespace math{ @@ -22,7 +41,7 @@ namespace boost{ namespace math{ typename Point::value_type alpha_distance(Point const & p1, Point const & p2, typename Point::value_type alpha) { using std::pow; - using std::size; + using std_workaround::size; typename Point::value_type dsq = 0; for (size_t i = 0; i < size(p1); ++i) { @@ -33,44 +52,45 @@ namespace boost{ namespace math{ } } -template +template > class catmull_rom { + typedef typename Point::value_type value_type; public: - catmull_rom(std::vector&& points, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2); + catmull_rom(RandomAccessContainer&& points, bool closed = false, value_type alpha = (value_type) 1/ (value_type) 2); - catmull_rom(std::initializer_list l, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2) : catmull_rom(std::vector(l), closed, alpha) {} + catmull_rom(std::initializer_list l, bool closed = false, value_type alpha = (value_type) 1/ (value_type) 2) : catmull_rom(RandomAccessContainer(l), closed, alpha) {} - typename Point::value_type max_parameter() const + value_type max_parameter() const { return m_max_s; } - typename Point::value_type parameter_at_point(size_t i) const + value_type parameter_at_point(size_t i) const { return m_s[i+1]; } - Point operator()(const typename Point::value_type s) const; + Point operator()(const value_type s) const; - Point prime(const typename Point::value_type s) const; + Point prime(const value_type s) const; - std::vector&& get_points() + RandomAccessContainer&& get_points() { return std::move(m_pnts); } private: - std::vector m_pnts; - std::vector m_s; - typename Point::value_type m_max_s; + RandomAccessContainer m_pnts; + std::vector m_s; + value_type m_max_s; }; -template -catmull_rom::catmull_rom(std::vector&& points, bool closed, typename Point::value_type alpha) : m_pnts(std::move(points)) +template +catmull_rom::catmull_rom(RandomAccessContainer&& points, bool closed, typename Point::value_type alpha) : m_pnts(std::move(points)) { - size_t num_pnts = m_pnts.size(); + std::size_t num_pnts = m_pnts.size(); //std::cout << "Number of points = " << num_pnts << "\n"; if (num_pnts < 4) { @@ -90,7 +110,7 @@ catmull_rom::catmull_rom(std::vector&& points, bool closed, typena m_pnts[num_pnts+2] = m_pnts[1]; auto tmp = m_pnts[num_pnts-1]; - for (int64_t i = num_pnts-1; i >= 0; --i) + for (std::ptrdiff_t i = num_pnts-1; i >= 0; --i) { m_pnts[i+1] = m_pnts[i]; } @@ -122,10 +142,10 @@ catmull_rom::catmull_rom(std::vector&& points, bool closed, typena } -template -Point catmull_rom::operator()(const typename Point::value_type s) const +template +Point catmull_rom::operator()(const typename Point::value_type s) const { - using std::size; + using std_workaround::size; if (s < 0 || s > m_max_s) { throw std::domain_error("Parameter outside bounds."); @@ -182,10 +202,10 @@ Point catmull_rom::operator()(const typename Point::value_type s) const return B1_or_C; } -template -Point catmull_rom::prime(const typename Point::value_type s) const +template +Point catmull_rom::prime(const typename Point::value_type s) const { - using std::size; + using std_workaround::size; // https://math.stackexchange.com/questions/843595/how-can-i-calculate-the-derivative-of-a-catmull-rom-spline-with-nonuniform-param // http://denkovacs.com/2016/02/catmull-rom-spline-derivatives/ if (s < 0 || s > m_max_s) diff --git a/boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp b/boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp new file mode 100644 index 0000000000..044ac72179 --- /dev/null +++ b/boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp @@ -0,0 +1,206 @@ +// Copyright Nick Thompson, 2019 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP +#include + +namespace boost{ namespace math{ namespace interpolators{ namespace detail{ + +template +Real b2_spline(Real x) { + using std::abs; + Real absx = abs(x); + if (absx < 1/Real(2)) + { + Real y = absx - 1/Real(2); + Real z = absx + 1/Real(2); + return (2-y*y-z*z)/2; + } + if (absx < Real(3)/Real(2)) + { + Real y = absx - Real(3)/Real(2); + return y*y/2; + } + return (Real) 0; +} + +template +Real b2_spline_prime(Real x) { + if (x < 0) { + return -b2_spline_prime(-x); + } + + if (x < 1/Real(2)) + { + return -2*x; + } + if (x < Real(3)/Real(2)) + { + return x - Real(3)/Real(2); + } + return (Real) 0; +} + + +template +class cardinal_quadratic_b_spline_detail +{ +public: + // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them. + // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1). + + cardinal_quadratic_b_spline_detail(const Real* const y, + size_t n, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), + Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()) + { + if (h <= 0) { + throw std::logic_error("Spacing must be > 0."); + } + m_inv_h = 1/h; + m_t0 = t0; + + if (n < 3) { + throw std::logic_error("The interpolator requires at least 3 points."); + } + + using std::isnan; + Real a; + if (isnan(left_endpoint_derivative)) { + // http://web.media.mit.edu/~crtaylor/calculator.html + a = -3*y[0] + 4*y[1] - y[2]; + } + else { + a = 2*h*left_endpoint_derivative; + } + + Real b; + if (isnan(right_endpoint_derivative)) { + b = 3*y[n-1] - 4*y[n-2] + y[n-3]; + } + else { + b = 2*h*right_endpoint_derivative; + } + + m_alpha.resize(n + 2); + + // Begin row reduction: + std::vector rhs(n + 2, std::numeric_limits::quiet_NaN()); + std::vector super_diagonal(n + 2, std::numeric_limits::quiet_NaN()); + + rhs[0] = -a; + rhs[rhs.size() - 1] = b; + + super_diagonal[0] = 0; + + for(size_t i = 1; i < rhs.size() - 1; ++i) { + rhs[i] = 8*y[i - 1]; + super_diagonal[i] = 1; + } + + // Patch up 5-diagonal problem: + rhs[1] = (rhs[1] - rhs[0])/6; + super_diagonal[1] = Real(1)/Real(3); + // First two rows are now: + // 1 0 -1 | -2hy0' + // 0 1 1/3| (8y0+2hy0')/6 + + + // Start traditional tridiagonal row reduction: + for (size_t i = 2; i < rhs.size() - 1; ++i) { + Real diagonal = 6 - super_diagonal[i - 1]; + rhs[i] = (rhs[i] - rhs[i - 1])/diagonal; + super_diagonal[i] /= diagonal; + } + + // 1 sd[n-1] 0 | rhs[n-1] + // 0 1 sd[n] | rhs[n] + // -1 0 1 | rhs[n+1] + + rhs[n+1] = rhs[n+1] + rhs[n-1]; + Real bottom_subdiagonal = super_diagonal[n-1]; + + // We're here: + // 1 sd[n-1] 0 | rhs[n-1] + // 0 1 sd[n] | rhs[n] + // 0 bs 1 | rhs[n+1] + + rhs[n+1] = (rhs[n+1]-bottom_subdiagonal*rhs[n])/(1-bottom_subdiagonal*super_diagonal[n]); + + m_alpha[n+1] = rhs[n+1]; + for (size_t i = n; i > 0; --i) { + m_alpha[i] = rhs[i] - m_alpha[i+1]*super_diagonal[i]; + } + m_alpha[0] = m_alpha[2] + rhs[0]; + } + + Real operator()(Real t) const { + if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + throw std::domain_error(err_msg); + } + // Let k, γ be defined via t = t0 + kh + γh. + // Now find all j: |k-j+1+γ|< 3/2, or, in other words + // j_min = ceil((t-t0)/h - 1/2) + // j_max = floor(t-t0)/h + 5/2) + using std::floor; + using std::ceil; + Real x = (t-m_t0)*m_inv_h; + size_t j_min = ceil(x - Real(1)/Real(2)); + size_t j_max = ceil(x + Real(5)/Real(2)); + if (j_max >= m_alpha.size()) { + j_max = m_alpha.size() - 1; + } + + Real y = 0; + x += 1; + for (size_t j = j_min; j <= j_max; ++j) { + y += m_alpha[j]*detail::b2_spline(x - j); + } + return y; + } + + Real prime(Real t) const { + if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + throw std::domain_error(err_msg); + } + // Let k, γ be defined via t = t0 + kh + γh. + // Now find all j: |k-j+1+γ|< 3/2, or, in other words + // j_min = ceil((t-t0)/h - 1/2) + // j_max = floor(t-t0)/h + 5/2) + using std::floor; + using std::ceil; + Real x = (t-m_t0)*m_inv_h; + size_t j_min = ceil(x - Real(1)/Real(2)); + size_t j_max = ceil(x + Real(5)/Real(2)); + if (j_max >= m_alpha.size()) { + j_max = m_alpha.size() - 1; + } + + Real y = 0; + x += 1; + for (size_t j = j_min; j <= j_max; ++j) { + y += m_alpha[j]*detail::b2_spline_prime(x - j); + } + return y*m_inv_h; + } + + Real t_max() const { + return m_t0 + (m_alpha.size()-3)/m_inv_h; + } + +private: + std::vector m_alpha; + Real m_inv_h; + Real m_t0; +}; + +}}}} +#endif diff --git a/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp b/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp new file mode 100644 index 0000000000..2f8d177a35 --- /dev/null +++ b/boost/math/interpolators/detail/vector_barycentric_rational_detail.hpp @@ -0,0 +1,193 @@ +/* + * Copyright Nick Thompson, 2019 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_DETAIL_HPP + +#include +#include // for std::move +#include + +namespace boost{ namespace math{ namespace detail{ + +template +class vector_barycentric_rational_imp +{ +public: + using Real = typename TimeContainer::value_type; + using Point = typename SpaceContainer::value_type; + + vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order); + + void operator()(Point& p, Real t) const; + + void eval_with_prime(Point& x, Point& dxdt, Real t) const; + + // The barycentric weights are only interesting to the unit tests: + Real weight(size_t i) const { return w_[i]; } + +private: + + void calculate_weights(size_t approximation_order); + + TimeContainer t_; + SpaceContainer y_; + TimeContainer w_; +}; + +template +vector_barycentric_rational_imp::vector_barycentric_rational_imp(TimeContainer&& t, SpaceContainer&& y, size_t approximation_order) +{ + using Real = typename TimeContainer::value_type; + using std::numeric_limits; + t_ = std::move(t); + y_ = std::move(y); + + BOOST_ASSERT_MSG(t_.size() == y_.size(), "There must be the same number of time points as space points."); + BOOST_ASSERT_MSG(approximation_order < y_.size(), "Approximation order must be < data length."); + for (size_t i = 1; i < t_.size(); ++i) + { + BOOST_ASSERT_MSG(t_[i] - t_[i-1] > (numeric_limits::min)(), "The abscissas must be listed in strictly increasing order t[0] < t[1] < ... < t[n-1]."); + } + calculate_weights(approximation_order); +} + + +template +void vector_barycentric_rational_imp::calculate_weights(size_t approximation_order) +{ + using Real = typename TimeContainer::value_type; + using std::abs; + int64_t n = t_.size(); + w_.resize(n, Real(0)); + for(int64_t k = 0; k < n; ++k) + { + int64_t i_min = (std::max)(k - (int64_t) approximation_order, (int64_t) 0); + int64_t i_max = k; + if (k >= n - (std::ptrdiff_t)approximation_order) + { + i_max = n - approximation_order - 1; + } + + for(int64_t i = i_min; i <= i_max; ++i) + { + Real inv_product = 1; + int64_t j_max = (std::min)(static_cast(i + approximation_order), static_cast(n - 1)); + for(int64_t j = i; j <= j_max; ++j) + { + if (j == k) + { + continue; + } + Real diff = t_[k] - t_[j]; + inv_product *= diff; + } + if (i % 2 == 0) + { + w_[k] += 1/inv_product; + } + else + { + w_[k] -= 1/inv_product; + } + } + } +} + + +template +void vector_barycentric_rational_imp::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const +{ + using Real = typename TimeContainer::value_type; + for (auto & x : p) + { + x = Real(0); + } + Real denominator = 0; + for(size_t i = 0; i < t_.size(); ++i) + { + // See associated commentary in the scalar version of this function. + if (t == t_[i]) + { + p = y_[i]; + return; + } + Real x = w_[i]/(t - t_[i]); + for (decltype(p.size()) j = 0; j < p.size(); ++j) + { + p[j] += x*y_[i][j]; + } + denominator += x; + } + for (decltype(p.size()) j = 0; j < p.size(); ++j) + { + p[j] /= denominator; + } + return; +} + +template +void vector_barycentric_rational_imp::eval_with_prime(typename SpaceContainer::value_type& x, typename SpaceContainer::value_type& dxdt, typename TimeContainer::value_type t) const +{ + using Point = typename SpaceContainer::value_type; + using Real = typename TimeContainer::value_type; + this->operator()(x, t); + Point numerator; + for (decltype(x.size()) i = 0; i < x.size(); ++i) + { + numerator[i] = 0; + } + Real denominator = 0; + for(decltype(t_.size()) i = 0; i < t_.size(); ++i) + { + if (t == t_[i]) + { + Point sum; + for (decltype(x.size()) i = 0; i < x.size(); ++i) + { + sum[i] = 0; + } + + for (decltype(t_.size()) j = 0; j < t_.size(); ++j) + { + if (j == i) + { + continue; + } + for (decltype(sum.size()) k = 0; k < sum.size(); ++k) + { + sum[k] += w_[j]*(y_[i][k] - y_[j][k])/(t_[i] - t_[j]); + } + } + for (decltype(sum.size()) k = 0; k < sum.size(); ++k) + { + dxdt[k] = -sum[k]/w_[i]; + } + return; + } + Real tw = w_[i]/(t - t_[i]); + Point diff; + for (decltype(diff.size()) j = 0; j < diff.size(); ++j) + { + diff[j] = (x[j] - y_[i][j])/(t-t_[i]); + } + for (decltype(diff.size()) j = 0; j < diff.size(); ++j) + { + numerator[j] += tw*diff[j]; + } + denominator += tw; + } + + for (decltype(dxdt.size()) j = 0; j < dxdt.size(); ++j) + { + dxdt[j] = numerator[j]/denominator; + } + return; +} + +}}} +#endif diff --git a/boost/math/interpolators/detail/whittaker_shannon_detail.hpp b/boost/math/interpolators/detail/whittaker_shannon_detail.hpp new file mode 100644 index 0000000000..3d4ace397e --- /dev/null +++ b/boost/math/interpolators/detail/whittaker_shannon_detail.hpp @@ -0,0 +1,126 @@ +// Copyright Nick Thompson, 2019 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_DETAIL_HPP +#include +#include +#include +#include + +namespace boost { namespace math { namespace interpolators { namespace detail { + +template +class whittaker_shannon_detail { +public: + + using Real = typename RandomAccessContainer::value_type; + whittaker_shannon_detail(RandomAccessContainer&& y, Real const & t0, Real const & h) : m_y{std::move(y)}, m_t0{t0}, m_h{h} + { + for (size_t i = 1; i < m_y.size(); i += 2) + { + m_y[i] = -m_y[i]; + } + } + + inline Real operator()(Real t) const { + using boost::math::constants::pi; + using std::isfinite; + using std::floor; + Real y = 0; + Real x = (t - m_t0)/m_h; + Real z = x; + auto it = m_y.begin(); + + // For some reason, neither clang nor g++ will cache the address of m_y.end() in a register. + // Hence make a copy of it: + auto end = m_y.end(); + while(it != end) + { + + y += *it++/z; + z -= 1; + } + + if (!isfinite(y)) + { + BOOST_ASSERT_MSG(floor(x) == ceil(x), "Floor and ceiling should be equal.\n"); + size_t i = static_cast(floor(x)); + if (i & 1) + { + return -m_y[i]; + } + return m_y[i]; + } + return y*boost::math::sin_pi(x)/pi(); + } + + Real prime(Real t) const { + using boost::math::constants::pi; + using std::isfinite; + using std::floor; + + Real x = (t - m_t0)/m_h; + if (ceil(x) == x) { + Real s = 0; + long j = static_cast(x); + long n = m_y.size(); + for (long i = 0; i < n; ++i) + { + if (j - i != 0) + { + s += m_y[i]/(j-i); + } + // else derivative of sinc at zero is zero. + } + if (j & 1) { + s /= -m_h; + } else { + s /= m_h; + } + return s; + } + Real z = x; + auto it = m_y.begin(); + Real cospix = boost::math::cos_pi(x); + Real sinpix_div_pi = boost::math::sin_pi(x)/pi(); + + Real s = 0; + auto end = m_y.end(); + while(it != end) + { + s += (*it++)*(z*cospix - sinpix_div_pi)/(z*z); + z -= 1; + } + + return s/m_h; + } + + + + Real operator[](size_t i) const { + if (i & 1) + { + return -m_y[i]; + } + return m_y[i]; + } + + RandomAccessContainer&& return_data() { + for (size_t i = 1; i < m_y.size(); i += 2) + { + m_y[i] = -m_y[i]; + } + return std::move(m_y); + } + + +private: + RandomAccessContainer m_y; + Real m_t0; + Real m_h; +}; +}}}} +#endif diff --git a/boost/math/interpolators/vector_barycentric_rational.hpp b/boost/math/interpolators/vector_barycentric_rational.hpp new file mode 100644 index 0000000000..8289799bc9 --- /dev/null +++ b/boost/math/interpolators/vector_barycentric_rational.hpp @@ -0,0 +1,82 @@ +/* + * Copyright Nick Thompson, 2019 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + * + * Exactly the same as barycentric_rational.hpp, but delivers values in $\mathbb{R}^n$. + * In some sense this is trivial, since each component of the vector is computed in exactly the same + * as would be computed by barycentric_rational.hpp. But this is a bit more efficient and convenient. + */ + +#ifndef BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_HPP +#define BOOST_MATH_INTERPOLATORS_VECTOR_BARYCENTRIC_RATIONAL_HPP + +#include +#include + +namespace boost{ namespace math{ + +template +class vector_barycentric_rational +{ +public: + using Real = typename TimeContainer::value_type; + using Point = typename SpaceContainer::value_type; + vector_barycentric_rational(TimeContainer&& times, SpaceContainer&& points, size_t approximation_order = 3); + + void operator()(Point& x, Real t) const; + + // I have validated using google benchmark that returning a value is no more expensive populating it, + // at least for Eigen vectors with known size at compile-time. + // This is kinda a weird thing to discover since it goes against the advice of basically every high-performance computing book. + Point operator()(Real t) const { + Point p; + this->operator()(p, t); + return p; + } + + void prime(Point& dxdt, Real t) const { + Point x; + m_imp->eval_with_prime(x, dxdt, t); + } + + Point prime(Real t) const { + Point p; + this->prime(p, t); + return p; + } + + void eval_with_prime(Point& x, Point& dxdt, Real t) const { + m_imp->eval_with_prime(x, dxdt, t); + return; + } + + std::pair eval_with_prime(Real t) const { + Point x; + Point dxdt; + m_imp->eval_with_prime(x, dxdt, t); + return {x, dxdt}; + } + +private: + std::shared_ptr> m_imp; +}; + + +template +vector_barycentric_rational::vector_barycentric_rational(TimeContainer&& times, SpaceContainer&& points, size_t approximation_order): + m_imp(std::make_shared>(std::move(times), std::move(points), approximation_order)) +{ + return; +} + +template +void vector_barycentric_rational::operator()(typename SpaceContainer::value_type& p, typename TimeContainer::value_type t) const +{ + m_imp->operator()(p, t); + return; +} + +}} +#endif diff --git a/boost/math/interpolators/whittaker_shannon.hpp b/boost/math/interpolators/whittaker_shannon.hpp new file mode 100644 index 0000000000..a5c5367a9d --- /dev/null +++ b/boost/math/interpolators/whittaker_shannon.hpp @@ -0,0 +1,47 @@ +// Copyright Nick Thompson, 2019 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP +#define BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP +#include +#include + +namespace boost { namespace math { namespace interpolators { + +template +class whittaker_shannon { +public: + + using Real = typename RandomAccessContainer::value_type; + whittaker_shannon(RandomAccessContainer&& y, Real const & t0, Real const & h) + : m_impl(std::make_shared>(std::move(y), t0, h)) + {} + + inline Real operator()(Real t) const + { + return m_impl->operator()(t); + } + + inline Real prime(Real t) const + { + return m_impl->prime(t); + } + + inline Real operator[](size_t i) const + { + return m_impl->operator[](i); + } + + RandomAccessContainer&& return_data() + { + return m_impl->return_data(); + } + + +private: + std::shared_ptr> m_impl; +}; +}}} +#endif diff --git a/boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp b/boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp new file mode 100644 index 0000000000..a449c29b58 --- /dev/null +++ b/boost/math/quadrature/detail/ooura_fourier_integrals_detail.hpp @@ -0,0 +1,651 @@ +// Copyright Nick Thompson, 2019 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef BOOST_MATH_QUADRATURE_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP +#define BOOST_MATH_QUADRATURE_DETAIL_OOURA_FOURIER_INTEGRALS_DETAIL_HPP +#include // for std::pair. +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { namespace math { namespace quadrature { namespace detail { + +// Ooura and Mori, A robust double exponential formula for Fourier-type integrals, +// eta is the argument to the exponential in equation 3.3: +template +std::pair ooura_eta(Real x, Real alpha) { + using std::expm1; + using std::exp; + using std::abs; + Real expx = exp(x); + Real eta_prime = 2 + alpha/expx + expx/4; + Real eta; + // This is the fast branch: + if (abs(x) > 0.125) { + eta = 2*x - alpha*(1/expx - 1) + (expx - 1)/4; + } + else {// this is the slow branch using expm1 for small x: + eta = 2*x - alpha*expm1(-x) + expm1(x)/4; + } + return {eta, eta_prime}; +} + +// Ooura and Mori, A robust double exponential formula for Fourier-type integrals, +// equation 3.6: +template +Real calculate_ooura_alpha(Real h) +{ + using boost::math::constants::pi; + using std::log1p; + using std::sqrt; + Real x = sqrt(16 + 4*log1p(pi()/h)/h); + return 1/x; +} + +template +std::pair ooura_sin_node_and_weight(long n, Real h, Real alpha) +{ + using std::expm1; + using std::exp; + using std::abs; + using boost::math::constants::pi; + using std::isnan; + + if (n == 0) { + // Equation 44 of https://arxiv.org/pdf/0911.4796.pdf + // Fourier Transform of the Stretched Exponential Function: Analytic Error Bounds, + // Double Exponential Transform, and Open-Source Implementation, + // Joachim Wuttke, + // The C library libkww provides functions to compute the Kohlrausch-Williams-Watts function, + // the Laplace-Fourier transform of the stretched (or compressed) exponential function exp(-t^beta) + // for exponent beta between 0.1 and 1.9 with sixteen decimal digits accuracy. + + Real eta_prime_0 = Real(2) + alpha + Real(1)/Real(4); + Real node = pi()/(eta_prime_0*h); + Real weight = pi()*boost::math::sin_pi(1/(eta_prime_0*h)); + Real eta_dbl_prime = -alpha + Real(1)/Real(4); + Real phi_prime_0 = (1 - eta_dbl_prime/(eta_prime_0*eta_prime_0))/2; + weight *= phi_prime_0; + return {node, weight}; + } + Real x = n*h; + auto p = ooura_eta(x, alpha); + auto eta = p.first; + auto eta_prime = p.second; + + Real expm1_meta = expm1(-eta); + Real exp_meta = exp(-eta); + Real node = -n*pi()/expm1_meta; + + + // I have verified that this is not a significant source of inaccuracy in the weight computation: + Real phi_prime = -(expm1_meta + x*exp_meta*eta_prime)/(expm1_meta*expm1_meta); + + // The main source of inaccuracy is in computation of sin_pi. + // But I've agonized over this, and I think it's as good as it can get: + Real s = pi(); + Real arg; + if(eta > 1) { + arg = n/( 1/exp_meta - 1 ); + s *= boost::math::sin_pi(arg); + if (n&1) { + s *= -1; + } + } + else if (eta < -1) { + arg = n/(1-exp_meta); + s *= boost::math::sin_pi(arg); + } + else { + arg = -n*exp_meta/expm1_meta; + s *= boost::math::sin_pi(arg); + if (n&1) { + s *= -1; + } + } + + Real weight = s*phi_prime; + return {node, weight}; +} + +#ifdef BOOST_MATH_INSTRUMENT_OOURA +template +void print_ooura_estimate(size_t i, Real I0, Real I1, Real omega) { + using std::abs; + std::cout << std::defaultfloat + << std::setprecision(std::numeric_limits::digits10) + << std::fixed; + std::cout << "h = " << Real(1)/Real(1<::epsilon() * 2) { + throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff."); + } + using std::abs; + requested_levels_ = levels; + starting_level_ = 0; + rel_err_goal_ = relative_error_goal; + big_nodes_.reserve(levels); + bweights_.reserve(levels); + little_nodes_.reserve(levels); + lweights_.reserve(levels); + + for (size_t i = 0; i < levels; ++i) { + if (std::is_same::value) { + add_level(i); + } + else if (std::is_same::value) { + add_level(i); + } + else { + add_level(i); + } + } + } + + std::vector> const & big_nodes() const { + return big_nodes_; + } + + std::vector> const & weights_for_big_nodes() const { + return bweights_; + } + + std::vector> const & little_nodes() const { + return little_nodes_; + } + + std::vector> const & weights_for_little_nodes() const { + return lweights_; + } + + template + std::pair integrate(F const & f, Real omega) { + using std::abs; + using std::max; + using boost::math::constants::pi; + + if (omega == 0) { + return {Real(0), Real(0)}; + } + if (omega < 0) { + auto p = this->integrate(f, -omega); + return {-p.first, p.second}; + } + + Real I1 = std::numeric_limits::quiet_NaN(); + Real relative_error_estimate = std::numeric_limits::quiet_NaN(); + // As we compute integrals, we learn about their structure. + // Assuming we compute f(t)sin(wt) for many different omega, this gives some + // a posteriori ability to choose a refinement level that is roughly appropriate. + size_t i = starting_level_; + do { + Real I0 = estimate_integral(f, omega, i); +#ifdef BOOST_MATH_INSTRUMENT_OOURA + print_ooura_estimate(i, I0, I1, omega); +#endif + Real absolute_error_estimate = abs(I0-I1); + Real scale = max(abs(I0), abs(I1)); + if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) { + starting_level_ = max(long(i) - 1, long(0)); + return {I0/omega, absolute_error_estimate/scale}; + } + I1 = I0; + } while(++i < big_nodes_.size()); + + // We've used up all our precomputed levels. + // Now we need to add more. + // It might seems reasonable to just keep adding levels indefinitely, if that's what the user wants. + // But in fact the nodes and weights just merge into each other and the error gets worse after a certain number. + // This value for max_additional_levels was chosen by observation of a slowly converging oscillatory integral: + // f(x) := cos(7cos(x))sin(x)/x + size_t max_additional_levels = 4; + while (big_nodes_.size() < requested_levels_ + max_additional_levels) { + size_t ii = big_nodes_.size(); + if (std::is_same::value) { + add_level(ii); + } + else if (std::is_same::value) { + add_level(ii); + } + else { + add_level(ii); + } + Real I0 = estimate_integral(f, omega, ii); + Real absolute_error_estimate = abs(I0-I1); + Real scale = max(abs(I0), abs(I1)); +#ifdef BOOST_MATH_INSTRUMENT_OOURA + print_ooura_estimate(ii, I0, I1, omega); +#endif + if (absolute_error_estimate <= rel_err_goal_*scale) { + starting_level_ = max(long(ii) - 1, long(0)); + return {I0/omega, absolute_error_estimate/scale}; + } + I1 = I0; + ++ii; + } + + starting_level_ = static_cast(big_nodes_.size() - 2); + return {I1/omega, relative_error_estimate}; + } + +private: + + template + void add_level(size_t i) { + using std::abs; + size_t current_num_levels = big_nodes_.size(); + Real unit_roundoff = std::numeric_limits::epsilon()/2; + // h0 = 1. Then all further levels have h_i = 1/2^i. + // Since the nodes don't nest, we could conceivably divide h by (say) 1.5, or 3. + // It's not clear how much benefit (or loss) would be obtained from this. + PreciseReal h = PreciseReal(1)/PreciseReal(1< bnode_row; + std::vector bweight_row; + + // This is a pretty good estimate for how many elements will be placed in the vector: + bnode_row.reserve((static_cast(1)<(1)< lnode_row; + std::vector lweight_row; + + lnode_row.reserve((static_cast(1)<(1)<(precise_nw.first); + Real weight = static_cast(precise_nw.second); + w = weight; + if (bnode_row.size() == bnode_row.capacity()) { + bnode_row.reserve(2*bnode_row.size()); + bweight_row.reserve(2*bnode_row.size()); + } + + bnode_row.push_back(node); + bweight_row.push_back(weight); + if (abs(weight) > max_weight) { + max_weight = abs(weight); + } + ++n; + // f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff. + } while(abs(w) > unit_roundoff*max_weight); + + // This class tends to consume a lot of memory; shrink the vectors back down to size: + bnode_row.shrink_to_fit(); + bweight_row.shrink_to_fit(); + // Why we are splitting the nodes into regimes where t_n >> 1 and t_n << 1? + // It will create the opportunity to sensibly truncate the quadrature sum to significant terms. + n = -1; + do { + auto precise_nw = ooura_sin_node_and_weight(n, h, alpha); + Real node = static_cast(precise_nw.first); + if (node <= 0) { + break; + } + Real weight = static_cast(precise_nw.second); + w = weight; + using std::isnan; + if (isnan(node)) { + // This occurs at n = -11 in quad precision: + break; + } + if (lnode_row.size() > 0) { + if (lnode_row[lnode_row.size()-1] == node) { + // The nodes have fused into each other: + break; + } + } + if (lnode_row.size() == lnode_row.capacity()) { + lnode_row.reserve(2*lnode_row.size()); + lweight_row.reserve(2*lnode_row.size()); + } + lnode_row.push_back(node); + lweight_row.push_back(weight); + if (abs(weight) > max_weight) { + max_weight = abs(weight); + } + --n; + // f(t)->infty is possible as t->0, hence compute up to the min. + } while(abs(w) > (std::numeric_limits::min)()*max_weight); + + lnode_row.shrink_to_fit(); + lweight_row.shrink_to_fit(); + + // std::scoped_lock once C++17 is more common? + std::lock_guard lock(node_weight_mutex_); + // Another thread might have already finished this calculation and appended it to the nodes/weights: + if (current_num_levels == big_nodes_.size()) { + big_nodes_.push_back(bnode_row); + bweights_.push_back(bweight_row); + + little_nodes_.push_back(lnode_row); + lweights_.push_back(lweight_row); + } + } + + template + Real estimate_integral(F const & f, Real omega, size_t i) { + // Because so few function evaluations are required to get high accuracy on the integrals in the tests, + // Kahan summation doesn't really help. + //auto cond = boost::math::tools::summation_condition_number(0); + Real I0 = 0; + auto const & b_nodes = big_nodes_[i]; + auto const & b_weights = bweights_[i]; + // Will benchmark if this is helpful: + Real inv_omega = 1/omega; + for(size_t j = 0 ; j < b_nodes.size(); ++j) { + I0 += f(b_nodes[j]*inv_omega)*b_weights[j]; + } + + auto const & l_nodes = little_nodes_[i]; + auto const & l_weights = lweights_[i]; + // If f decays rapidly as |t|->infty, not all of these calls are necessary. + for (size_t j = 0; j < l_nodes.size(); ++j) { + I0 += f(l_nodes[j]*inv_omega)*l_weights[j]; + } + return I0; + } + + std::mutex node_weight_mutex_; + // Nodes for n >= 0, giving t_n = pi*phi(nh)/h. Generally t_n >> 1. + std::vector> big_nodes_; + // The term bweights_ will indicate that these are weights corresponding + // to the big nodes: + std::vector> bweights_; + + // Nodes for n < 0: Generally t_n << 1, and an invariant is that t_n > 0. + std::vector> little_nodes_; + std::vector> lweights_; + Real rel_err_goal_; + std::atomic starting_level_; + size_t requested_levels_; +}; + +template +class ooura_fourier_cos_detail { +public: + ooura_fourier_cos_detail(const Real relative_error_goal, size_t levels) { +#ifdef BOOST_MATH_INSTRUMENT_OOURA + std::cout << "ooura_fourier_cos with relative error goal " << relative_error_goal + << " & " << levels << " levels." << std::endl; + std::cout << "epsilon for type = " << std::numeric_limits::epsilon() << std::endl; +#endif // BOOST_MATH_INSTRUMENT_OOURA + if (relative_error_goal < std::numeric_limits::epsilon() * 2) { + throw std::domain_error("The relative error goal cannot be smaller than the unit roundoff!"); + } + + using std::abs; + requested_levels_ = levels; + starting_level_ = 0; + rel_err_goal_ = relative_error_goal; + big_nodes_.reserve(levels); + bweights_.reserve(levels); + little_nodes_.reserve(levels); + lweights_.reserve(levels); + + for (size_t i = 0; i < levels; ++i) { + if (std::is_same::value) { + add_level(i); + } + else if (std::is_same::value) { + add_level(i); + } + else { + add_level(i); + } + } + + } + + template + std::pair integrate(F const & f, Real omega) { + using std::abs; + using std::max; + using boost::math::constants::pi; + + if (omega == 0) { + throw std::domain_error("At omega = 0, the integral is not oscillatory. The user must choose an appropriate method for this case.\n"); + } + + if (omega < 0) { + return this->integrate(f, -omega); + } + + Real I1 = std::numeric_limits::quiet_NaN(); + Real absolute_error_estimate = std::numeric_limits::quiet_NaN(); + Real scale = std::numeric_limits::quiet_NaN(); + size_t i = starting_level_; + do { + Real I0 = estimate_integral(f, omega, i); +#ifdef BOOST_MATH_INSTRUMENT_OOURA + print_ooura_estimate(i, I0, I1, omega); +#endif + absolute_error_estimate = abs(I0-I1); + scale = max(abs(I0), abs(I1)); + if (!isnan(I1) && absolute_error_estimate <= rel_err_goal_*scale) { + starting_level_ = max(long(i) - 1, long(0)); + return {I0/omega, absolute_error_estimate/scale}; + } + I1 = I0; + } while(++i < big_nodes_.size()); + + size_t max_additional_levels = 4; + while (big_nodes_.size() < requested_levels_ + max_additional_levels) { + size_t ii = big_nodes_.size(); + if (std::is_same::value) { + add_level(ii); + } + else if (std::is_same::value) { + add_level(ii); + } + else { + add_level(ii); + } + Real I0 = estimate_integral(f, omega, ii); +#ifdef BOOST_MATH_INSTRUMENT_OOURA + print_ooura_estimate(ii, I0, I1, omega); +#endif + absolute_error_estimate = abs(I0-I1); + scale = max(abs(I0), abs(I1)); + if (absolute_error_estimate <= rel_err_goal_*scale) { + starting_level_ = max(long(ii) - 1, long(0)); + return {I0/omega, absolute_error_estimate/scale}; + } + I1 = I0; + ++ii; + } + + starting_level_ = static_cast(big_nodes_.size() - 2); + return {I1/omega, absolute_error_estimate/scale}; + } + +private: + + template + void add_level(size_t i) { + using std::abs; + size_t current_num_levels = big_nodes_.size(); + Real unit_roundoff = std::numeric_limits::epsilon()/2; + PreciseReal h = PreciseReal(1)/PreciseReal(1< bnode_row; + std::vector bweight_row; + bnode_row.reserve((static_cast(1)<(1)< lnode_row; + std::vector lweight_row; + + lnode_row.reserve((static_cast(1)<(1)<(precise_nw.first); + Real weight = static_cast(precise_nw.second); + w = weight; + if (bnode_row.size() == bnode_row.capacity()) { + bnode_row.reserve(2*bnode_row.size()); + bweight_row.reserve(2*bnode_row.size()); + } + + bnode_row.push_back(node); + bweight_row.push_back(weight); + if (abs(weight) > max_weight) { + max_weight = abs(weight); + } + ++n; + // f(t)->0 as t->infty, which is why the weights are computed up to the unit roundoff. + } while(abs(w) > unit_roundoff*max_weight); + + bnode_row.shrink_to_fit(); + bweight_row.shrink_to_fit(); + n = -1; + do { + auto precise_nw = ooura_cos_node_and_weight(n, h, alpha); + Real node = static_cast(precise_nw.first); + // The function cannot be singular at zero, + // so zero is not a unreasonable node, + // unlike in the case of the Fourier Sine. + // Hence only break if the node is negative. + if (node < 0) { + break; + } + Real weight = static_cast(precise_nw.second); + w = weight; + if (lnode_row.size() > 0) { + if (lnode_row.back() == node) { + // The nodes have fused into each other: + break; + } + } + if (lnode_row.size() == lnode_row.capacity()) { + lnode_row.reserve(2*lnode_row.size()); + lweight_row.reserve(2*lnode_row.size()); + } + + lnode_row.push_back(node); + lweight_row.push_back(weight); + if (abs(weight) > max_weight) { + max_weight = abs(weight); + } + --n; + } while(abs(w) > (std::numeric_limits::min)()*max_weight); + + lnode_row.shrink_to_fit(); + lweight_row.shrink_to_fit(); + + std::lock_guard lock(node_weight_mutex_); + // Another thread might have already finished this calculation and appended it to the nodes/weights: + if (current_num_levels == big_nodes_.size()) { + big_nodes_.push_back(bnode_row); + bweights_.push_back(bweight_row); + + little_nodes_.push_back(lnode_row); + lweights_.push_back(lweight_row); + } + } + + template + Real estimate_integral(F const & f, Real omega, size_t i) { + Real I0 = 0; + auto const & b_nodes = big_nodes_[i]; + auto const & b_weights = bweights_[i]; + Real inv_omega = 1/omega; + for(size_t j = 0 ; j < b_nodes.size(); ++j) { + I0 += f(b_nodes[j]*inv_omega)*b_weights[j]; + } + + auto const & l_nodes = little_nodes_[i]; + auto const & l_weights = lweights_[i]; + for (size_t j = 0; j < l_nodes.size(); ++j) { + I0 += f(l_nodes[j]*inv_omega)*l_weights[j]; + } + return I0; + } + + std::mutex node_weight_mutex_; + std::vector> big_nodes_; + std::vector> bweights_; + + std::vector> little_nodes_; + std::vector> lweights_; + Real rel_err_goal_; + std::atomic starting_level_; + size_t requested_levels_; +}; + + +}}}} +#endif diff --git a/boost/math/quadrature/ooura_fourier_integrals.hpp b/boost/math/quadrature/ooura_fourier_integrals.hpp new file mode 100644 index 0000000000..b31996c2bc --- /dev/null +++ b/boost/math/quadrature/ooura_fourier_integrals.hpp @@ -0,0 +1,68 @@ +// Copyright Nick Thompson, 2019 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +/* + * References: + * Ooura, Takuya, and Masatake Mori. "A robust double exponential formula for Fourier-type integrals." Journal of computational and applied mathematics 112.1-2 (1999): 229-241. + * http://www.kurims.kyoto-u.ac.jp/~ooura/intde.html + */ +#ifndef BOOST_MATH_QUADRATURE_OOURA_FOURIER_INTEGRALS_HPP +#define BOOST_MATH_QUADRATURE_OOURA_FOURIER_INTEGRALS_HPP +#include +#include + +namespace boost { namespace math { namespace quadrature { + +template +class ooura_fourier_sin { +public: + ooura_fourier_sin(const Real relative_error_tolerance = tools::root_epsilon(), size_t levels = sizeof(Real)) : impl_(std::make_shared>(relative_error_tolerance, levels)) + {} + + template + std::pair integrate(F const & f, Real omega) { + return impl_->integrate(f, omega); + } + + // These are just for debugging/unit tests: + std::vector> const & big_nodes() const { + return impl_->big_nodes(); + } + + std::vector> const & weights_for_big_nodes() const { + return impl_->weights_for_big_nodes(); + } + + std::vector> const & little_nodes() const { + return impl_->little_nodes(); + } + + std::vector> const & weights_for_little_nodes() const { + return impl_->weights_for_little_nodes(); + } + +private: + std::shared_ptr> impl_; +}; + + +template +class ooura_fourier_cos { +public: + ooura_fourier_cos(const Real relative_error_tolerance = tools::root_epsilon(), size_t levels = sizeof(Real)) : impl_(std::make_shared>(relative_error_tolerance, levels)) + {} + + template + std::pair integrate(F const & f, Real omega) { + return impl_->integrate(f, omega); + } +private: + std::shared_ptr> impl_; +}; + + +}}} +#endif diff --git a/boost/math/special_functions/beta.hpp b/boost/math/special_functions/beta.hpp index 35b114ef15..c4acbcb4b0 100644 --- a/boost/math/special_functions/beta.hpp +++ b/boost/math/special_functions/beta.hpp @@ -790,8 +790,11 @@ struct Pn_size { // This is likely to be enough for ~35-50 digit accuracy // but it's hard to quantify exactly: - BOOST_STATIC_CONSTANT(unsigned, value = 50); - BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= 100); + BOOST_STATIC_CONSTANT(unsigned, value = + ::boost::math::max_factorial::value >= 100 ? 50 + : ::boost::math::max_factorial::value >= ::boost::math::max_factorial::value ? 30 + : ::boost::math::max_factorial::value >= ::boost::math::max_factorial::value ? 15 : 1); + BOOST_STATIC_ASSERT(::boost::math::max_factorial::value >= ::boost::math::max_factorial::value); }; template <> struct Pn_size diff --git a/boost/math/special_functions/cos_pi.hpp b/boost/math/special_functions/cos_pi.hpp index 669a2c87ae..f9aeeca993 100644 --- a/boost/math/special_functions/cos_pi.hpp +++ b/boost/math/special_functions/cos_pi.hpp @@ -33,7 +33,7 @@ T cos_pi_imp(T x, const Policy& pol) x = -x; } T rem = floor(x); - if(itrunc(rem, pol) & 1) + if(iconvert(rem, pol) & 1) invert = !invert; rem = x - rem; if(rem > 0.5f) @@ -66,7 +66,10 @@ inline typename tools::promote_args::type cos_pi(T x, const Policy&) policies::promote_float, policies::promote_double, policies::discrete_quantile<>, - policies::assert_undefined<> >::type forwarding_policy; + policies::assert_undefined<>, + // We want to igore overflows since the result is in [-1,1] and the + // check slows the code down considerably. + policies::overflow_error >::type forwarding_policy; return policies::checked_narrowing_cast(boost::math::detail::cos_pi_imp(x, forwarding_policy()), "cos_pi"); } diff --git a/boost/math/special_functions/detail/unchecked_factorial.hpp b/boost/math/special_functions/detail/unchecked_factorial.hpp index 642df7ce50..aa5b484df0 100644 --- a/boost/math/special_functions/detail/unchecked_factorial.hpp +++ b/boost/math/special_functions/detail/unchecked_factorial.hpp @@ -774,6 +774,33 @@ inline T unchecked_factorial_imp(unsigned i, const mpl::int_<0>&) return factorials[i]; } +template +inline T unchecked_factorial_imp(unsigned i, const mpl::int_::digits>&) +{ + return unchecked_factorial(i); +} + +template +inline T unchecked_factorial_imp(unsigned i, const mpl::int_::digits>&) +{ + return unchecked_factorial(i); +} + +#if DBL_MANT_DIG != LDBL_MANT_DIG +template +inline T unchecked_factorial_imp(unsigned i, const mpl::int_&) +{ + return unchecked_factorial(i); +} +#endif +#ifdef BOOST_MATH_USE_FLOAT128 +template +inline T unchecked_factorial_imp(unsigned i, const mpl::int_<113>&) +{ + return unchecked_factorial(i); +} +#endif + template inline T unchecked_factorial(unsigned i) { @@ -781,12 +808,25 @@ inline T unchecked_factorial(unsigned i) return unchecked_factorial_imp(i, tag_type()); } +#ifdef BOOST_MATH_USE_FLOAT128 +#define BOOST_MATH_DETAIL_FLOAT128_MAX_FACTORIAL : std::numeric_limits::digits == 113 ? max_factorial::value +#else +#define BOOST_MATH_DETAIL_FLOAT128_MAX_FACTORIAL +#endif + template struct max_factorial { - BOOST_STATIC_CONSTANT(unsigned, value = 100); + BOOST_STATIC_CONSTANT(unsigned, value = + std::numeric_limits::digits == std::numeric_limits::digits ? max_factorial::value + : std::numeric_limits::digits == std::numeric_limits::digits ? max_factorial::value + : std::numeric_limits::digits == std::numeric_limits::digits ? max_factorial::value + BOOST_MATH_DETAIL_FLOAT128_MAX_FACTORIAL + : 100); }; +#undef BOOST_MATH_DETAIL_FLOAT128_MAX_FACTORIAL + #else // BOOST_MATH_NO_LEXICAL_CAST template diff --git a/boost/math/special_functions/ellint_1.hpp b/boost/math/special_functions/ellint_1.hpp index d1d9d72e30..334c0ad808 100644 --- a/boost/math/special_functions/ellint_1.hpp +++ b/boost/math/special_functions/ellint_1.hpp @@ -51,12 +51,6 @@ T ellint_f_imp(T phi, T k, const Policy& pol) BOOST_MATH_INSTRUMENT_VARIABLE(k); BOOST_MATH_INSTRUMENT_VARIABLE(function); - if (abs(k) > 1) - { - return policies::raise_domain_error(function, - "Got k = %1%, function requires |k| <= 1", k, pol); - } - bool invert = false; if(phi < 0) { @@ -104,18 +98,24 @@ T ellint_f_imp(T phi, T k, const Policy& pol) } T sinp = sin(rphi); sinp *= sinp; + if (sinp * k * k >= 1) + { + return policies::raise_domain_error(function, + "Got k^2 * sin^2(phi) = %1%, but the function requires this < 1", sinp * k * k, pol); + } T cosp = cos(rphi); cosp *= cosp; BOOST_MATH_INSTRUMENT_VARIABLE(sinp); BOOST_MATH_INSTRUMENT_VARIABLE(cosp); if(sinp > tools::min_value()) { + BOOST_ASSERT(rphi != 0); // precondition, can't be true if sin(rphi) != 0. // // Use http://dlmf.nist.gov/19.25#E5, note that // c-1 simplifies to cot^2(rphi) which avoid cancellation: // T c = 1 / sinp; - result = rphi == 0 ? static_cast(0) : static_cast(s * ellint_rf_imp(T(cosp / sinp), T(c - k * k), c, pol)); + result = static_cast(s * ellint_rf_imp(T(cosp / sinp), T(c - k * k), c, pol)); } else result = s * sin(rphi); diff --git a/boost/math/special_functions/ellint_2.hpp b/boost/math/special_functions/ellint_2.hpp index 9ee6b63821..2687e9c8ce 100644 --- a/boost/math/special_functions/ellint_2.hpp +++ b/boost/math/special_functions/ellint_2.hpp @@ -49,6 +49,9 @@ T ellint_e_imp(T phi, T k, const Policy& pol) using namespace boost::math::constants; bool invert = false; + if (phi == 0) + return 0; + if(phi < 0) { phi = fabs(phi); @@ -95,11 +98,7 @@ T ellint_e_imp(T phi, T k, const Policy& pol) rphi = constants::half_pi() - rphi; } T k2 = k * k; - if(k2 > 1) - { - return policies::raise_domain_error("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol); - } - else if(rphi < tools::root_epsilon()) + if(boost::math::pow<3>(rphi) * k2 / 6 < tools::epsilon() * fabs(rphi)) { // See http://functions.wolfram.com/EllipticIntegrals/EllipticE2/06/01/03/0001/ result = s * rphi; @@ -108,6 +107,10 @@ T ellint_e_imp(T phi, T k, const Policy& pol) { // http://dlmf.nist.gov/19.25#E10 T sinp = sin(rphi); + if (k2 * sinp * sinp >= 1) + { + return policies::raise_domain_error("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol); + } T cosp = cos(rphi); T c = 1 / (sinp * sinp); T cm1 = cosp * cosp / (sinp * sinp); // c - 1 diff --git a/boost/math/special_functions/ellint_3.hpp b/boost/math/special_functions/ellint_3.hpp index b8b36729cf..dfc104514e 100644 --- a/boost/math/special_functions/ellint_3.hpp +++ b/boost/math/special_functions/ellint_3.hpp @@ -49,15 +49,15 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol) static const char* function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)"; - if(abs(k) > 1) - { - return policies::raise_domain_error(function, - "Got k = %1%, function requires |k| <= 1", k, pol); - } T sphi = sin(fabs(phi)); T result = 0; + if (k * k * sphi * sphi > 1) + { + return policies::raise_domain_error(function, + "Got k = %1%, function requires |k| <= 1", k, pol); + } // Special cases first: if(v == 0) { @@ -73,6 +73,9 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol) if(v == 1) { + if (k == 0) + return tan(phi); + // http://functions.wolfram.com/08.06.03.0008.01 T m = k * k; result = sqrt(1 - m * sphi * sphi) * tan(phi) - ellint_e_imp(phi, k, pol); @@ -143,10 +146,6 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol) T vcr = sqrt(vc); return atan(vcr * tan(phi)) / vcr; } - else if(v == 1) - { - return tan(phi); - } else { // v > 1: @@ -155,7 +154,7 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol) return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr); } } - if(v < 0) + if((v < 0) && fabs(k) <= 1) { // // If we don't shift to 0 <= v <= 1 we get diff --git a/boost/math/special_functions/ellint_d.hpp b/boost/math/special_functions/ellint_d.hpp index fa5c53db18..b482922d5d 100644 --- a/boost/math/special_functions/ellint_d.hpp +++ b/boost/math/special_functions/ellint_d.hpp @@ -91,7 +91,7 @@ T ellint_d_imp(T phi, T k, const Policy& pol) T c = 1 / (sinp * sinp); T cm1 = cosp * cosp / (sinp * sinp); // c - 1 T k2 = k * k; - if(k2 > 1) + if(k2 * sinp * sinp > 1) { return policies::raise_domain_error("boost::math::ellint_d<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol); } diff --git a/boost/math/special_functions/lambert_w.hpp b/boost/math/special_functions/lambert_w.hpp index 6b3dcebe2f..65d3be5f5c 100644 --- a/boost/math/special_functions/lambert_w.hpp +++ b/boost/math/special_functions/lambert_w.hpp @@ -61,7 +61,7 @@ BOOST_MATH_INSTRUMENT_LAMBERT_W_SMALL_Z_SERIES_ITERATIONS // Show evaluation of #include #include #include // boost::math::tools::max_value(). -#include // For create_test_value and macro BOOST_MATH_TEST_VALUE. +#include #include #include @@ -296,27 +296,27 @@ T lambert_w_singularity_series(const T p) // -T(0.000672061631156136204L), j14 //+T(1003663334225097487uLL) / 234281684403486720000uLL, // j15 0.00044247306181462090993020760858473726479232802068800 error C2177: constant too big //+T(0.000442473061814620910L, // j15 - BOOST_MATH_TEST_VALUE(T, +0.000442473061814620910), // j15 + BOOST_MATH_BIG_CONSTANT(T, 64, +0.000442473061814620910), // j15 // -T(0.000292677224729627445L), // j16 - BOOST_MATH_TEST_VALUE(T, -0.000292677224729627445), // j16 + BOOST_MATH_BIG_CONSTANT(T, 64, -0.000292677224729627445), // j16 //+T(0.000194387276054539318L), // j17 - BOOST_MATH_TEST_VALUE(T, 0.000194387276054539318), // j17 + BOOST_MATH_BIG_CONSTANT(T, 64, 0.000194387276054539318), // j17 //-T(0.000129574266852748819L), // j18 - BOOST_MATH_TEST_VALUE(T, -0.000129574266852748819), // j18 + BOOST_MATH_BIG_CONSTANT(T, 64, -0.000129574266852748819), // j18 //+T(0.0000866503580520812717L), // j19 N[+1150497127780071399782389/13277465363600276402995200000, 50] 0.000086650358052081271660451590462390293190597827783288 - BOOST_MATH_TEST_VALUE(T, +0.0000866503580520812717), // j19 + BOOST_MATH_BIG_CONSTANT(T, 64, +0.0000866503580520812717), // j19 //-T(0.0000581136075044138168L) // j20 N[2853534237182741069/49102686267859224000000, 50] 0.000058113607504413816772205464778828177256611844221913 // -T(2853534237182741069uLL) / 49102686267859224000000uLL // j20 // error C2177: constant too big, - // so must use BOOST_MATH_TEST_VALUE(T, ) format in hope of using suffix Q for quad or decimal digits string for others. + // so must use BOOST_MATH_BIG_CONSTANT(T, ) format in hope of using suffix Q for quad or decimal digits string for others. //-T(0.000058113607504413816772205464778828177256611844221913L), // j20 N[2853534237182741069/49102686267859224000000, 50] 0.000058113607504413816772205464778828177256611844221913 - BOOST_MATH_TEST_VALUE(T, -0.000058113607504413816772205464778828177256611844221913) // j20 - last used by Fukushima + BOOST_MATH_BIG_CONSTANT(T, 113, -0.000058113607504413816772205464778828177256611844221913) // j20 - last used by Fukushima // More terms don't seem to give any improvement (worse in fact) and are not use for many z values. - //BOOST_MATH_TEST_VALUE(T, +0.000039076684867439051635395583044527492132109160553593), // j21 - //BOOST_MATH_TEST_VALUE(T, -0.000026338064747231098738584082718649443078703982217219), // j22 - //BOOST_MATH_TEST_VALUE(T, +0.000017790345805079585400736282075184540383274460464169), // j23 - //BOOST_MATH_TEST_VALUE(T, -0.000012040352739559976942274116578992585158113153190354), // j24 - //BOOST_MATH_TEST_VALUE(T, +8.1635319824966121713827512573558687050675701559448E-6), // j25 - //BOOST_MATH_TEST_VALUE(T, -5.5442032085673591366657251660804575198155559225316E-6) // j26 + //BOOST_MATH_BIG_CONSTANT(T, +0.000039076684867439051635395583044527492132109160553593), // j21 + //BOOST_MATH_BIG_CONSTANT(T, -0.000026338064747231098738584082718649443078703982217219), // j22 + //BOOST_MATH_BIG_CONSTANT(T, +0.000017790345805079585400736282075184540383274460464169), // j23 + //BOOST_MATH_BIG_CONSTANT(T, -0.000012040352739559976942274116578992585158113153190354), // j24 + //BOOST_MATH_BIG_CONSTANT(T, +8.1635319824966121713827512573558687050675701559448E-6), // j25 + //BOOST_MATH_BIG_CONSTANT(T, -5.5442032085673591366657251660804575198155559225316E-6) // j26 // -T(5.5442032085673591366657251660804575198155559225316E-6L) // j26 // 21 to 26 Added for long double. }; // static const T q[] diff --git a/boost/math/special_functions/prime.hpp b/boost/math/special_functions/prime.hpp index 858d96d10b..31a9cba868 100644 --- a/boost/math/special_functions/prime.hpp +++ b/boost/math/special_functions/prime.hpp @@ -1232,7 +1232,7 @@ namespace boost{ namespace math{ return boost::math::prime(n, boost::math::policies::policy<>()); } - static const unsigned max_prime = 10000; + static const unsigned max_prime = 9999; }} // namespace boost and math diff --git a/boost/math/special_functions/sin_pi.hpp b/boost/math/special_functions/sin_pi.hpp index ae6b3e7442..de80387320 100644 --- a/boost/math/special_functions/sin_pi.hpp +++ b/boost/math/special_functions/sin_pi.hpp @@ -20,11 +20,11 @@ namespace boost{ namespace math{ namespace detail{ template -T sin_pi_imp(T x, const Policy& pol) +inline T sin_pi_imp(T x, const Policy& pol) { BOOST_MATH_STD_USING // ADL of std names if(x < 0) - return -sin_pi(-x); + return -sin_pi_imp(T(-x), pol); // sin of pi*x: bool invert; if(x < 0.5) @@ -38,7 +38,7 @@ T sin_pi_imp(T x, const Policy& pol) invert = false; T rem = floor(x); - if(itrunc(rem, pol) & 1) + if(iconvert(rem, pol) & 1) invert = !invert; rem = x - rem; if(rem > 0.5f) @@ -62,8 +62,11 @@ inline typename tools::promote_args::type sin_pi(T x, const Policy&) policies::promote_float, policies::promote_double, policies::discrete_quantile<>, - policies::assert_undefined<> >::type forwarding_policy; - return policies::checked_narrowing_cast(boost::math::detail::sin_pi_imp(x, forwarding_policy()), "cos_pi"); + policies::assert_undefined<>, + // We want to igore overflows since the result is in [-1,1] and the + // check slows the code down considerably. + policies::overflow_error >::type forwarding_policy; + return policies::checked_narrowing_cast(boost::math::detail::sin_pi_imp(x, forwarding_policy()), "sin_pi"); } template diff --git a/boost/math/special_functions/trunc.hpp b/boost/math/special_functions/trunc.hpp index 8a7b042449..2671a0dc88 100644 --- a/boost/math/special_functions/trunc.hpp +++ b/boost/math/special_functions/trunc.hpp @@ -14,6 +14,8 @@ #include #include #include +#include +#include namespace boost{ namespace math{ namespace detail{ @@ -106,6 +108,55 @@ inline boost::long_long_type lltrunc(const T& v) #endif +template +inline typename boost::enable_if_c::value, int>::type + iconvert(const T& v, const Policy&) +{ + return static_cast(v); +} + +template +inline typename boost::disable_if_c::value, int>::type + iconvert(const T& v, const Policy& pol) +{ + using boost::math::itrunc; + return itrunc(v, pol); +} + +template +inline typename boost::enable_if_c::value, long>::type + lconvert(const T& v, const Policy&) +{ + return static_cast(v); +} + +template +inline typename boost::disable_if_c::value, long>::type + lconvert(const T& v, const Policy& pol) +{ + using boost::math::ltrunc; + return ltrunc(v, pol); +} + +#ifdef BOOST_HAS_LONG_LONG + +template +inline typename boost::enable_if_c::value, boost::long_long_type>::type + llconvertert(const T& v, const Policy&) +{ + return static_cast(v); +} + +template +inline typename boost::disable_if_c::value, boost::long_long_type>::type + llconvertert(const T& v, const Policy& pol) +{ + using boost::math::lltrunc; + return lltrunc(v, pol); +} + +#endif + }} // namespaces #endif // BOOST_MATH_TRUNC_HPP diff --git a/boost/math/tools/roots.hpp b/boost/math/tools/roots.hpp index e16294dc15..552f7f7ae3 100644 --- a/boost/math/tools/roots.hpp +++ b/boost/math/tools/roots.hpp @@ -32,6 +32,7 @@ #endif #include +#include #include #include @@ -245,6 +246,11 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_ boost::uintmax_t count(max_iter); +#ifdef BOOST_MATH_INSTRUMENT + std::cout << "Newton_raphson_iterate, guess = " << guess << ", min = " << min << ", max = " << max + << ", digits = " << digits << ", max_iter = " << max_iter << std::endl; +#endif + do{ last_f0 = f0; delta2 = delta1; @@ -257,7 +263,7 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_ { // Oops zero derivative!!! #ifdef BOOST_MATH_INSTRUMENT - std::cout << "Newton iteration, zero derivative found" << std::endl; + std::cout << "Newton iteration, zero derivative found!" << std::endl; #endif detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); } @@ -266,15 +272,16 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_ delta = f0 / f1; } #ifdef BOOST_MATH_INSTRUMENT - std::cout << "Newton iteration, delta = " << delta << std::endl; + std::cout << "Newton iteration " << max_iter - count << ", delta = " << delta << std::endl; #endif if(fabs(delta * 2) > fabs(delta2)) { - // last two steps haven't converged. + // Last two steps haven't converged. T shift = (delta > 0) ? (result - min) / 2 : (result - max) / 2; if ((result != 0) && (fabs(shift) > fabs(result))) { - delta = sign(delta) * result; // protect against huge jumps! + delta = sign(delta) * fabs(result) * 0.9; // Protect against huge jumps! + //delta = sign(delta) * result; // Protect against huge jumps! Failed for negative result. https://github.com/boostorg/math/issues/216 } else delta = shift; @@ -298,7 +305,7 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_ if((result == min) || (result == max)) break; } - // update brackets: + // Update brackets: if(delta > 0) max = guess; else @@ -308,13 +315,14 @@ T newton_raphson_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_ max_iter -= count; #ifdef BOOST_MATH_INSTRUMENT - std::cout << "Newton Raphson iteration, final count = " << max_iter << std::endl; + std::cout << "Newton Raphson final iteration count = " << max_iter << std::endl; static boost::uintmax_t max_count = 0; if(max_iter > max_count) { max_count = max_iter; - std::cout << "Maximum iterations: " << max_iter << std::endl; + // std::cout << "Maximum iterations: " << max_iter << std::endl; + // Puzzled what this tells us, so commented out for now? } #endif @@ -354,12 +362,131 @@ namespace detail{ } }; + template + T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count); + + template + T bracket_root_towards_max(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count) + { + using std::fabs; + // + // Move guess towards max until we bracket the root, updating min and max as we go: + // + T guess0 = guess; + T multiplier = 2; + T f_current = f0; + if (fabs(min) < fabs(max)) + { + while (--count && ((f_current < 0) == (f0 < 0))) + { + min = guess; + guess *= multiplier; + if (guess > max) + { + guess = max; + f_current = -f_current; // There must be a change of sign! + break; + } + multiplier *= 2; + unpack_0(f(guess), f_current); + } + } + else + { + // + // If min and max are negative we have to divide to head towards max: + // + while (--count && ((f_current < 0) == (f0 < 0))) + { + min = guess; + guess /= multiplier; + if (guess > max) + { + guess = max; + f_current = -f_current; // There must be a change of sign! + break; + } + multiplier *= 2; + unpack_0(f(guess), f_current); + } + } + + if (count) + { + max = guess; + if (multiplier > 16) + return (guess0 - guess) + bracket_root_towards_min(f, guess, f_current, min, max, count); + } + return guess0 - (max + min) / 2; + } + + template + T bracket_root_towards_min(F f, T guess, const T& f0, T& min, T& max, boost::uintmax_t& count) + { + using std::fabs; + // + // Move guess towards min until we bracket the root, updating min and max as we go: + // + T guess0 = guess; + T multiplier = 2; + T f_current = f0; + + if (fabs(min) < fabs(max)) + { + while (--count && ((f_current < 0) == (f0 < 0))) + { + max = guess; + guess /= multiplier; + if (guess < min) + { + guess = min; + f_current = -f_current; // There must be a change of sign! + break; + } + multiplier *= 2; + unpack_0(f(guess), f_current); + } + } + else + { + // + // If min and max are negative we have to multiply to head towards min: + // + while (--count && ((f_current < 0) == (f0 < 0))) + { + max = guess; + guess *= multiplier; + if (guess < min) + { + guess = min; + f_current = -f_current; // There must be a change of sign! + break; + } + multiplier *= 2; + unpack_0(f(guess), f_current); + } + } + + if (count) + { + min = guess; + if (multiplier > 16) + return (guess0 - guess) + bracket_root_towards_max(f, guess, f_current, min, max, count); + } + return guess0 - (max + min) / 2; + } + template T second_order_root_finder(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { BOOST_MATH_STD_USING - T f0(0), f1, f2; +#ifdef BOOST_MATH_INSTRUMENT + std::cout << "Second order root iteration, guess = " << guess << ", min = " << min << ", max = " << max + << ", digits = " << digits << ", max_iter = " << max_iter << std::endl; +#endif + + T f0(0), f1, f2; T result = guess; T factor = ldexp(static_cast(1.0), 1 - digits); @@ -392,7 +519,7 @@ namespace detail{ { // Oops zero derivative!!! #ifdef BOOST_MATH_INSTRUMENT - std::cout << "Second order root iteration, zero derivative found" << std::endl; + std::cout << "Second order root iteration, zero derivative found!" << std::endl; #endif detail::handle_zero_derivative(f, last_f0, f0, delta, result, guess, min, max); } @@ -429,7 +556,7 @@ namespace detail{ // last two steps haven't converged. delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2; if ((result != 0) && (fabs(delta) > result)) - delta = sign(delta) * result; // protect against huge jumps! + delta = sign(delta) * fabs(result) * 0.9; // protect against huge jumps! // reset delta2 so that this branch will *not* be taken on the // next iteration: delta2 = delta * 3; @@ -459,10 +586,15 @@ namespace detail{ } else { - delta = (guess - min) / 2; - result = guess - delta; - if((result == min) || (result == max)) + if (fabs(float_distance(min, max)) < 2) + { + result = guess = (min + max) / 2; break; + } + delta = bracket_root_towards_min(f, guess, f0, min, max, count); + result = guess - delta; + guess = min; + continue; } } else if(result > max) @@ -480,10 +612,15 @@ namespace detail{ } else { - delta = (guess - max) / 2; - result = guess - delta; - if((result == min) || (result == max)) + if (fabs(float_distance(min, max)) < 2) + { + result = guess = (min + max) / 2; break; + } + delta = bracket_root_towards_max(f, guess, f0, min, max, count); + result = guess - delta; + guess = min; + continue; } } // update brackets: @@ -496,14 +633,12 @@ namespace detail{ max_iter -= count; #ifdef BOOST_MATH_INSTRUMENT - std::cout << "Second order root iteration, final count = " << max_iter << std::endl; + std::cout << "Second order root finder, final iteration count = " << max_iter << std::endl; #endif return result; } - -} - +} // T second_order_root_finder template T halley_iterate(F f, T guess, T min, T max, int digits, boost::uintmax_t& max_iter) BOOST_NOEXCEPT_IF(BOOST_MATH_IS_FLOAT(T) && noexcept(std::declval()(std::declval()))) { diff --git a/boost/math/tools/univariate_statistics.hpp b/boost/math/tools/univariate_statistics.hpp index 226fdf46d2..387ff859ee 100644 --- a/boost/math/tools/univariate_statistics.hpp +++ b/boost/math/tools/univariate_statistics.hpp @@ -30,12 +30,47 @@ auto mean(ForwardIterator first, ForwardIterator last) } return mu; } - else + else if constexpr (std::is_same_v::iterator_category, std::random_access_iterator_tag>) { - Real mu = 0; + size_t elements = std::distance(first, last); + Real mu0 = 0; + Real mu1 = 0; + Real mu2 = 0; + Real mu3 = 0; Real i = 1; - for(auto it = first; it != last; ++it) { - mu = mu + (*it - mu)/i; + auto end = last - (elements % 4); + for(auto it = first; it != end; it += 4) { + Real inv = Real(1)/i; + Real tmp0 = (*it - mu0); + Real tmp1 = (*(it+1) - mu1); + Real tmp2 = (*(it+2) - mu2); + Real tmp3 = (*(it+3) - mu3); + // please generate a vectorized fma here + mu0 += tmp0*inv; + mu1 += tmp1*inv; + mu2 += tmp2*inv; + mu3 += tmp3*inv; + i += 1; + } + Real num1 = Real(elements - (elements %4))/Real(4); + Real num2 = num1 + Real(elements % 4); + + for (auto it = end; it != last; ++it) + { + mu3 += (*it-mu3)/i; + i += 1; + } + + return (num1*(mu0+mu1+mu2) + num2*mu3)/Real(elements); + } + else + { + auto it = first; + Real mu = *it; + Real i = 2; + while(++it != last) + { + mu += (*it - mu)/i; i += 1; } return mu; -- cgit v1.2.3