X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=boost%2Fmath%2Fspecial_functions%2Fellint_1.hpp;fp=boost%2Fmath%2Fspecial_functions%2Fellint_1.hpp;h=e266b2e29055a1c194ae756951d119d937269285;hp=0000000000000000000000000000000000000000;hb=2d71eb92104693ca9baa5a2e1c23eeca776d8fd3;hpb=da57529b92adbb7ae74a89861cb39fb35ac7c62d diff --git a/boost/math/special_functions/ellint_1.hpp b/boost/math/special_functions/ellint_1.hpp new file mode 100644 index 0000000..e266b2e --- /dev/null +++ b/boost/math/special_functions/ellint_1.hpp @@ -0,0 +1,188 @@ +// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2006 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 ensure +// that the code continues to work no matter how many digits +// type T has. + +#ifndef BOOST_MATH_ELLINT_1_HPP +#define BOOST_MATH_ELLINT_1_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include +#include +#include +#include +#include + +// Elliptic integrals (complete and incomplete) of the first kind +// Carlson, Numerische Mathematik, vol 33, 1 (1979) + +namespace boost { namespace math { + +template +typename tools::promote_args::type ellint_1(T1 k, T2 phi, const Policy& pol); + +namespace detail{ + +template +T ellint_k_imp(T k, const Policy& pol); + +// Elliptic integral (Legendre form) of the first kind +template +T ellint_f_imp(T phi, T k, const Policy& pol) +{ + BOOST_MATH_STD_USING + using namespace boost::math::tools; + using namespace boost::math::constants; + + static const char* function = "boost::math::ellint_f<%1%>(%1%,%1%)"; + BOOST_MATH_INSTRUMENT_VARIABLE(phi); + BOOST_MATH_INSTRUMENT_VARIABLE(k); + BOOST_MATH_INSTRUMENT_VARIABLE(function); + + if (abs(k) > 1) + { + return policies::raise_domain_error(function, + "Got k = %1%, function requires |k| <= 1", k, pol); + } + + bool invert = false; + if(phi < 0) + { + BOOST_MATH_INSTRUMENT_VARIABLE(phi); + phi = fabs(phi); + invert = true; + } + + T result; + + if(phi >= tools::max_value()) + { + // Need to handle infinity as a special case: + result = policies::raise_overflow_error(function, 0, pol); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + } + else if(phi > 1 / tools::epsilon()) + { + // Phi is so large that phi%pi is necessarily zero (or garbage), + // just return the second part of the duplication formula: + result = 2 * phi * ellint_k_imp(k, pol) / constants::pi(); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + } + else + { + // Carlson's algorithm works only for |phi| <= pi/2, + // use the integrand's periodicity to normalize phi + // + // Xiaogang's original code used a cast to long long here + // but that fails if T has more digits than a long long, + // so rewritten to use fmod instead: + // + BOOST_MATH_INSTRUMENT_CODE("pi/2 = " << constants::pi() / 2); + T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi())); + BOOST_MATH_INSTRUMENT_VARIABLE(rphi); + T m = boost::math::round((phi - rphi) / constants::half_pi()); + BOOST_MATH_INSTRUMENT_VARIABLE(m); + int s = 1; + if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5) + { + m += 1; + s = -1; + rphi = constants::half_pi() - rphi; + BOOST_MATH_INSTRUMENT_VARIABLE(rphi); + } + T sinp = sin(rphi); + T cosp = cos(rphi); + BOOST_MATH_INSTRUMENT_VARIABLE(sinp); + BOOST_MATH_INSTRUMENT_VARIABLE(cosp); + result = s * sinp * ellint_rf_imp(T(cosp * cosp), T(1 - k * k * sinp * sinp), T(1), pol); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + if(m != 0) + { + result += m * ellint_k_imp(k, pol); + BOOST_MATH_INSTRUMENT_VARIABLE(result); + } + } + return invert ? T(-result) : result; +} + +// Complete elliptic integral (Legendre form) of the first kind +template +T ellint_k_imp(T k, const Policy& pol) +{ + BOOST_MATH_STD_USING + using namespace boost::math::tools; + + static const char* function = "boost::math::ellint_k<%1%>(%1%)"; + + if (abs(k) > 1) + { + return policies::raise_domain_error(function, + "Got k = %1%, function requires |k| <= 1", k, pol); + } + if (abs(k) == 1) + { + return policies::raise_overflow_error(function, 0, pol); + } + + T x = 0; + T y = 1 - k * k; + T z = 1; + T value = ellint_rf_imp(x, y, z, pol); + + return value; +} + +template +inline typename tools::promote_args::type ellint_1(T k, const Policy& pol, const mpl::true_&) +{ + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + return policies::checked_narrowing_cast(detail::ellint_k_imp(static_cast(k), pol), "boost::math::ellint_1<%1%>(%1%)"); +} + +template +inline typename tools::promote_args::type ellint_1(T1 k, T2 phi, const mpl::false_&) +{ + return boost::math::ellint_1(k, phi, policies::policy<>()); +} + +} + +// Complete elliptic integral (Legendre form) of the first kind +template +inline typename tools::promote_args::type ellint_1(T k) +{ + return ellint_1(k, policies::policy<>()); +} + +// Elliptic integral (Legendre form) of the first kind +template +inline typename tools::promote_args::type ellint_1(T1 k, T2 phi, 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_f_imp(static_cast(phi), static_cast(k), pol), "boost::math::ellint_1<%1%>(%1%,%1%)"); +} + +template +inline typename tools::promote_args::type ellint_1(T1 k, T2 phi) +{ + typedef typename policies::is_policy::type tag_type; + return detail::ellint_1(k, phi, tag_type()); +} + +}} // namespaces + +#endif // BOOST_MATH_ELLINT_1_HPP +