131 lines · plain
1// Copyright (c) 2015 John Maddock2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5//6#ifndef BOOST_MATH_ELLINT_RG_HPP7#define BOOST_MATH_ELLINT_RG_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/special_functions/math_fwd.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 21namespace boost { namespace math { namespace detail{22 23 template <typename T, typename Policy>24 BOOST_MATH_GPU_ENABLED T ellint_rg_imp(T x, T y, T z, const Policy& pol)25 {26 BOOST_MATH_STD_USING27 constexpr auto 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, "domain error, all arguments must be non-negative, only sensible result is %1%.", boost::math::numeric_limits<T>::quiet_NaN(), pol);32 }33 //34 // Function is symmetric in x, y and z, but we require35 // (x - z)(y - z) >= 0 to avoid cancellation error in the result36 // which implies (for example) x >= z >= y37 //38 if(x < y)39 BOOST_MATH_GPU_SAFE_SWAP(x, y);40 if(x < z)41 BOOST_MATH_GPU_SAFE_SWAP(x, z);42 if(y > z)43 BOOST_MATH_GPU_SAFE_SWAP(y, z);44 45 BOOST_MATH_ASSERT(x >= z);46 BOOST_MATH_ASSERT(z >= y);47 //48 // Special cases from http://dlmf.nist.gov/19.20#ii49 //50 if(x == z)51 {52 if(y == z)53 {54 // x = y = z55 // This also works for x = y = z = 0 presumably.56 return sqrt(x);57 }58 else if(y == 0)59 {60 // x = y, z = 061 return constants::pi<T>() * sqrt(x) / 4;62 }63 else64 {65 // x = z, y != 066 BOOST_MATH_GPU_SAFE_SWAP(x, y);67 return (x == 0) ? T(sqrt(z) / 2) : T((z * ellint_rc_imp(x, z, pol) + sqrt(x)) / 2);68 }69 }70 else if(y == z)71 {72 BOOST_MATH_ASSERT(x > 0); // Ordering of x,y,z above takes care of x == 0 case.73 return (y == 0) ? T(sqrt(x) / 2) : T((y * ellint_rc_imp(x, y, pol) + sqrt(x)) / 2);74 }75 else if(y == 0)76 {77 BOOST_MATH_GPU_SAFE_SWAP(y, z);78 //79 // Special handling for common case, from80 // Numerical Computation of Real or Complex Elliptic Integrals, eq.4681 //82 T xn = sqrt(x);83 T yn = sqrt(y);84 T x0 = xn;85 T y0 = yn;86 T sum = 0;87 T sum_pow = 0.25f;88 89 while(fabs(xn - yn) >= T(2.7) * tools::root_epsilon<T>() * fabs(xn))90 {91 T t = sqrt(xn * yn);92 xn = (xn + yn) / 2;93 yn = t;94 sum_pow *= 2;95 sum += sum_pow * boost::math::pow<2>(xn - yn);96 }97 T RF = constants::pi<T>() / (xn + yn);98 return ((boost::math::pow<2>((x0 + y0) / 2) - sum) * RF) / 2;99 }100 return (z * ellint_rf_imp(x, y, z, pol)101 - (x - z) * (y - z) * ellint_rd_imp(x, y, z, pol) / 3102 + sqrt(x * y / z)) / 2;103 }104 105} // namespace detail106 107template <class T1, class T2, class T3, class Policy>108BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type 109 ellint_rg(T1 x, T2 y, T3 z, const Policy& pol)110{111 typedef typename tools::promote_args<T1, T2, T3>::type result_type;112 typedef typename policies::evaluation<result_type, Policy>::type value_type;113 return policies::checked_narrowing_cast<result_type, Policy>(114 detail::ellint_rg_imp(115 static_cast<value_type>(x),116 static_cast<value_type>(y),117 static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)");118}119 120template <class T1, class T2, class T3>121BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type 122 ellint_rg(T1 x, T2 y, T3 z)123{124 return ellint_rg(x, y, z, policies::policy<>());125}126 127}} // namespaces128 129#endif // BOOST_MATH_ELLINT_RG_HPP130 131