/usr/include/boost/math/special_functions
Edit: /usr/include/boost/math/special_functions/ellint_rc.hpp (3135B)
// Copyright (c) 2006 Xiaogang Zhang
// 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)
//
// History:
// XZ wrote the original of this file as part of the Google
// Summer of Code 2006. JM modified it to fit into the
// Boost.Math conceptual framework better, and to correctly
// handle the y < 0 case.
// Updated 2015 to use Carlson's latest methods.
//
#ifndef BOOST_MATH_ELLINT_RC_HPP
#define BOOST_MATH_ELLINT_RC_HPP
#ifdef _MSC_VER
#pragma once
#endif
#include
#include
#include
#include
#include
#include
// Carlson's degenerate elliptic integral
// R_C(x, y) = R_F(x, y, y) = 0.5 * \int_{0}^{\infty} (t+x)^{-1/2} (t+y)^{-1} dt
// Carlson, Numerische Mathematik, vol 33, 1 (1979)
namespace boost { namespace math { namespace detail{
template
T ellint_rc_imp(T x, T y, const Policy& pol)
{
BOOST_MATH_STD_USING
static const char* function = "boost::math::ellint_rc<%1%>(%1%,%1%)";
if(x < 0)
{
return policies::raise_domain_error(function,
"Argument x must be non-negative but got %1%", x, pol);
}
if(y == 0)
{
return policies::raise_domain_error(function,
"Argument y must not be zero but got %1%", y, pol);
}
// for y < 0, the integral is singular, return Cauchy principal value
T prefix, result;
if(y < 0)
{
prefix = sqrt(x / (x - y));
x = x - y;
y = -y;
}
else
prefix = 1;
if(x == 0)
{
result = constants::half_pi() / sqrt(y);
}
else if(x == y)
{
result = 1 / sqrt(x);
}
else if(y > x)
{
result = atan(sqrt((y - x) / x)) / sqrt(y - x);
}
else
{
if(y / x > 0.5)
{
T arg = sqrt((x - y) / x);
result = (boost::math::log1p(arg) - boost::math::log1p(-arg)) / (2 * sqrt(x - y));
}
else
{
result = log((sqrt(x) + sqrt(x - y)) / sqrt(y)) / sqrt(x - y);
}
}
return prefix * result;
}
} // namespace detail
template
inline typename tools::promote_args::type
ellint_rc(T1 x, T2 y, const Policy& pol)
{
typedef typename tools::promote_args::type result_type;
typedef typename policies::evaluation::type value_type;
return policies::checked_narrowing_cast(
detail::ellint_rc_imp(
static_cast(x),
static_cast(y), pol), "boost::math::ellint_rc<%1%>(%1%,%1%)");
}
template
inline typename tools::promote_args::type
ellint_rc(T1 x, T2 y)
{
return ellint_rc(x, y, policies::policy<>());
}
}} // namespaces
#endif // BOOST_MATH_ELLINT_RC_HPP