/usr/include/boost/math/special_functions
Edit: /usr/include/boost/math/special_functions/ellint_rg.hpp (4125B)
// 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
#include
#include
#include
#include
#include
#include
namespace boost { namespace math { namespace detail{
template
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(function,
"domain error, all arguments must be non-negative, "
"only sensible result is %1%.",
std::numeric_limits::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() * 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() * 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() * 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() / (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
inline typename tools::promote_args::type
ellint_rg(T1 x, T2 y, T3 z, 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_rg_imp(
static_cast(x),
static_cast(y),
static_cast(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)");
}
template
inline typename tools::promote_args::type
ellint_rg(T1 x, T2 y, T3 z)
{
return ellint_rg(x, y, z, policies::policy<>());
}
}} // namespaces
#endif // BOOST_MATH_ELLINT_RG_HPP