184 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang2// Copyright (c) 2006 John Maddock3// Copyright (c) 2024 Matt Borland4// Use, modification and distribution are subject to the5// Boost Software License, Version 1.0. (See accompanying file6// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)7//8// History:9// XZ wrote the original of this file as part of the Google10// Summer of Code 2006. JM modified it to fit into the11// Boost.Math conceptual framework better, and to ensure12// that the code continues to work no matter how many digits13// type T has.14 15#ifndef BOOST_MATH_ELLINT_D_HPP16#define BOOST_MATH_ELLINT_D_HPP17 18#ifdef _MSC_VER19#pragma once20#endif21 22#include <boost/math/tools/config.hpp>23#include <boost/math/tools/type_traits.hpp>24#include <boost/math/special_functions/math_fwd.hpp>25#include <boost/math/special_functions/ellint_rf.hpp>26#include <boost/math/special_functions/ellint_rd.hpp>27#include <boost/math/special_functions/ellint_rg.hpp>28#include <boost/math/constants/constants.hpp>29#include <boost/math/policies/error_handling.hpp>30#include <boost/math/tools/workaround.hpp>31#include <boost/math/special_functions/round.hpp>32 33// Elliptic integrals (complete and incomplete) of the second kind34// Carlson, Numerische Mathematik, vol 33, 1 (1979)35 36namespace boost { namespace math {37 38template <class T1, class T2, class Policy>39BOOST_MATH_GPU_ENABLED typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol);40 41namespace detail{42 43template <typename T, typename Policy>44BOOST_MATH_GPU_ENABLED T ellint_d_imp(T k, const Policy& pol);45 46// Elliptic integral (Legendre form) of the second kind47template <typename T, typename Policy>48BOOST_MATH_GPU_ENABLED T ellint_d_imp(T phi, T k, const Policy& pol)49{50 BOOST_MATH_STD_USING51 using namespace boost::math::tools;52 using namespace boost::math::constants;53 54 bool invert = false;55 if(phi < 0)56 {57 phi = fabs(phi);58 invert = true;59 }60 61 T result;62 63 if(phi >= tools::max_value<T>())64 {65 // Need to handle infinity as a special case:66 result = policies::raise_overflow_error<T>("boost::math::ellint_d<%1%>(%1%,%1%)", nullptr, pol);67 }68 else if(phi > 1 / tools::epsilon<T>())69 {70 // Phi is so large that phi%pi is necessarily zero (or garbage),71 // just return the second part of the duplication formula:72 result = 2 * phi * ellint_d_imp(k, pol) / constants::pi<T>();73 }74 else75 {76 // Carlson's algorithm works only for |phi| <= pi/2,77 // use the integrand's periodicity to normalize phi78 //79 T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>()));80 T m = boost::math::round((phi - rphi) / constants::half_pi<T>());81 int s = 1;82 if(boost::math::tools::fmod_workaround(m, T(2)) > T(0.5))83 {84 m += 1;85 s = -1;86 rphi = constants::half_pi<T>() - rphi;87 }88 BOOST_MATH_INSTRUMENT_VARIABLE(rphi);89 BOOST_MATH_INSTRUMENT_VARIABLE(m);90 T sinp = sin(rphi);91 T cosp = cos(rphi);92 BOOST_MATH_INSTRUMENT_VARIABLE(sinp);93 BOOST_MATH_INSTRUMENT_VARIABLE(cosp);94 T c = 1 / (sinp * sinp);95 T cm1 = cosp * cosp / (sinp * sinp); // c - 196 T k2 = k * k;97 if(k2 * sinp * sinp > 1)98 {99 return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol);100 }101 else if(rphi == 0)102 {103 result = 0;104 }105 else106 {107 // http://dlmf.nist.gov/19.25#E10108 result = s * ellint_rd_imp(cm1, T(c - k2), c, pol) / 3;109 BOOST_MATH_INSTRUMENT_VARIABLE(result);110 }111 if(m != 0)112 result += m * ellint_d_imp(k, pol);113 }114 return invert ? T(-result) : result;115}116 117// Complete elliptic integral (Legendre form) of the second kind118template <typename T, typename Policy>119BOOST_MATH_GPU_ENABLED T ellint_d_imp(T k, const Policy& pol)120{121 BOOST_MATH_STD_USING122 using namespace boost::math::tools;123 124 if (abs(k) >= 1)125 {126 return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%)", "Got k = %1%, function requires |k| <= 1", k, pol);127 }128 if(fabs(k) <= tools::root_epsilon<T>())129 return constants::pi<T>() / 4;130 131 T x = 0;132 T t = k * k;133 T y = 1 - t;134 T z = 1;135 T value = ellint_rd_imp(x, y, z, pol) / 3;136 137 return value;138}139 140template <typename T, typename Policy>141BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type ellint_d(T k, const Policy& pol, const boost::math::true_type&)142{143 typedef typename tools::promote_args<T>::type result_type;144 typedef typename policies::evaluation<result_type, Policy>::type value_type;145 return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_d_imp(static_cast<value_type>(k), pol), "boost::math::ellint_d<%1%>(%1%)");146}147 148// Elliptic integral (Legendre form) of the second kind149template <class T1, class T2>150BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const boost::math::false_type&)151{152 return boost::math::ellint_d(k, phi, policies::policy<>());153}154 155} // detail156 157// Complete elliptic integral (Legendre form) of the second kind158template <typename T>159BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type ellint_d(T k)160{161 return ellint_d(k, policies::policy<>());162}163 164// Elliptic integral (Legendre form) of the second kind165template <class T1, class T2>166BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi)167{168 typedef typename policies::is_policy<T2>::type tag_type;169 return detail::ellint_d(k, phi, tag_type());170}171 172template <class T1, class T2, class Policy>173BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol)174{175 typedef typename tools::promote_args<T1, T2>::type result_type;176 typedef typename policies::evaluation<result_type, Policy>::type value_type;177 return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_d_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%,%1%)");178}179 180}} // namespaces181 182#endif // BOOST_MATH_ELLINT_D_HPP183 184