diff options
Diffstat (limited to 'boost/math/special_functions/ellint_rg.hpp')
| -rw-r--r-- | boost/math/special_functions/ellint_rg.hpp | 136 |
1 files changed, 136 insertions, 0 deletions
diff --git a/boost/math/special_functions/ellint_rg.hpp b/boost/math/special_functions/ellint_rg.hpp new file mode 100644 index 0000000000..bb5b7c376a --- /dev/null +++ b/boost/math/special_functions/ellint_rg.hpp @@ -0,0 +1,136 @@ +// Copyright (c) 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) +// +#ifndef BOOST_MATH_ELLINT_RG_HPP +#define BOOST_MATH_ELLINT_RG_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/tools/config.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/policies/error_handling.hpp> +#include <boost/math/special_functions/ellint_rd.hpp> +#include <boost/math/special_functions/ellint_rf.hpp> +#include <boost/math/special_functions/pow.hpp> + +namespace boost { namespace math { namespace detail{ + + template <typename T, typename Policy> + T ellint_rg_imp(T x, T y, T z, const Policy& pol) + { + BOOST_MATH_STD_USING + static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; + + if(x < 0 || y < 0 || z < 0) + { + return policies::raise_domain_error<T>(function, + "domain error, all arguments must be non-negative, " + "only sensible result is %1%.", + std::numeric_limits<T>::quiet_NaN(), pol); + } + // + // Function is symmetric in x, y and z, but we require + // (x - z)(y - z) >= 0 to avoid cancellation error in the result + // which implies (for example) x >= z >= y + // + using std::swap; + if(x < y) + swap(x, y); + if(x < z) + swap(x, z); + if(y > z) + swap(y, z); + + BOOST_ASSERT(x >= z); + BOOST_ASSERT(z >= y); + // + // Special cases from http://dlmf.nist.gov/19.20#ii + // + if(x == z) + { + if(y == z) + { + // x = y = z + // This also works for x = y = z = 0 presumably. + return sqrt(x); + } + else if(y == 0) + { + // x = y, z = 0 + return constants::pi<T>() * sqrt(x) / 4; + } + else + { + // x = z, y != 0 + swap(x, y); + return (x == 0) ? T(sqrt(z) / 2) : T((z * ellint_rc_imp(x, z, pol) + sqrt(x)) / 2); + } + } + else if(y == z) + { + if(x == 0) + return constants::pi<T>() * sqrt(y) / 4; + else + return (y == 0) ? T(sqrt(x) / 2) : T((y * ellint_rc_imp(x, y, pol) + sqrt(x)) / 2); + } + else if(y == 0) + { + swap(y, z); + // + // Special handling for common case, from + // Numerical Computation of Real or Complex Elliptic Integrals, eq.46 + // + T xn = sqrt(x); + T yn = sqrt(y); + T x0 = xn; + T y0 = yn; + T sum = 0; + T sum_pow = 0.25f; + + while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn)) + { + T t = sqrt(xn * yn); + xn = (xn + yn) / 2; + yn = t; + sum_pow *= 2; + sum += sum_pow * boost::math::pow<2>(xn - yn); + } + T RF = constants::pi<T>() / (xn + yn); + return ((boost::math::pow<2>((x0 + y0) / 2) - sum) * RF) / 2; + } + return (z * ellint_rf_imp(x, y, z, pol) + - (x - z) * (y - z) * ellint_rd_imp(x, y, z, pol) / 3 + + sqrt(x * y / z)) / 2; + } + +} // namespace detail + +template <class T1, class T2, class T3, class Policy> +inline typename tools::promote_args<T1, T2, T3>::type + ellint_rg(T1 x, T2 y, T3 z, const Policy& pol) +{ + typedef typename tools::promote_args<T1, T2, T3>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<result_type, Policy>( + detail::ellint_rg_imp( + static_cast<value_type>(x), + static_cast<value_type>(y), + static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"); +} + +template <class T1, class T2, class T3> +inline typename tools::promote_args<T1, T2, T3>::type + ellint_rg(T1 x, T2 y, T3 z) +{ + return ellint_rg(x, y, z, policies::policy<>()); +} + +}} // namespaces + +#endif // BOOST_MATH_ELLINT_RG_HPP + |