1 // Copyright (c) 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 #ifndef BOOST_MATH_ELLINT_RG_HPP 7 #define BOOST_MATH_ELLINT_RG_HPP 8 9 #ifdef _MSC_VER 10 #pragma once 11 #endif 12 13 #include <boost/math/special_functions/math_fwd.hpp> 14 #include <boost/math/tools/config.hpp> 15 #include <boost/math/constants/constants.hpp> 16 #include <boost/math/policies/error_handling.hpp> 17 #include <boost/math/special_functions/ellint_rd.hpp> 18 #include <boost/math/special_functions/ellint_rf.hpp> 19 #include <boost/math/special_functions/pow.hpp> 20 21 namespace boost { namespace math { namespace detail{ 22 23 template <typename T, typename Policy> 24 T ellint_rg_imp(T x, T y, T z, const Policy& pol) 25 { 26 BOOST_MATH_STD_USING 27 static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; 28 29 if(x < 0 || y < 0 || z < 0) 30 { 31 return policies::raise_domain_error<T>(function, 32 "domain error, all arguments must be non-negative, " 33 "only sensible result is %1%.", 34 std::numeric_limits<T>::quiet_NaN(), pol); 35 } 36 // 37 // Function is symmetric in x, y and z, but we require 38 // (x - z)(y - z) >= 0 to avoid cancellation error in the result 39 // which implies (for example) x >= z >= y 40 // 41 using std::swap; 42 if(x < y) 43 swap(x, y); 44 if(x < z) 45 swap(x, z); 46 if(y > z) 47 swap(y, z); 48 49 BOOST_ASSERT(x >= z); 50 BOOST_ASSERT(z >= y); 51 // 52 // Special cases from http://dlmf.nist.gov/19.20#ii 53 // 54 if(x == z) 55 { 56 if(y == z) 57 { 58 // x = y = z 59 // This also works for x = y = z = 0 presumably. 60 return sqrt(x); 61 } 62 else if(y == 0) 63 { 64 // x = y, z = 0 65 return constants::pi<T>() * sqrt(x) / 4; 66 } 67 else 68 { 69 // x = z, y != 0 70 swap(x, y); 71 return (x == 0) ? T(sqrt(z) / 2) : T((z * ellint_rc_imp(x, z, pol) + sqrt(x)) / 2); 72 } 73 } 74 else if(y == z) 75 { 76 if(x == 0) 77 return constants::pi<T>() * sqrt(y) / 4; 78 else 79 return (y == 0) ? T(sqrt(x) / 2) : T((y * ellint_rc_imp(x, y, pol) + sqrt(x)) / 2); 80 } 81 else if(y == 0) 82 { 83 swap(y, z); 84 // 85 // Special handling for common case, from 86 // Numerical Computation of Real or Complex Elliptic Integrals, eq.46 87 // 88 T xn = sqrt(x); 89 T yn = sqrt(y); 90 T x0 = xn; 91 T y0 = yn; 92 T sum = 0; 93 T sum_pow = 0.25f; 94 95 while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn)) 96 { 97 T t = sqrt(xn * yn); 98 xn = (xn + yn) / 2; 99 yn = t; 100 sum_pow *= 2; 101 sum += sum_pow * boost::math::pow<2>(xn - yn); 102 } 103 T RF = constants::pi<T>() / (xn + yn); 104 return ((boost::math::pow<2>((x0 + y0) / 2) - sum) * RF) / 2; 105 } 106 return (z * ellint_rf_imp(x, y, z, pol) 107 - (x - z) * (y - z) * ellint_rd_imp(x, y, z, pol) / 3 108 + sqrt(x * y / z)) / 2; 109 } 110 111 } // namespace detail 112 113 template <class T1, class T2, class T3, class Policy> 114 inline typename tools::promote_args<T1, T2, T3>::type ellint_rg(T1 x,T2 y,T3 z,const Policy & pol)115 ellint_rg(T1 x, T2 y, T3 z, const Policy& pol) 116 { 117 typedef typename tools::promote_args<T1, T2, T3>::type result_type; 118 typedef typename policies::evaluation<result_type, Policy>::type value_type; 119 return policies::checked_narrowing_cast<result_type, Policy>( 120 detail::ellint_rg_imp( 121 static_cast<value_type>(x), 122 static_cast<value_type>(y), 123 static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"); 124 } 125 126 template <class T1, class T2, class T3> 127 inline typename tools::promote_args<T1, T2, T3>::type ellint_rg(T1 x,T2 y,T3 z)128 ellint_rg(T1 x, T2 y, T3 z) 129 { 130 return ellint_rg(x, y, z, policies::policy<>()); 131 } 132 133 }} // namespaces 134 135 #endif // BOOST_MATH_ELLINT_RG_HPP 136 137