summaryrefslogtreecommitdiff log msg author committer range
diff options
 context: 12345678910152025303540 space: includeignore mode: unifiedssdiffstat only
Diffstat (limited to 'boost/math/special_functions/ellint_rf.hpp')
-rw-r--r--boost/math/special_functions/ellint_rf.hpp182
1 files changed, 112 insertions, 70 deletions
 diff --git a/boost/math/special_functions/ellint_rf.hpp b/boost/math/special_functions/ellint_rf.hppindex ac57257..a8a7b4b 100644--- a/boost/math/special_functions/ellint_rf.hpp+++ b/boost/math/special_functions/ellint_rf.hpp@@ -1,4 +1,4 @@-// Copyright (c) 2006 Xiaogang Zhang+// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock // 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)@@ -8,6 +8,7 @@ // Summer of Code 2006. JM modified it to fit into the // Boost.Math conceptual framework better, and to handle // types longer than 80-bit reals.+// Updated 2015 to use Carlson's latest methods. // #ifndef BOOST_MATH_ELLINT_RF_HPP #define BOOST_MATH_ELLINT_RF_HPP@@ -18,8 +19,9 @@ #include #include -+#include #include +#include // Carlson's elliptic integral of the first kind // R_F(x, y, z) = 0.5 * \int_{0}^{\infty} [(t+x)(t+y)(t+z)]^{-1/2} dt@@ -27,82 +29,122 @@ namespace boost { namespace math { namespace detail{ -template -T ellint_rf_imp(T x, T y, T z, const Policy& pol)-{- T value, X, Y, Z, E2, E3, u, lambda, tolerance;- unsigned long k;-- BOOST_MATH_STD_USING- using namespace boost::math::tools;+ template + T ellint_rf_imp(T x, T y, T z, const Policy& pol)+ {+ BOOST_MATH_STD_USING+ using namespace boost::math;+ using std::swap; - static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)";+ static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; - if (x < 0 || y < 0 || z < 0)- {- return policies::raise_domain_error(function,+ if(x < 0 || y < 0 || z < 0)+ {+ return policies::raise_domain_error(function, "domain error, all arguments must be non-negative, " "only sensible result is %1%.", std::numeric_limits::quiet_NaN(), pol);- }- if (x + y == 0 || y + z == 0 || z + x == 0)- {- return policies::raise_domain_error(function,+ }+ if(x + y == 0 || y + z == 0 || z + x == 0)+ {+ return policies::raise_domain_error(function, "domain error, at most one argument can be zero, " "only sensible result is %1%.", std::numeric_limits::quiet_NaN(), pol);- }-- // Carlson scales error as the 6th power of tolerance,- // but this seems not to work for types larger than- // 80-bit reals, this heuristic seems to work OK:- if(policies::digits() > 64)- {- tolerance = pow(tools::epsilon(), T(1)/4.25f);- BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);- }- else- {- tolerance = pow(4*tools::epsilon(), T(1)/6);- BOOST_MATH_INSTRUMENT_VARIABLE(tolerance);- }-- // duplication- k = 1;- do- {- u = (x + y + z) / 3;- X = (u - x) / u;- Y = (u - y) / u;- Z = (u - z) / u;-- // Termination condition: - if ((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance) - break; -- T sx = sqrt(x);- T sy = sqrt(y);- T sz = sqrt(z);- lambda = sy * (sx + sz) + sz * sx;- x = (x + lambda) / 4;- y = (y + lambda) / 4;- z = (z + lambda) / 4;- ++k;- }- while(k < policies::get_max_series_iterations());-- // Check to see if we gave up too soon:- policies::check_series_iterations(function, k, pol);- BOOST_MATH_INSTRUMENT_VARIABLE(k);-- // Taylor series expansion to the 5th order- E2 = X * Y - Z * Z;- E3 = X * Y * Z;- value = (1 + E2*(E2/24 - E3*T(3)/44 - T(0.1)) + E3/14) / sqrt(u);- BOOST_MATH_INSTRUMENT_VARIABLE(value);-- return value;-}+ }+ //+ // Special cases from http://dlmf.nist.gov/19.20#i+ //+ if(x == y)+ {+ if(x == z)+ {+ // x, y, z equal:+ return 1 / sqrt(x);+ }+ else+ {+ // 2 equal, x and y:+ if(z == 0)+ return constants::pi() / (2 * sqrt(x));+ else+ return ellint_rc_imp(z, x, pol);+ }+ }+ if(x == z)+ {+ if(y == 0)+ return constants::pi() / (2 * sqrt(x));+ else+ return ellint_rc_imp(y, x, pol);+ }+ if(y == z)+ {+ if(x == 0)+ return constants::pi() / (2 * sqrt(y));+ else+ return ellint_rc_imp(x, y, pol);+ }+ if(x == 0)+ swap(x, z);+ else if(y == 0)+ swap(y, z);+ if(z == 0)+ {+ //+ // Special case for one value zero:+ //+ T xn = sqrt(x);+ T yn = sqrt(y);++ while(fabs(xn - yn) >= 2.7 * tools::root_epsilon() * fabs(xn))+ {+ T t = sqrt(xn * yn);+ xn = (xn + yn) / 2;+ yn = t;+ }+ return constants::pi() / (xn + yn);+ }++ T xn = x;+ T yn = y;+ T zn = z;+ T An = (x + y + z) / 3;+ T A0 = An;+ T Q = pow(3 * boost::math::tools::epsilon(), T(-1) / 8) * (std::max)((std::max)(fabs(An - xn), fabs(An - yn)), fabs(An - zn));+ T fn = 1;+++ // duplication+ unsigned k = 1;+ for(; k < boost::math::policies::get_max_series_iterations(); ++k)+ {+ T root_x = sqrt(xn);+ T root_y = sqrt(yn);+ T root_z = sqrt(zn);+ T lambda = root_x * root_y + root_x * root_z + root_y * root_z;+ An = (An + lambda) / 4;+ xn = (xn + lambda) / 4;+ yn = (yn + lambda) / 4;+ zn = (zn + lambda) / 4;+ Q /= 4;+ fn *= 4;+ if(Q < fabs(An))+ break;+ }+ // Check to see if we gave up too soon:+ policies::check_series_iterations(function, k, pol);+ BOOST_MATH_INSTRUMENT_VARIABLE(k);++ T X = (A0 - x) / (An * fn);+ T Y = (A0 - y) / (An * fn);+ T Z = -X - Y;++ // Taylor series expansion to the 7th order+ T E2 = X * Y - Z * Z;+ T E3 = X * Y * Z;+ return (1 + E3 * (T(1) / 14 + 3 * E3 / 104) + E2 * (T(-1) / 10 + E2 / 24 - (3 * E3) / 44 - 5 * E2 * E2 / 208 + E2 * E3 / 16)) / sqrt(An);+ } } // namespace detail