1 // Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock 2 // Use, modification and distribution are subject to the 3 // Boost Software License, Version 1.0. (See accompanying file 4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 5 // 6 // History: 7 // XZ wrote the original of this file as part of the Google 8 // Summer of Code 2006. JM modified it to fit into the 9 // Boost.Math conceptual framework better, and to handle 10 // types longer than 80-bit reals. 11 // Updated 2015 to use Carlson's latest methods. 12 // 13 #ifndef BOOST_MATH_ELLINT_RF_HPP 14 #define BOOST_MATH_ELLINT_RF_HPP 15 16 #ifdef _MSC_VER 17 #pragma once 18 #endif 19 20 #include <boost/math/special_functions/math_fwd.hpp> 21 #include <boost/math/tools/config.hpp> 22 #include <boost/math/constants/constants.hpp> 23 #include <boost/math/policies/error_handling.hpp> 24 #include <boost/math/special_functions/ellint_rc.hpp> 25 26 // Carlson's elliptic integral of the first kind 27 // R_F(x, y, z) = 0.5 * \int_{0}^{\infty} [(t+x)(t+y)(t+z)]^{-1/2} dt 28 // Carlson, Numerische Mathematik, vol 33, 1 (1979) 29 30 namespace boost { namespace math { namespace detail{ 31 32 template <typename T, typename Policy> ellint_rf_imp(T x,T y,T z,const Policy & pol)33 T ellint_rf_imp(T x, T y, T z, const Policy& pol) 34 { 35 BOOST_MATH_STD_USING 36 using namespace boost::math; 37 using std::swap; 38 39 static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; 40 41 if(x < 0 || y < 0 || z < 0) 42 { 43 return policies::raise_domain_error<T>(function, 44 "domain error, all arguments must be non-negative, " 45 "only sensible result is %1%.", 46 std::numeric_limits<T>::quiet_NaN(), pol); 47 } 48 if(x + y == 0 || y + z == 0 || z + x == 0) 49 { 50 return policies::raise_domain_error<T>(function, 51 "domain error, at most one argument can be zero, " 52 "only sensible result is %1%.", 53 std::numeric_limits<T>::quiet_NaN(), pol); 54 } 55 // 56 // Special cases from http://dlmf.nist.gov/19.20#i 57 // 58 if(x == y) 59 { 60 if(x == z) 61 { 62 // x, y, z equal: 63 return 1 / sqrt(x); 64 } 65 else 66 { 67 // 2 equal, x and y: 68 if(z == 0) 69 return constants::pi<T>() / (2 * sqrt(x)); 70 else 71 return ellint_rc_imp(z, x, pol); 72 } 73 } 74 if(x == z) 75 { 76 if(y == 0) 77 return constants::pi<T>() / (2 * sqrt(x)); 78 else 79 return ellint_rc_imp(y, x, pol); 80 } 81 if(y == z) 82 { 83 if(x == 0) 84 return constants::pi<T>() / (2 * sqrt(y)); 85 else 86 return ellint_rc_imp(x, y, pol); 87 } 88 if(x == 0) 89 swap(x, z); 90 else if(y == 0) 91 swap(y, z); 92 if(z == 0) 93 { 94 // 95 // Special case for one value zero: 96 // 97 T xn = sqrt(x); 98 T yn = sqrt(y); 99 100 while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn)) 101 { 102 T t = sqrt(xn * yn); 103 xn = (xn + yn) / 2; 104 yn = t; 105 } 106 return constants::pi<T>() / (xn + yn); 107 } 108 109 T xn = x; 110 T yn = y; 111 T zn = z; 112 T An = (x + y + z) / 3; 113 T A0 = An; 114 T Q = pow(3 * boost::math::tools::epsilon<T>(), T(-1) / 8) * (std::max)((std::max)(fabs(An - xn), fabs(An - yn)), fabs(An - zn)); 115 T fn = 1; 116 117 118 // duplication 119 unsigned k = 1; 120 for(; k < boost::math::policies::get_max_series_iterations<Policy>(); ++k) 121 { 122 T root_x = sqrt(xn); 123 T root_y = sqrt(yn); 124 T root_z = sqrt(zn); 125 T lambda = root_x * root_y + root_x * root_z + root_y * root_z; 126 An = (An + lambda) / 4; 127 xn = (xn + lambda) / 4; 128 yn = (yn + lambda) / 4; 129 zn = (zn + lambda) / 4; 130 Q /= 4; 131 fn *= 4; 132 if(Q < fabs(An)) 133 break; 134 } 135 // Check to see if we gave up too soon: 136 policies::check_series_iterations<T>(function, k, pol); 137 BOOST_MATH_INSTRUMENT_VARIABLE(k); 138 139 T X = (A0 - x) / (An * fn); 140 T Y = (A0 - y) / (An * fn); 141 T Z = -X - Y; 142 143 // Taylor series expansion to the 7th order 144 T E2 = X * Y - Z * Z; 145 T E3 = X * Y * Z; 146 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); 147 } 148 149 } // namespace detail 150 151 template <class T1, class T2, class T3, class Policy> 152 inline typename tools::promote_args<T1, T2, T3>::type ellint_rf(T1 x,T2 y,T3 z,const Policy & pol)153 ellint_rf(T1 x, T2 y, T3 z, const Policy& pol) 154 { 155 typedef typename tools::promote_args<T1, T2, T3>::type result_type; 156 typedef typename policies::evaluation<result_type, Policy>::type value_type; 157 return policies::checked_narrowing_cast<result_type, Policy>( 158 detail::ellint_rf_imp( 159 static_cast<value_type>(x), 160 static_cast<value_type>(y), 161 static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"); 162 } 163 164 template <class T1, class T2, class T3> 165 inline typename tools::promote_args<T1, T2, T3>::type ellint_rf(T1 x,T2 y,T3 z)166 ellint_rf(T1 x, T2 y, T3 z) 167 { 168 return ellint_rf(x, y, z, policies::policy<>()); 169 } 170 171 }} // namespaces 172 173 #endif // BOOST_MATH_ELLINT_RF_HPP 174 175