113 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock2// Copyright (c) 2024 Matt Borland3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6//7// History:8// XZ wrote the original of this file as part of the Google9// Summer of Code 2006. JM modified it to fit into the10// Boost.Math conceptual framework better, and to correctly11// handle the y < 0 case.12// Updated 2015 to use Carlson's latest methods.13//14 15#ifndef BOOST_MATH_ELLINT_RC_HPP16#define BOOST_MATH_ELLINT_RC_HPP17 18#ifdef _MSC_VER19#pragma once20#endif21 22#include <boost/math/tools/config.hpp>23#include <boost/math/policies/error_handling.hpp>24#include <boost/math/special_functions/math_fwd.hpp>25#include <boost/math/special_functions/log1p.hpp>26#include <boost/math/constants/constants.hpp>27 28// Carlson's degenerate elliptic integral29// R_C(x, y) = R_F(x, y, y) = 0.5 * \int_{0}^{\infty} (t+x)^{-1/2} (t+y)^{-1} dt30// Carlson, Numerische Mathematik, vol 33, 1 (1979)31 32namespace boost { namespace math { namespace detail{33 34template <typename T, typename Policy>35BOOST_MATH_GPU_ENABLED T ellint_rc_imp(T x, T y, const Policy& pol)36{37 BOOST_MATH_STD_USING38 39 constexpr auto function = "boost::math::ellint_rc<%1%>(%1%,%1%)";40 41 if(x < 0)42 {43 return policies::raise_domain_error<T>(function, "Argument x must be non-negative but got %1%", x, pol);44 }45 if(y == 0)46 {47 return policies::raise_domain_error<T>(function, "Argument y must not be zero but got %1%", y, pol);48 }49 50 // for y < 0, the integral is singular, return Cauchy principal value51 T prefix, result;52 if(y < 0)53 {54 prefix = sqrt(x / (x - y));55 x = x - y;56 y = -y;57 }58 else59 prefix = 1;60 61 if(x == 0)62 {63 result = constants::half_pi<T>() / sqrt(y);64 }65 else if(x == y)66 {67 result = 1 / sqrt(x);68 }69 else if(y > x)70 {71 result = atan(sqrt((y - x) / x)) / sqrt(y - x);72 }73 else74 {75 if(y / x > T(0.5))76 {77 T arg = sqrt((x - y) / x);78 result = (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * sqrt(x - y));79 }80 else81 {82 result = log((sqrt(x) + sqrt(x - y)) / sqrt(y)) / sqrt(x - y);83 }84 }85 return prefix * result;86}87 88} // namespace detail89 90template <class T1, class T2, class Policy>91BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type 92 ellint_rc(T1 x, T2 y, const Policy& pol)93{94 typedef typename tools::promote_args<T1, T2>::type result_type;95 typedef typename policies::evaluation<result_type, Policy>::type value_type;96 return policies::checked_narrowing_cast<result_type, Policy>(97 detail::ellint_rc_imp(98 static_cast<value_type>(x),99 static_cast<value_type>(y), pol), "boost::math::ellint_rc<%1%>(%1%,%1%)");100}101 102template <class T1, class T2>103BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type 104 ellint_rc(T1 x, T2 y)105{106 return ellint_rc(x, y, policies::policy<>());107}108 109}} // namespaces110 111#endif // BOOST_MATH_ELLINT_RC_HPP112 113