diff --git a/doc/graphs/elliptic_integral_e__80_bit_long_double.svg b/doc/graphs/elliptic_integral_e__80_bit_long_double.svg index 9dfe3e317d..42d446dcc9 100644 --- a/doc/graphs/elliptic_integral_e__80_bit_long_double.svg +++ b/doc/graphs/elliptic_integral_e__80_bit_long_double.svg @@ -12,38 +12,38 @@ xmlns ="http://www.w3.org/2000/svg" - + - + - - - - - - + + + + + + -0 -1 -0 --1 +0 +1 +0 +-1 -0 -2 -0 --2 +0 +2 +0 +-2 - + -x - - +x + + - - + + Errors in Elliptic Integral E, 80 bit long double diff --git a/doc/graphs/elliptic_integral_e__double.svg b/doc/graphs/elliptic_integral_e__double.svg index 84e71db211..80fce0a802 100644 --- a/doc/graphs/elliptic_integral_e__double.svg +++ b/doc/graphs/elliptic_integral_e__double.svg @@ -12,70 +12,39 @@ xmlns ="http://www.w3.org/2000/svg" - + - + - - - - - - + + + + + + -0 +0 1 -0 --1 +0 +-1 -0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 -1.8 -2 -0 --0.2 --0.4 --0.6 --0.8 --1 --1.2 --1.4 --1.6 --1.8 --2 --2.2 --2.4 --2.6 --2.8 --3 --3.2 --3.4 --3.6 --3.8 --4 --4.2 --4.4 --4.6 --4.8 +0 +2 +0 +-2 +-4 -x - - +x + + - - + + Errors in Elliptic Integral E, double diff --git a/doc/graphs/elliptic_integral_k__80_bit_long_double.svg b/doc/graphs/elliptic_integral_k__80_bit_long_double.svg index 133d9a77ed..3c4d184e39 100644 --- a/doc/graphs/elliptic_integral_k__80_bit_long_double.svg +++ b/doc/graphs/elliptic_integral_k__80_bit_long_double.svg @@ -12,39 +12,39 @@ xmlns ="http://www.w3.org/2000/svg" - + - + - - - - - - + + + + + + -0 -1 -0 --1 +0 +1 +0 +-1 -0 -20 -0 --20 --40 +0 +20 +0 +-20 +-40 - + -x - - +x + + - - + + Errors in Elliptic Integral K, 80 bit long double diff --git a/doc/graphs/elliptic_integral_k__double.svg b/doc/graphs/elliptic_integral_k__double.svg index b7731b5808..9b9244addb 100644 --- a/doc/graphs/elliptic_integral_k__double.svg +++ b/doc/graphs/elliptic_integral_k__double.svg @@ -39,11 +39,11 @@ xmlns ="http://www.w3.org/2000/svg" x - - + + - - + + Errors in Elliptic Integral K, double diff --git a/doc/graphs/plot_1d_errors.cpp b/doc/graphs/plot_1d_errors.cpp index 6ccb813505..0217b66ee2 100644 --- a/doc/graphs/plot_1d_errors.cpp +++ b/doc/graphs/plot_1d_errors.cpp @@ -43,7 +43,7 @@ void plot_errors_1d(F f, Real start, Real end, unsigned points, const char* func Real pos = start; Real interval = (end - start) / points; - std::map points_upper, points_lower; + std::map points_upper, points_lower; Real max_distance(0), min_distance(0), max_error(0), max_error_location(0); @@ -60,19 +60,19 @@ void plot_errors_1d(F f, Real start, Real end, unsigned points, const char* func Real exact_value = static_cast(f(mp_type(pos))); Real distance = boost::math::sign(found_value - exact_value) * boost::math::epsilon_difference(found_value, exact_value); Real bin = start + ((end - start) / num_bins) * boost::math::itrunc(num_bins * (pos - start) / (end - start)); - if (points_lower.find(bin) == points_lower.end()) - points_lower[bin] = 0; - if (points_upper.find(bin) == points_upper.end()) - points_upper[bin] = 0; + if (points_lower.find((double)bin) == points_lower.end()) + points_lower[(double)bin] = 0; + if (points_upper.find((double)bin) == points_upper.end()) + points_upper[(double)bin] = 0; if (distance > 0) { - if (points_upper[bin] < distance) - points_upper[bin] = (std::min)(distance, max_y_scale); + if (points_upper[(double)bin] < distance) + points_upper[(double)bin] = (double)(std::min)(distance, max_y_scale); } else { - if (points_lower[bin] > distance) - points_lower[bin] = (std::max)(distance, -max_y_scale); + if (points_lower[(double)bin] > distance) + points_lower[(double)bin] = (double)(std::max)(distance, -max_y_scale); } if (max_distance < distance) max_distance = (std::min)(distance, max_y_scale); diff --git a/doc/sf/ellint_legendre.qbk b/doc/sf/ellint_legendre.qbk index 82a0e64d09..c780a9b019 100644 --- a/doc/sf/ellint_legendre.qbk +++ b/doc/sf/ellint_legendre.qbk @@ -97,7 +97,12 @@ NTL::RR at 1000-bit precision and this implementation. [heading Implementation] -These functions are implemented in terms of Carlson's integrals using the relations: +For up to 80-bit long double precision the complete versions of these functions +are implemented as Taylor series expansions as in: +"Fast computation of complete elliptic integrals and Jacobian elliptic functions", +Celestial Mechanics and Dynamical Astronomy, April 2012. + +Otherwise these functions are implemented in terms of Carlson's integrals using the relations: [equation ellint19] @@ -198,7 +203,12 @@ NTL::RR at 1000-bit precision and this implementation. [heading Implementation] -These functions are implemented in terms of Carlson's integrals +For up to 80-bit long double precision the complete versions of these functions +are implemented as Taylor series expansions as in: +"Fast computation of complete elliptic integrals and Jacobian elliptic functions", +Celestial Mechanics and Dynamical Astronomy, April 2012. + +Otherwise these functions are implemented in terms of Carlson's integrals using the relations: [equation ellint21] diff --git a/include/boost/math/policies/error_handling.hpp b/include/boost/math/policies/error_handling.hpp index c4120d4482..38fd949e4d 100644 --- a/include/boost/math/policies/error_handling.hpp +++ b/include/boost/math/policies/error_handling.hpp @@ -751,7 +751,7 @@ namespace detail { template -inline bool check_overflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) +BOOST_FORCEINLINE bool check_overflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) { BOOST_MATH_STD_USING if(fabs(val) > tools::max_value()) @@ -763,7 +763,7 @@ inline bool check_overflow(T val, R* result, const char* function, const Policy& return false; } template -inline bool check_overflow(std::complex val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) +BOOST_FORCEINLINE bool check_overflow(std::complex val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) { typedef typename R::value_type r_type; r_type re, im; @@ -773,7 +773,7 @@ inline bool check_overflow(std::complex val, R* result, const char* function, return r; } template -inline bool check_underflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) +BOOST_FORCEINLINE bool check_underflow(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) { if((val != 0) && (static_cast(val) == 0)) { @@ -783,7 +783,7 @@ inline bool check_underflow(T val, R* result, const char* function, const Policy return false; } template -inline bool check_underflow(std::complex val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) +BOOST_FORCEINLINE bool check_underflow(std::complex val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) { typedef typename R::value_type r_type; r_type re, im; @@ -793,7 +793,7 @@ inline bool check_underflow(std::complex val, R* result, const char* function return r; } template -inline bool check_denorm(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) +BOOST_FORCEINLINE bool check_denorm(T val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) { BOOST_MATH_STD_USING if((fabs(val) < static_cast(tools::min_value())) && (static_cast(val) != 0)) @@ -804,7 +804,7 @@ inline bool check_denorm(T val, R* result, const char* function, const Policy& p return false; } template -inline bool check_denorm(std::complex val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) +BOOST_FORCEINLINE bool check_denorm(std::complex val, R* result, const char* function, const Policy& pol) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && (Policy::value != throw_on_error) && (Policy::value != user_error)) { typedef typename R::value_type r_type; r_type re, im; @@ -816,28 +816,28 @@ inline bool check_denorm(std::complex val, R* result, const char* function, c // Default instantiations with ignore_error policy. template -inline constexpr bool check_overflow(T /* val */, R* /* result */, const char* /* function */, const overflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) +BOOST_FORCEINLINE constexpr bool check_overflow(T /* val */, R* /* result */, const char* /* function */, const overflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) { return false; } template -inline constexpr bool check_overflow(std::complex /* val */, R* /* result */, const char* /* function */, const overflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) +BOOST_FORCEINLINE constexpr bool check_overflow(std::complex /* val */, R* /* result */, const char* /* function */, const overflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) { return false; } template -inline constexpr bool check_underflow(T /* val */, R* /* result */, const char* /* function */, const underflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) +BOOST_FORCEINLINE constexpr bool check_underflow(T /* val */, R* /* result */, const char* /* function */, const underflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) { return false; } template -inline constexpr bool check_underflow(std::complex /* val */, R* /* result */, const char* /* function */, const underflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) +BOOST_FORCEINLINE constexpr bool check_underflow(std::complex /* val */, R* /* result */, const char* /* function */, const underflow_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) { return false; } template -inline constexpr bool check_denorm(T /* val */, R* /* result*/, const char* /* function */, const denorm_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) +BOOST_FORCEINLINE constexpr bool check_denorm(T /* val */, R* /* result*/, const char* /* function */, const denorm_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) { return false; } template -inline constexpr bool check_denorm(std::complex /* val */, R* /* result*/, const char* /* function */, const denorm_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) +BOOST_FORCEINLINE constexpr bool check_denorm(std::complex /* val */, R* /* result*/, const char* /* function */, const denorm_error&) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T)) { return false; } } // namespace detail template -inline R checked_narrowing_cast(T val, const char* function) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && is_noexcept_error_policy::value) +BOOST_FORCEINLINE R checked_narrowing_cast(T val, const char* function) noexcept(BOOST_MATH_IS_FLOAT(R) && BOOST_MATH_IS_FLOAT(T) && is_noexcept_error_policy::value) { typedef typename Policy::overflow_error_type overflow_type; typedef typename Policy::underflow_error_type underflow_type; diff --git a/include/boost/math/special_functions/ellint_1.hpp b/include/boost/math/special_functions/ellint_1.hpp index cc634ad4ef..a47ad76766 100644 --- a/include/boost/math/special_functions/ellint_1.hpp +++ b/include/boost/math/special_functions/ellint_1.hpp @@ -36,7 +36,11 @@ typename tools::promote_args::type ellint_1(T1 k, T2 phi, const Policy& namespace detail{ template -T ellint_k_imp(T k, const Policy& pol); +T ellint_k_imp(T k, const Policy& pol, std::integral_constant const&); +template +T ellint_k_imp(T k, const Policy& pol, std::integral_constant const&); +template +T ellint_k_imp(T k, const Policy& pol, std::integral_constant const&); // Elliptic integral (Legendre form) of the first kind template @@ -71,7 +75,12 @@ T ellint_f_imp(T phi, T k, const Policy& pol) { // Phi is so large that phi%pi is necessarily zero (or garbage), // just return the second part of the duplication formula: - result = 2 * phi * ellint_k_imp(k, pol) / constants::pi(); + typedef std::integral_constant::value&& std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 + > precision_tag_type; + + result = 2 * phi * ellint_k_imp(k, pol, precision_tag_type()) / constants::pi(); BOOST_MATH_INSTRUMENT_VARIABLE(result); } else @@ -122,7 +131,12 @@ T ellint_f_imp(T phi, T k, const Policy& pol) BOOST_MATH_INSTRUMENT_VARIABLE(result); if(m != 0) { - result += m * ellint_k_imp(k, pol); + typedef std::integral_constant::value&& std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 + > precision_tag_type; + + result += m * ellint_k_imp(k, pol, precision_tag_type()); BOOST_MATH_INSTRUMENT_VARIABLE(result); } } @@ -131,7 +145,7 @@ T ellint_f_imp(T phi, T k, const Policy& pol) // Complete elliptic integral (Legendre form) of the first kind template -T ellint_k_imp(T k, const Policy& pol) +T ellint_k_imp(T k, const Policy& pol, std::integral_constant const&) { BOOST_MATH_STD_USING using namespace boost::math::tools; @@ -156,16 +170,599 @@ T ellint_k_imp(T k, const Policy& pol) return value; } +// +// Special versions for double and 80-bit long double precision, +// double precision versions use the coefficients from: +// "Fast computation of complete elliptic integrals and Jacobian elliptic functions", +// Celestial Mechanics and Dynamical Astronomy, April 2012. +// +// Higher precision coefficients for 80-bit long doubles can be calculated +// using for example: +// Table[N[SeriesCoefficient[ EllipticK [ m ], { m, 875/1000, i} ], 20], {i, 0, 24}] +// and checking the value of the first neglected term with: +// N[SeriesCoefficient[ EllipticK [ m ], { m, 875/1000, 24} ], 20] * (2.5/100)^24 +// +// For m > 0.9 we don't use the method of the paper above, but simply call our +// existing routines. The routine used in the above paper was tried (and is +// archived in the code below), but was found to have slightly higher error rates. +// +template +BOOST_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, std::integral_constant const&) +{ + using std::abs; + using namespace boost::math::tools; + + T m = k * k; + + switch (static_cast(m * 20)) + { + case 0: + case 1: + //if (m < 0.1) + { + constexpr T coef[] = + { + 1.591003453790792180, + 0.416000743991786912, + 0.245791514264103415, + 0.179481482914906162, + 0.144556057087555150, + 0.123200993312427711, + 0.108938811574293531, + 0.098853409871592910, + 0.091439629201749751, + 0.085842591595413900, + 0.081541118718303215, + 0.078199656811256481910 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.05); + } + case 2: + case 3: + //else if (m < 0.2) + { + constexpr T coef[] = + { + 1.635256732264579992, + 0.471190626148732291, + 0.309728410831499587, + 0.252208311773135699, + 0.226725623219684650, + 0.215774446729585976, + 0.213108771877348910, + 0.216029124605188282, + 0.223255831633057896, + 0.234180501294209925, + 0.248557682972264071, + 0.266363809892617521 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.15); + } + case 4: + case 5: + //else if (m < 0.3) + { + constexpr T coef[] = + { + 1.685750354812596043, + 0.541731848613280329, + 0.401524438390690257, + 0.369642473420889090, + 0.376060715354583645, + 0.405235887085125919, + 0.453294381753999079, + 0.520518947651184205, + 0.609426039204995055, + 0.724263522282908870, + 0.871013847709812357, + 1.057652872753547036 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.25); + } + case 6: + case 7: + //else if (m < 0.4) + { + constexpr T coef[] = + { + 1.744350597225613243, + 0.634864275371935304, + 0.539842564164445538, + 0.571892705193787391, + 0.670295136265406100, + 0.832586590010977199, + 1.073857448247933265, + 1.422091460675497751, + 1.920387183402304829, + 2.632552548331654201, + 3.652109747319039160, + 5.115867135558865806, + 7.224080007363877411 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.35); + } + case 8: + case 9: + //else if (m < 0.5) + { + constexpr T coef[] = + { + 1.813883936816982644, + 0.763163245700557246, + 0.761928605321595831, + 0.951074653668427927, + 1.315180671703161215, + 1.928560693477410941, + 2.937509342531378755, + 4.594894405442878062, + 7.330071221881720772, + 11.87151259742530180, + 19.45851374822937738, + 32.20638657246426863, + 53.73749198700554656, + 90.27388602940998849 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.45); + } + case 10: + case 11: + //else if (m < 0.6) + { + constexpr T coef[] = + { + 1.898924910271553526, + 0.950521794618244435, + 1.151077589959015808, + 1.750239106986300540, + 2.952676812636875180, + 5.285800396121450889, + 9.832485716659979747, + 18.78714868327559562, + 36.61468615273698145, + 72.45292395127771801, + 145.1079577347069102, + 293.4786396308497026, + 598.3851815055010179, + 1228.420013075863451, + 2536.529755382764488 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.55); + } + case 12: + case 13: + //else if (m < 0.7) + { + constexpr T coef[] = + { + 2.007598398424376302, + 1.248457231212347337, + 1.926234657076479729, + 3.751289640087587680, + 8.119944554932045802, + 18.66572130873555361, + 44.60392484291437063, + 109.5092054309498377, + 274.2779548232413480, + 697.5598008606326163, + 1795.716014500247129, + 4668.381716790389910, + 12235.76246813664335, + 32290.17809718320818, + 85713.07608195964685, + 228672.1890493117096, + 612757.2711915852774 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.65); + } + case 14: + case 15: + //else if (m < 0.8) + { + constexpr T coef[] = + { + 2.156515647499643235, + 1.791805641849463243, + 3.826751287465713147, + 10.38672468363797208, + 31.40331405468070290, + 100.9237039498695416, + 337.3268282632272897, + 1158.707930567827917, + 4060.990742193632092, + 14454.00184034344795, + 52076.66107599404803, + 189493.6591462156887, + 695184.5762413896145, + 2567994.048255284686, + 9541921.966748386322, + 35634927.44218076174, + 133669298.4612040871, + 503352186.6866284541, + 1901975729.538660119, + 7208915015.330103756 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.75); + } + case 16: + //else if (m < 0.85) + { + constexpr T coef[] = + { + 2.318122621712510589, + 2.616920150291232841, + 7.897935075731355823, + 30.50239715446672327, + 131.4869365523528456, + 602.9847637356491617, + 2877.024617809972641, + 14110.51991915180325, + 70621.44088156540229, + 358977.2665825309926, + 1847238.263723971684, + 9600515.416049214109, + 50307677.08502366879, + 265444188.6527127967, + 1408862325.028702687, + 7515687935.373774627 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.825); + } + case 17: + //else if (m < 0.90) + { + constexpr T coef[] = + { + 2.473596173751343912, + 3.727624244118099310, + 15.60739303554930496, + 84.12850842805887747, + 506.9818197040613935, + 3252.277058145123644, + 21713.24241957434256, + 149037.0451890932766, + 1043999.331089990839, + 7427974.817042038995, + 53503839.67558661151, + 389249886.9948708474, + 2855288351.100810619, + 21090077038.76684053, + 156699833947.7902014, + 1170222242422.439893, + 8777948323668.937971, + 66101242752484.95041, + 499488053713388.7989, + 37859743397240299.20 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.875); + } + default: + // + // This handles all cases where m > 0.9, + // including all error handling: + // + return ellint_k_imp(k, pol, std::integral_constant()); +#if 0 + else + { + T lambda_prime = (1 - sqrt(k)) / (2 * (1 + sqrt(k))); + T k_prime = ellint_k(sqrt((1 - k) * (1 + k))); // K(m') + T lambda_prime_4th = boost::math::pow<4>(lambda_prime); + T q_prime = ((((((20910 * lambda_prime_4th) + 1707) * lambda_prime_4th + 150) * lambda_prime_4th + 15) * lambda_prime_4th + 2) * lambda_prime_4th + 1) * lambda_prime; + /*T q_prime_2 = lambda_prime + + 2 * boost::math::pow<5>(lambda_prime) + + 15 * boost::math::pow<9>(lambda_prime) + + 150 * boost::math::pow<13>(lambda_prime) + + 1707 * boost::math::pow<17>(lambda_prime) + + 20910 * boost::math::pow<21>(lambda_prime);*/ + return -log(q_prime) * k_prime / boost::math::constants::pi(); + } +#endif + } +} +template +BOOST_FORCEINLINE T ellint_k_imp(T k, const Policy& pol, std::integral_constant const&) +{ + using std::abs; + using namespace boost::math::tools; + + T m = k * k; + switch (static_cast(m * 20)) + { + case 0: + case 1: + { + constexpr T coef[] = + { + 1.5910034537907921801L, + 0.41600074399178691174L, + 0.24579151426410341536L, + 0.17948148291490616181L, + 0.14455605708755514976L, + 0.12320099331242771115L, + 0.10893881157429353105L, + 0.098853409871592910399L, + 0.091439629201749751268L, + 0.085842591595413899672L, + 0.081541118718303214749L, + 0.078199656811256481910L, + 0.075592617535422415648L, + 0.073562939365441925050L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.05L); + } + case 2: + case 3: + { + constexpr T coef[] = + { + 1.6352567322645799924L, + 0.47119062614873229055L, + 0.30972841083149958708L, + 0.25220831177313569923L, + 0.22672562321968464974L, + 0.21577444672958597588L, + 0.21310877187734890963L, + 0.21602912460518828154L, + 0.22325583163305789567L, + 0.23418050129420992492L, + 0.24855768297226407136L, + 0.26636380989261752077L, + 0.28772845215611466775L, + 0.31290024539780334906L, + 0.34223105446381299902L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.15L); + } + case 4: + case 5: + { + constexpr T coef[] = + { + 1.6857503548125960429L, + 0.54173184861328032882L, + 0.40152443839069025682L, + 0.36964247342088908995L, + 0.37606071535458364462L, + 0.40523588708512591863L, + 0.45329438175399907924L, + 0.52051894765118420473L, + 0.60942603920499505544L, + 0.72426352228290886975L, + 0.87101384770981235737L, + 1.0576528727535470365L, + 1.2945970872087764321L, + 1.5953368253888783747L, + 1.9772844873556364793L, + 2.4628890581910021287L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.25L); + } + case 6: + case 7: + { + constexpr T coef[] = + { + 1.7443505972256132429L, + 0.63486427537193530383L, + 0.53984256416444553751L, + 0.57189270519378739093L, + 0.67029513626540610034L, + 0.83258659001097719939L, + 1.0738574482479332654L, + 1.4220914606754977514L, + 1.9203871834023048288L, + 2.6325525483316542006L, + 3.6521097473190391602L, + 5.1158671355588658061L, + 7.2240800073638774108L, + 10.270306349944787227L, + 14.685616935355757348L, + 21.104114212004582734L, + 30.460132808575799413L, + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.35L); + } + case 8: + case 9: + { + constexpr T coef[] = + { + 1.8138839368169826437L, + 0.76316324570055724607L, + 0.76192860532159583095L, + 0.95107465366842792679L, + 1.3151806717031612153L, + 1.9285606934774109412L, + 2.9375093425313787550L, + 4.5948944054428780618L, + 7.3300712218817207718L, + 11.871512597425301798L, + 19.458513748229377383L, + 32.206386572464268628L, + 53.737491987005546559L, + 90.273886029409988491L, + 152.53312130253275268L, + 259.02388747148299086L, + 441.78537518096201946L, + 756.39903981567380952L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.45L); + } + case 10: + case 11: + { + constexpr T coef[] = + { + 1.8989249102715535257L, + 0.95052179461824443490L, + 1.1510775899590158079L, + 1.7502391069863005399L, + 2.9526768126368751802L, + 5.2858003961214508892L, + 9.8324857166599797471L, + 18.787148683275595622L, + 36.614686152736981447L, + 72.452923951277718013L, + 145.10795773470691023L, + 293.47863963084970259L, + 598.38518150550101790L, + 1228.4200130758634505L, + 2536.5297553827644880L, + 5263.9832725075189576L, + 10972.138126273491753L, + 22958.388550988306870L, + 48203.103373625406989L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.55L); + } + case 12: + case 13: + { + constexpr T coef[] = + { + 2.0075983984243763017L, + 1.2484572312123473371L, + 1.9262346570764797287L, + 3.7512896400875876798L, + 8.1199445549320458022L, + 18.665721308735553611L, + 44.603924842914370633L, + 109.50920543094983774L, + 274.27795482324134804L, + 697.55980086063261629L, + 1795.7160145002471293L, + 4668.3817167903899100L, + 12235.762468136643348L, + 32290.178097183208178L, + 85713.076081959646847L, + 228672.18904931170958L, + 612757.27119158527740L, + 1.6483233976504668314e6L, + 4.4492251046211960936e6L, + 1.2046317340783185238e7L, + 3.2705187507963254185e7L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.65L); + } + case 14: + case 15: + { + constexpr T coef[] = + { + 2.1565156474996432354L, + 1.7918056418494632425L, + 3.8267512874657131470L, + 10.386724683637972080L, + 31.403314054680702901L, + 100.92370394986954165L, + 337.32682826322728966L, + 1158.7079305678279173L, + 4060.9907421936320917L, + 14454.001840343447947L, + 52076.661075994048028L, + 189493.65914621568866L, + 695184.57624138961450L, + 2.5679940482552846861e6L, + 9.5419219667483863221e6L, + 3.5634927442180761743e7L, + 1.3366929846120408712e8L, + 5.0335218668662845411e8L, + 1.9019757295386601192e9L, + 7.2089150153301037563e9L, + 2.7398741806339510931e10L, + 1.0439286724885300495e11L, + 3.9864875581513728207e11L, + 1.5254661585564745591e12L, + 5.8483259088850315936e12 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.75L); + } + case 16: + { + constexpr T coef[] = + { + 2.3181226217125105894L, + 2.6169201502912328409L, + 7.8979350757313558232L, + 30.502397154466723270L, + 131.48693655235284561L, + 602.98476373564916170L, + 2877.0246178099726410L, + 14110.519919151803247L, + 70621.440881565402289L, + 358977.26658253099258L, + 1.8472382637239716844e6L, + 9.6005154160492141090e6L, + 5.0307677085023668786e7L, + 2.6544418865271279673e8L, + 1.4088623250287026866e9L, + 7.5156879353737746270e9L, + 4.0270783964955246149e10L, + 2.1662089325801126339e11L, + 1.1692489201929996116e12L, + 6.3306543358985679881e12 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.825L); + } + case 17: + { + constexpr T coef[] = + { + 2.4735961737513439120L, + 3.7276242441180993105L, + 15.607393035549304964L, + 84.128508428058877470L, + 506.98181970406139349L, + 3252.2770581451236438L, + 21713.242419574342564L, + 149037.04518909327662L, + 1.0439993310899908390e6L, + 7.4279748170420389947e6L, + 5.3503839675586611510e7L, + 3.8924988699487084738e8L, + 2.8552883511008106195e9L, + 2.1090077038766840525e10L, + 1.5669983394779020136e11L, + 1.1702222424224398927e12L, + 8.7779483236689379709e12L, + 6.6101242752484950408e13L, + 4.9948805371338879891e14L, + 3.7859743397240299201e15L, + 2.8775996123036112296e16L, + 2.1926346839925760143e17L, + 1.6744985438468349361e18L, + 1.2814410112866546052e19L, + 9.8249807041031260167e19 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.875L); + } + default: + // + // All cases where m > 0.9 + // including all error handling: + // + return ellint_k_imp(k, pol, std::integral_constant()); + } +} + template -inline typename tools::promote_args::type ellint_1(T k, const Policy& pol, const std::true_type&) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_1(T k, const Policy& pol, const std::true_type&) { typedef typename tools::promote_args::type result_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::ellint_k_imp(static_cast(k), pol), "boost::math::ellint_1<%1%>(%1%)"); + typedef std::integral_constant::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 +#endif + > precision_tag_type; + return policies::checked_narrowing_cast(detail::ellint_k_imp(static_cast(k), pol, precision_tag_type()), "boost::math::ellint_1<%1%>(%1%)"); } template -inline typename tools::promote_args::type ellint_1(T1 k, T2 phi, const std::false_type&) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_1(T1 k, T2 phi, const std::false_type&) { return boost::math::ellint_1(k, phi, policies::policy<>()); } @@ -174,14 +771,14 @@ inline typename tools::promote_args::type ellint_1(T1 k, T2 phi, const s // Complete elliptic integral (Legendre form) of the first kind template -inline typename tools::promote_args::type ellint_1(T k) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_1(T k) { return ellint_1(k, policies::policy<>()); } // Elliptic integral (Legendre form) of the first kind template -inline typename tools::promote_args::type ellint_1(T1 k, T2 phi, const Policy& pol) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_1(T1 k, T2 phi, const Policy& pol) { typedef typename tools::promote_args::type result_type; typedef typename policies::evaluation::type value_type; @@ -189,7 +786,7 @@ inline typename tools::promote_args::type ellint_1(T1 k, T2 phi, const P } template -inline typename tools::promote_args::type ellint_1(T1 k, T2 phi) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_1(T1 k, T2 phi) { typedef typename policies::is_policy::type tag_type; return detail::ellint_1(k, phi, tag_type()); diff --git a/include/boost/math/special_functions/ellint_2.hpp b/include/boost/math/special_functions/ellint_2.hpp index d5cd094d97..e135264ccf 100644 --- a/include/boost/math/special_functions/ellint_2.hpp +++ b/include/boost/math/special_functions/ellint_2.hpp @@ -38,7 +38,11 @@ typename tools::promote_args::type ellint_2(T1 k, T2 phi, const Policy& namespace detail{ template -T ellint_e_imp(T k, const Policy& pol); +T ellint_e_imp(T k, const Policy& pol, const std::integral_constant&); +template +T ellint_e_imp(T k, const Policy& pol, const std::integral_constant&); +template +T ellint_e_imp(T k, const Policy& pol, const std::integral_constant&); // Elliptic integral (Legendre form) of the second kind template @@ -67,9 +71,13 @@ T ellint_e_imp(T phi, T k, const Policy& pol) } else if(phi > 1 / tools::epsilon()) { + typedef std::integral_constant::value&& std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 + > precision_tag_type; // Phi is so large that phi%pi is necessarily zero (or garbage), // just return the second part of the duplication formula: - result = 2 * phi * ellint_e_imp(k, pol) / constants::pi(); + result = 2 * phi * ellint_e_imp(k, pol, precision_tag_type()) / constants::pi(); } else if(k == 0) { @@ -128,15 +136,21 @@ T ellint_e_imp(T phi, T k, const Policy& pol) T cm1 = cosp * cosp / (sinp * sinp); // c - 1 result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); } - if(m != 0) - result += m * ellint_e_imp(k, pol); + if (m != 0) + { + typedef std::integral_constant::value&& std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 + > precision_tag_type; + result += m * ellint_e_imp(k, pol, precision_tag_type()); + } } return invert ? T(-result) : result; } // Complete elliptic integral (Legendre form) of the second kind template -T ellint_e_imp(T k, const Policy& pol) +T ellint_e_imp(T k, const Policy& pol, std::integral_constant const&) { BOOST_MATH_STD_USING using namespace boost::math::tools; @@ -159,18 +173,549 @@ T ellint_e_imp(T k, const Policy& pol) return value; } +// +// Special versions for double and 80-bit long double precision, +// double precision versions use the coefficients from: +// "Fast computation of complete elliptic integrals and Jacobian elliptic functions", +// Celestial Mechanics and Dynamical Astronomy, April 2012. +// +// Higher precision coefficients for 80-bit long doubles can be calculated +// using for example: +// Table[N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, i} ], 20], {i, 0, 24}] +// and checking the value of the first neglected term with: +// N[SeriesCoefficient[ EllipticE [ m ], { m, 875/1000, 24} ], 20] * (2.5/100)^24 +// +// For m > 0.9 we don't use the method of the paper above, but simply call our +// existing routines. +// +template +BOOST_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant const&) +{ + using std::abs; + using namespace boost::math::tools; + + T m = k * k; + switch (static_cast(20 * m)) + { + case 0: + case 1: + //if (m < 0.1) + { + constexpr T coef[] = + { + 1.550973351780472328, + -0.400301020103198524, + -0.078498619442941939, + -0.034318853117591992, + -0.019718043317365499, + -0.013059507731993309, + -0.009442372874146547, + -0.007246728512402157, + -0.005807424012956090, + -0.004809187786009338, + -0.004086399233255150 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.05); + } + case 2: + case 3: + //else if (m < 0.2) + { + constexpr T coef[] = + { + 1.510121832092819728, + -0.417116333905867549, + -0.090123820404774569, + -0.043729944019084312, + -0.027965493064761785, + -0.020644781177568105, + -0.016650786739707238, + -0.014261960828842520, + -0.012759847429264803, + -0.011799303775587354, + -0.011197445703074968 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.15); + } + case 4: + case 5: + //else if (m < 0.3) + { + constexpr T coef[] = + { + 1.467462209339427155, + -0.436576290946337775, + -0.105155557666942554, + -0.057371843593241730, + -0.041391627727340220, + -0.034527728505280841, + -0.031495443512532783, + -0.030527000890325277, + -0.030916984019238900, + -0.032371395314758122, + -0.034789960386404158 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.25); + } + case 6: + case 7: + //else if (m < 0.4) + { + constexpr T coef[] = + { + 1.422691133490879171, + -0.459513519621048674, + -0.125250539822061878, + -0.078138545094409477, + -0.064714278472050002, + -0.062084339131730311, + -0.065197032815572477, + -0.072793895362578779, + -0.084959075171781003, + -0.102539850131045997, + -0.127053585157696036, + -0.160791120691274606 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.35); + } + case 8: + case 9: + //else if (m < 0.5) + { + constexpr T coef[] = + { + 1.375401971871116291, + -0.487202183273184837, + -0.153311701348540228, + -0.111849444917027833, + -0.108840952523135768, + -0.122954223120269076, + -0.152217163962035047, + -0.200495323642697339, + -0.276174333067751758, + -0.393513114304375851, + -0.575754406027879147, + -0.860523235727239756, + -1.308833205758540162 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.45); + } + case 10: + case 11: + //else if (m < 0.6) + { + constexpr T coef[] = + { + 1.325024497958230082, + -0.521727647557566767, + -0.194906430482126213, + -0.171623726822011264, + -0.202754652926419141, + -0.278798953118534762, + -0.420698457281005762, + -0.675948400853106021, + -1.136343121839229244, + -1.976721143954398261, + -3.531696773095722506, + -6.446753640156048150, + -11.97703130208884026 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.55); + } + case 12: + case 13: + //else if (m < 0.7) + { + constexpr T coef[] = + { + 1.270707479650149744, + -0.566839168287866583, + -0.262160793432492598, + -0.292244173533077419, + -0.440397840850423189, + -0.774947641381397458, + -1.498870837987561088, + -3.089708310445186667, + -6.667595903381001064, + -14.89436036517319078, + -34.18120574251449024, + -80.15895841905397306, + -191.3489480762984920, + -463.5938853480342030, + -1137.380822169360061 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.65); + } + case 14: + case 15: + //else if (m < 0.8) + { + constexpr T coef[] = + { + 1.211056027568459525, + -0.630306413287455807, + -0.387166409520669145, + -0.592278235311934603, + -1.237555584513049844, + -3.032056661745247199, + -8.181688221573590762, + -23.55507217389693250, + -71.04099935893064956, + -221.8796853192349888, + -712.1364793277635425, + -2336.125331440396407, + -7801.945954775964673, + -26448.19586059191933, + -90799.48341621365251, + -315126.0406449163424, + -1104011.344311591159 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.75); + } + case 16: + //else if (m < 0.85) + { + constexpr T coef[] = + { + 1.161307152196282836, + -0.701100284555289548, + -0.580551474465437362, + -1.243693061077786614, + -3.679383613496634879, + -12.81590924337895775, + -49.25672530759985272, + -202.1818735434090269, + -869.8602699308701437, + -3877.005847313289571, + -17761.70710170939814, + -83182.69029154232061, + -396650.4505013548170, + -1920033.413682634405 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.825); + } + case 17: + //else if (m < 0.90) + { + constexpr T coef[] = + { + 1.124617325119752213, + -0.770845056360909542, + -0.844794053644911362, + -2.490097309450394453, + -10.23971741154384360, + -49.74900546551479866, + -267.0986675195705196, + -1532.665883825229947, + -9222.313478526091951, + -57502.51612140314030, + -368596.1167416106063, + -2415611.088701091428, + -16120097.81581656797, + -109209938.5203089915, + -749380758.1942496220, + -5198725846.725541393, + -36409256888.12139973 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.875); + } + default: + // + // All cases where m > 0.9 + // including all error handling: + // + return ellint_e_imp(k, pol, std::integral_constant()); + } +} +template +BOOST_FORCEINLINE T ellint_e_imp(T k, const Policy& pol, std::integral_constant const&) +{ + using std::abs; + using namespace boost::math::tools; + + T m = k * k; + switch (static_cast(20 * m)) + { + case 0: + case 1: + //if (m < 0.1) + { + constexpr T coef[] = + { + 1.5509733517804723277L, + -0.40030102010319852390L, + -0.078498619442941939212L, + -0.034318853117591992417L, + -0.019718043317365499309L, + -0.013059507731993309191L, + -0.0094423728741465473894L, + -0.0072467285124021568126L, + -0.0058074240129560897940L, + -0.0048091877860093381762L, + -0.0040863992332551506768L, + -0.0035450302604139562644L, + -0.0031283511188028336315L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.05L); + } + case 2: + case 3: + //else if (m < 0.2) + { + constexpr T coef[] = + { + 1.5101218320928197276L, + -0.41711633390586754922L, + -0.090123820404774568894L, + -0.043729944019084311555L, + -0.027965493064761784548L, + -0.020644781177568105268L, + -0.016650786739707238037L, + -0.014261960828842519634L, + -0.012759847429264802627L, + -0.011799303775587354169L, + -0.011197445703074968018L, + -0.010850368064799902735L, + -0.010696133481060989818L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.15L); + } + case 4: + case 5: + //else if (m < 0.3L) + { + constexpr T coef[] = + { + 1.4674622093394271555L, + -0.43657629094633777482L, + -0.10515555766694255399L, + -0.057371843593241729895L, + -0.041391627727340220236L, + -0.034527728505280841188L, + -0.031495443512532782647L, + -0.030527000890325277179L, + -0.030916984019238900349L, + -0.032371395314758122268L, + -0.034789960386404158240L, + -0.038182654612387881967L, + -0.042636187648900252525L, + -0.048302272505241634467 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.25L); + } + case 6: + case 7: + //else if (m < 0.4L) + { + constexpr T coef[] = + { + 1.4226911334908791711L, + -0.45951351962104867394L, + -0.12525053982206187849L, + -0.078138545094409477156L, + -0.064714278472050001838L, + -0.062084339131730310707L, + -0.065197032815572476910L, + -0.072793895362578779473L, + -0.084959075171781003264L, + -0.10253985013104599679L, + -0.12705358515769603644L, + -0.16079112069127460621L, + -0.20705400012405941376L, + -0.27053164884730888948L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.35L); + } + case 8: + case 9: + //else if (m < 0.5L) + { + constexpr T coef[] = + { + 1.3754019718711162908L, + -0.48720218327318483652L, + -0.15331170134854022753L, + -0.11184944491702783273L, + -0.10884095252313576755L, + -0.12295422312026907610L, + -0.15221716396203504746L, + -0.20049532364269733857L, + -0.27617433306775175837L, + -0.39351311430437585139L, + -0.57575440602787914711L, + -0.86052323572723975634L, + -1.3088332057585401616L, + -2.0200280559452241745L, + -3.1566019548237606451L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.45L); + } + case 10: + case 11: + //else if (m < 0.6L) + { + constexpr T coef[] = + { + 1.3250244979582300818L, + -0.52172764755756676713L, + -0.19490643048212621262L, + -0.17162372682201126365L, + -0.20275465292641914128L, + -0.27879895311853476205L, + -0.42069845728100576224L, + -0.67594840085310602110L, + -1.1363431218392292440L, + -1.9767211439543982613L, + -3.5316967730957225064L, + -6.4467536401560481499L, + -11.977031302088840261L, + -22.581360948073964469L, + -43.109479829481450573L, + -83.186290908288807424L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.55L); + } + case 12: + case 13: + //else if (m < 0.7L) + { + constexpr T coef[] = + { + 1.2707074796501497440L, + -0.56683916828786658286L, + -0.26216079343249259779L, + -0.29224417353307741931L, + -0.44039784085042318909L, + -0.77494764138139745824L, + -1.4988708379875610880L, + -3.0897083104451866665L, + -6.6675959033810010645L, + -14.894360365173190775L, + -34.181205742514490240L, + -80.158958419053973056L, + -191.34894807629849204L, + -463.59388534803420301L, + -1137.3808221693600606L, + -2820.7073786352269339L, + -7061.1382244658715621L, + -17821.809331816437058L, + -45307.849987201897801L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.65L); + } + case 14: + case 15: + //else if (m < 0.8L) + { + constexpr T coef[] = + { + 1.2110560275684595248L, + -0.63030641328745580709L, + -0.38716640952066914514L, + -0.59227823531193460257L, + -1.2375555845130498445L, + -3.0320566617452471986L, + -8.1816882215735907624L, + -23.555072173896932503L, + -71.040999358930649565L, + -221.87968531923498875L, + -712.13647932776354253L, + -2336.1253314403964072L, + -7801.9459547759646726L, + -26448.195860591919335L, + -90799.483416213652512L, + -315126.04064491634241L, + -1.1040113443115911589e6L, + -3.8998018348056769095e6L, + -1.3876249116223745041e7L, + -4.9694982823537861149e7L, + -1.7900668836197342979e8L, + -6.4817399873722371964e8L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.75L); + } + case 16: + //else if (m < 0.85L) + { + constexpr T coef[] = + { + 1.1613071521962828360L, + -0.70110028455528954752L, + -0.58055147446543736163L, + -1.2436930610777866138L, + -3.6793836134966348789L, + -12.815909243378957753L, + -49.256725307599852720L, + -202.18187354340902693L, + -869.86026993087014372L, + -3877.0058473132895713L, + -17761.707101709398174L, + -83182.690291542320614L, + -396650.45050135481698L, + -1.9200334136826344054e6L, + -9.4131321779500838352e6L, + -4.6654858837335370627e7L, + -2.3343549352617609390e8L, + -1.1776928223045913454e9L, + -5.9850851892915740401e9L, + -3.0614702984618644983e10L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.825L); + } + case 17: + //else if (m < 0.90L) + { + constexpr T coef[] = + { + 1.1246173251197522132L, + -0.77084505636090954218L, + -0.84479405364491136236L, + -2.4900973094503944527L, + -10.239717411543843601L, + -49.749005465514798660L, + -267.09866751957051961L, + -1532.6658838252299468L, + -9222.3134785260919507L, + -57502.516121403140303L, + -368596.11674161060626L, + -2.4156110887010914281e6L, + -1.6120097815816567971e7L, + -1.0920993852030899148e8L, + -7.4938075819424962198e8L, + -5.1987258467255413931e9L, + -3.6409256888121399726e10L, + -2.5711802891217393544e11L, + -1.8290904062978796996e12L, + -1.3096838781743248404e13L, + -9.4325465851415135118e13L, + -6.8291980829471896669e14L + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.875L); + } + default: + // + // All cases where m > 0.9 + // including all error handling: + // + return ellint_e_imp(k, pol, std::integral_constant()); + } +} template -inline typename tools::promote_args::type ellint_2(T k, const Policy& pol, const std::true_type&) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_2(T k, const Policy& pol, const std::true_type&) { typedef typename tools::promote_args::type result_type; typedef typename policies::evaluation::type value_type; - return policies::checked_narrowing_cast(detail::ellint_e_imp(static_cast(k), pol), "boost::math::ellint_2<%1%>(%1%)"); + typedef std::integral_constant::value&& std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 + > precision_tag_type; + return policies::checked_narrowing_cast(detail::ellint_e_imp(static_cast(k), pol, precision_tag_type()), "boost::math::ellint_2<%1%>(%1%)"); } // Elliptic integral (Legendre form) of the second kind template -inline typename tools::promote_args::type ellint_2(T1 k, T2 phi, const std::false_type&) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_2(T1 k, T2 phi, const std::false_type&) { return boost::math::ellint_2(k, phi, policies::policy<>()); } @@ -179,21 +724,21 @@ inline typename tools::promote_args::type ellint_2(T1 k, T2 phi, const s // Complete elliptic integral (Legendre form) of the second kind template -inline typename tools::promote_args::type ellint_2(T k) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_2(T k) { return ellint_2(k, policies::policy<>()); } // Elliptic integral (Legendre form) of the second kind template -inline typename tools::promote_args::type ellint_2(T1 k, T2 phi) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_2(T1 k, T2 phi) { typedef typename policies::is_policy::type tag_type; return detail::ellint_2(k, phi, tag_type()); } template -inline typename tools::promote_args::type ellint_2(T1 k, T2 phi, const Policy& pol) +BOOST_FORCEINLINE typename tools::promote_args::type ellint_2(T1 k, T2 phi, const Policy& pol) { typedef typename tools::promote_args::type result_type; typedef typename policies::evaluation::type value_type; diff --git a/include/boost/math/special_functions/ellint_3.hpp b/include/boost/math/special_functions/ellint_3.hpp index 2c22bf309d..5397f671cb 100644 --- a/include/boost/math/special_functions/ellint_3.hpp +++ b/include/boost/math/special_functions/ellint_3.hpp @@ -294,7 +294,7 @@ T ellint_pi_imp(T v, T k, T vc, const Policy& pol) if(v == 0) { - return (k == 0) ? boost::math::constants::pi() / 2 : ellint_k_imp(k, pol); + return (k == 0) ? boost::math::constants::pi() / 2 : boost::math::ellint_1(k, pol); } if(v < 0) @@ -308,7 +308,7 @@ T ellint_pi_imp(T v, T k, T vc, const Policy& pol) // This next part is split in two to avoid spurious over/underflow: result *= -v / (1 - v); result *= (1 - k2) / (k2 - v); - result += ellint_k_imp(k, pol) * k2 / (k2 - v); + result += boost::math::ellint_1(k, pol) * k2 / (k2 - v); return result; } @@ -343,17 +343,18 @@ inline typename tools::promote_args::type ellint_3(T1 k, T2 v, const Pol } // namespace detail template -inline typename tools::promote_args::type ellint_3(T1 k, T2 v, T3 phi, const Policy& pol) +inline typename tools::promote_args::type ellint_3(T1 k, T2 v, T3 phi, const Policy&) { typedef typename tools::promote_args::type result_type; typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise, policies::promote_double >::type forwarding_policy; return policies::checked_narrowing_cast( detail::ellint_pi_imp( static_cast(v), static_cast(phi), static_cast(k), static_cast(1-v), - pol), "boost::math::ellint_3<%1%>(%1%,%1%,%1%)"); + forwarding_policy()), "boost::math::ellint_3<%1%>(%1%,%1%,%1%)"); } template diff --git a/include/boost/math/special_functions/heuman_lambda.hpp b/include/boost/math/special_functions/heuman_lambda.hpp index 6389443ee4..a926745eb1 100644 --- a/include/boost/math/special_functions/heuman_lambda.hpp +++ b/include/boost/math/special_functions/heuman_lambda.hpp @@ -52,6 +52,11 @@ T heuman_lambda_imp(T phi, T k, const Policy& pol) } else { + typedef std::integral_constant::value&& std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 + > precision_tag_type; + T rkp = sqrt(kp); T ratio; if(rkp == 1) @@ -59,8 +64,8 @@ T heuman_lambda_imp(T phi, T k, const Policy& pol) return policies::raise_domain_error(function, "When 1-k^2 == 1 then phi must be < Pi/2, but got phi = %1%", phi, pol); } else - ratio = ellint_f_imp(phi, rkp, pol) / ellint_k_imp(rkp, pol); - result = ratio + ellint_k_imp(k, pol) * jacobi_zeta_imp(phi, rkp, pol) / constants::half_pi(); + ratio = ellint_f_imp(phi, rkp, pol) / ellint_k_imp(rkp, pol, precision_tag_type()); + result = ratio + ellint_k_imp(k, pol, precision_tag_type()) * jacobi_zeta_imp(phi, rkp, pol) / constants::half_pi(); } return result; } diff --git a/include/boost/math/special_functions/jacobi_zeta.hpp b/include/boost/math/special_functions/jacobi_zeta.hpp index 638049bc25..aca9ff6d44 100644 --- a/include/boost/math/special_functions/jacobi_zeta.hpp +++ b/include/boost/math/special_functions/jacobi_zeta.hpp @@ -49,7 +49,13 @@ T jacobi_zeta_imp(T phi, T k, const Policy& pol) if(k == 1) result = sinp * (boost::math::sign)(cosp); // We get here by simplifying JacobiZeta[w, 1] in Mathematica, and the fact that 0 <= phi. else - result = k2 * sinp * cosp * sqrt(1 - k2 * s2) * ellint_rj_imp(T(0), kp, T(1), T(1 - k2 * s2), pol) / (3 * ellint_k_imp(k, pol)); + { + typedef std::integral_constant::value&& std::numeric_limits::digits && (std::numeric_limits::digits <= 54) ? 0 : + std::is_floating_point::value && std::numeric_limits::digits && (std::numeric_limits::digits <= 64) ? 1 : 2 + > precision_tag_type; + result = k2 * sinp * cosp * sqrt(1 - k2 * s2) * ellint_rj_imp(T(0), kp, T(1), T(1 - k2 * s2), pol) / (3 * ellint_k_imp(k, pol, precision_tag_type())); + } return invert ? T(-result) : result; } diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index 035428c253..bded3acd8a 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -127,6 +127,17 @@ # endif #endif +#if !defined(BOOST_FORCEINLINE) +# if defined(_MSC_VER) +# define BOOST_FORCEINLINE __forceinline +# elif defined(__GNUC__) && __GNUC__ > 3 + // Clang also defines __GNUC__ (as 4) +# define BOOST_FORCEINLINE inline __attribute__ ((__always_inline__)) +# else +# define BOOST_FORCEINLINE inline +# endif +#endif + #endif // BOOST_MATH_STANDALONE // Support compilers with P0024R2 implemented without linking TBB diff --git a/reporting/performance/ellint_performance.cpp b/reporting/performance/ellint_performance.cpp new file mode 100644 index 0000000000..b975409483 --- /dev/null +++ b/reporting/performance/ellint_performance.cpp @@ -0,0 +1,432 @@ +// (C) Copyright John Maddock 2022. +// 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) + +// +// The purpose of this performance test is to probe the performance +// the of polynomial approximation to Elliptic integrals K and E. +// +// In the original paper: +// "Fast computation of complete elliptic integrals and Jacobian elliptic functions", +// Celestial Mechanics and Dynamical Astronomy, April 2012. +// they claim performance comparable to std::log, so we add comparison to log here +// as a reference "known good". +// +// We also measure the effect of disabling overflow errors, and of a "bare-bones" +// implementation which lacks a lot of our usual boilerplate. +// +// Note in the performance test routines, we had to sum the results of the function calls +// otherwise some msvc results were unrealistically fast, despite the use of +// benchmark::DoNotOptimize. +// +// Some sample results from msvc, and taking log(double) as a score of 1: +// +// ellint_1_performance 1.7 +// ellint_1_performance_no_overflow_check 1.6 +// ellint_2_performance 1.8 +// ellint_2_performance_no_overflow_check 1.6 +// basic_ellint_rational_performance 1.6 +// +// We can in fact get basic_ellint_rational_performance to much the same performance as log +// ONLY if we remove all error hangling for cases with m > 0.9. In particular the code appears +// to be ultra-sensitive to the presence of "if" statements which significantly hamper optimisation. +// +// Performance with gcc-cygwin appears to be broasly similar. +// +#include +#include +#include +#include +#include + +using boost::math::generate_random_uniform_vector; + +template +const std::vector& get_test_data() +{ + static const std::vector data = generate_random_uniform_vector(1024, 12345, 0.0, 0.9); + return data; +} +template +Real& get_result_data() +{ + static Real data = 0; + return data; +} + + +template +BOOST_FORCEINLINE T basic_ellint_rational(T k) +{ + using namespace boost::math::tools; + + static const char* function = "boost::math::ellint_k<%1%>(%1%)"; + + T m = k * k; + switch (static_cast(m * 20)) + { + case 0: + case 1: + //if (m < 0.1) + { + constexpr T coef[] = + { + 1.591003453790792180, + 0.416000743991786912, + 0.245791514264103415, + 0.179481482914906162, + 0.144556057087555150, + 0.123200993312427711, + 0.108938811574293531, + 0.098853409871592910, + 0.091439629201749751, + 0.085842591595413900, + 0.081541118718303215, + 0.078199656811256481910 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.05); + } + case 2: + case 3: + //else if (m < 0.2) + { + constexpr T coef[] = + { + 1.635256732264579992, + 0.471190626148732291, + 0.309728410831499587, + 0.252208311773135699, + 0.226725623219684650, + 0.215774446729585976, + 0.213108771877348910, + 0.216029124605188282, + 0.223255831633057896, + 0.234180501294209925, + 0.248557682972264071, + 0.266363809892617521 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.15); + } + case 4: + case 5: + //else if (m < 0.3) + { + constexpr T coef[] = + { + 1.685750354812596043, + 0.541731848613280329, + 0.401524438390690257, + 0.369642473420889090, + 0.376060715354583645, + 0.405235887085125919, + 0.453294381753999079, + 0.520518947651184205, + 0.609426039204995055, + 0.724263522282908870, + 0.871013847709812357, + 1.057652872753547036 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.25); + } + case 6: + case 7: + //else if (m < 0.4) + { + constexpr T coef[] = + { + 1.744350597225613243, + 0.634864275371935304, + 0.539842564164445538, + 0.571892705193787391, + 0.670295136265406100, + 0.832586590010977199, + 1.073857448247933265, + 1.422091460675497751, + 1.920387183402304829, + 2.632552548331654201, + 3.652109747319039160, + 5.115867135558865806, + 7.224080007363877411 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.35); + } + case 8: + case 9: + //else if (m < 0.5) + { + constexpr T coef[] = + { + 1.813883936816982644, + 0.763163245700557246, + 0.761928605321595831, + 0.951074653668427927, + 1.315180671703161215, + 1.928560693477410941, + 2.937509342531378755, + 4.594894405442878062, + 7.330071221881720772, + 11.87151259742530180, + 19.45851374822937738, + 32.20638657246426863, + 53.73749198700554656, + 90.27388602940998849 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.45); + } + case 10: + case 11: + //else if (m < 0.6) + { + constexpr T coef[] = + { + 1.898924910271553526, + 0.950521794618244435, + 1.151077589959015808, + 1.750239106986300540, + 2.952676812636875180, + 5.285800396121450889, + 9.832485716659979747, + 18.78714868327559562, + 36.61468615273698145, + 72.45292395127771801, + 145.1079577347069102, + 293.4786396308497026, + 598.3851815055010179, + 1228.420013075863451, + 2536.529755382764488 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.55); + } + case 12: + case 13: + //else if (m < 0.7) + { + constexpr T coef[] = + { + 2.007598398424376302, + 1.248457231212347337, + 1.926234657076479729, + 3.751289640087587680, + 8.119944554932045802, + 18.66572130873555361, + 44.60392484291437063, + 109.5092054309498377, + 274.2779548232413480, + 697.5598008606326163, + 1795.716014500247129, + 4668.381716790389910, + 12235.76246813664335, + 32290.17809718320818, + 85713.07608195964685, + 228672.1890493117096, + 612757.2711915852774 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.65); + } + case 14: + case 15: + //else if (m < 0.8) + { + constexpr T coef[] = + { + 2.156515647499643235, + 1.791805641849463243, + 3.826751287465713147, + 10.38672468363797208, + 31.40331405468070290, + 100.9237039498695416, + 337.3268282632272897, + 1158.707930567827917, + 4060.990742193632092, + 14454.00184034344795, + 52076.66107599404803, + 189493.6591462156887, + 695184.5762413896145, + 2567994.048255284686, + 9541921.966748386322, + 35634927.44218076174, + 133669298.4612040871, + 503352186.6866284541, + 1901975729.538660119, + 7208915015.330103756 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.75); + } + case 16: + //else if (m < 0.85) + { + constexpr T coef[] = + { + 2.318122621712510589, + 2.616920150291232841, + 7.897935075731355823, + 30.50239715446672327, + 131.4869365523528456, + 602.9847637356491617, + 2877.024617809972641, + 14110.51991915180325, + 70621.44088156540229, + 358977.2665825309926, + 1847238.263723971684, + 9600515.416049214109, + 50307677.08502366879, + 265444188.6527127967, + 1408862325.028702687, + 7515687935.373774627 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.825); + } + case 17: + //else if (m < 0.90) + { + constexpr T coef[] = + { + 2.473596173751343912, + 3.727624244118099310, + 15.60739303554930496, + 84.12850842805887747, + 506.9818197040613935, + 3252.277058145123644, + 21713.24241957434256, + 149037.0451890932766, + 1043999.331089990839, + 7427974.817042038995, + 53503839.67558661151, + 389249886.9948708474, + 2855288351.100810619, + 21090077038.76684053, + 156699833947.7902014, + 1170222242422.439893, + 8777948323668.937971, + 66101242752484.95041, + 499488053713388.7989, + 37859743397240299.20 + }; + return boost::math::tools::evaluate_polynomial(coef, m - 0.875); + } + case 18: + case 19: + return boost::math::detail::ellint_k_imp(k, boost::math::policies::policy<>(), std::integral_constant()); + default: + if (m == 1) + { + return boost::math::policies::raise_overflow_error(function, nullptr, boost::math::policies::policy<>()); + } + else + return boost::math::policies::raise_domain_error(function, + "Got k = %1%, function requires |k| <= 1", k, boost::math::policies::policy<>()); + } +} + +template +void basic_ellint_rational_performance(benchmark::State& state) +{ + std::vectorconst& test_set = get_test_data(); + Real& r = get_result_data(); + + for(auto _ : state) + { + for(unsigned i = 0; i < test_set.size(); ++i) + benchmark::DoNotOptimize(r += basic_ellint_rational(test_set[i])); + } +} + +template +void basic_ellint_rational_performance_no_error_check(benchmark::State& state) +{ + std::vectorconst& test_set = get_test_data(); + Real& r = get_result_data(); + + for(auto _ : state) + { + for(unsigned i = 0; i < test_set.size(); ++i) + benchmark::DoNotOptimize(r += basic_ellint_rational_no_error_checks(test_set[i])); + } +} + +template +void ellint_1_performance(benchmark::State& state) +{ + std::vectorconst& test_set = get_test_data(); + Real& r = get_result_data(); + + for(auto _ : state) + { + for(unsigned i = 0; i < test_set.size(); ++i) + benchmark::DoNotOptimize(r += boost::math::ellint_1(test_set[i])); + } +} + +template +void ellint_1_performance_no_overflow_check(benchmark::State& state) +{ + std::vectorconst& test_set = get_test_data(); + Real& r = get_result_data(); + + for(auto _ : state) + { + for(unsigned i = 0; i < test_set.size(); ++i) + benchmark::DoNotOptimize(r += boost::math::ellint_1(test_set[i], boost::math::policies::make_policy(boost::math::policies::overflow_error()))); + } +} + +template +void ellint_2_performance(benchmark::State& state) +{ + std::vectorconst& test_set = get_test_data(); + Real& r = get_result_data(); + + for(auto _ : state) + { + for(unsigned i = 0; i < test_set.size(); ++i) + benchmark::DoNotOptimize(r += boost::math::ellint_2(test_set[i])); + } +} + +template +void ellint_2_performance_no_overflow_check(benchmark::State& state) +{ + std::vectorconst& test_set = get_test_data(); + Real& r = get_result_data(); + + for(auto _ : state) + { + for(unsigned i = 0; i < test_set.size(); ++i) + benchmark::DoNotOptimize(r += boost::math::ellint_2(test_set[i], boost::math::policies::make_policy(boost::math::policies::overflow_error()))); + } +} + +template +void log_performance(benchmark::State& state) +{ + std::vectorconst& test_set = get_test_data(); + Real& r = get_result_data(); + + for(auto _ : state) + { + for(unsigned i = 0; i < test_set.size(); ++i) + benchmark::DoNotOptimize(r += log(test_set[i])); + } +} + +BENCHMARK_TEMPLATE(ellint_1_performance, float); +BENCHMARK_TEMPLATE(ellint_1_performance, double); +BENCHMARK_TEMPLATE(ellint_1_performance, long double); +BENCHMARK_TEMPLATE(ellint_1_performance_no_overflow_check, float); +BENCHMARK_TEMPLATE(ellint_1_performance_no_overflow_check, double); +BENCHMARK_TEMPLATE(ellint_1_performance_no_overflow_check, long double); +BENCHMARK_TEMPLATE(ellint_2_performance, float); +BENCHMARK_TEMPLATE(ellint_2_performance, double); +BENCHMARK_TEMPLATE(ellint_2_performance, long double); +BENCHMARK_TEMPLATE(ellint_2_performance_no_overflow_check, float); +BENCHMARK_TEMPLATE(ellint_2_performance_no_overflow_check, double); +BENCHMARK_TEMPLATE(ellint_2_performance_no_overflow_check, long double); +BENCHMARK_TEMPLATE(log_performance, float); +BENCHMARK_TEMPLATE(log_performance, double); +BENCHMARK_TEMPLATE(log_performance, long double); +BENCHMARK_TEMPLATE(basic_ellint_rational_performance, float); +BENCHMARK_TEMPLATE(basic_ellint_rational_performance, double); +BENCHMARK_TEMPLATE(basic_ellint_rational_performance, long double); + +BENCHMARK_MAIN();