208 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock.2// 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 slightly to fit into the10// Boost.Math conceptual framework better.11// Updated 2015 to use Carlson's latest methods.12 13#ifndef BOOST_MATH_ELLINT_RD_HPP14#define BOOST_MATH_ELLINT_RD_HPP15 16#ifdef _MSC_VER17#pragma once18#endif19 20#include <boost/math/tools/config.hpp>21#include <boost/math/tools/promotion.hpp>22#include <boost/math/special_functions/math_fwd.hpp>23#include <boost/math/special_functions/ellint_rc.hpp>24#include <boost/math/policies/error_handling.hpp>25 26// Carlson's elliptic integral of the second kind27// R_D(x, y, z) = R_J(x, y, z, z) = 1.5 * \int_{0}^{\infty} [(t+x)(t+y)]^{-1/2} (t+z)^{-3/2} dt28// Carlson, Numerische Mathematik, vol 33, 1 (1979)29 30namespace boost { namespace math { namespace detail{31 32template <typename T, typename Policy>33BOOST_MATH_GPU_ENABLED T ellint_rd_imp(T x, T y, T z, const Policy& pol)34{35 BOOST_MATH_STD_USING36 37 constexpr auto function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)";38 39 if(x < 0)40 {41 return policies::raise_domain_error<T>(function, "Argument x must be >= 0, but got %1%", x, pol);42 }43 if(y < 0)44 {45 return policies::raise_domain_error<T>(function, "Argument y must be >= 0, but got %1%", y, pol);46 }47 if(z <= 0)48 {49 return policies::raise_domain_error<T>(function, "Argument z must be > 0, but got %1%", z, pol);50 }51 if(x + y == 0)52 {53 return policies::raise_domain_error<T>(function, "At most one argument can be zero, but got, x + y = %1%", x + y, pol);54 }55 //56 // Special cases from http://dlmf.nist.gov/19.20#iv57 //58 59 if(x == z)60 {61 BOOST_MATH_GPU_SAFE_SWAP(x, y);62 }63 if(y == z)64 {65 if(x == y)66 {67 return 1 / (x * sqrt(x));68 }69 else if(x == 0)70 {71 return 3 * constants::pi<T>() / (4 * y * sqrt(y));72 }73 else74 {75 if(BOOST_MATH_GPU_SAFE_MAX(x, y) / BOOST_MATH_GPU_SAFE_MIN(x, y) > T(1.3))76 return 3 * (ellint_rc_imp(x, y, pol) - sqrt(x) / y) / (2 * (y - x));77 // Otherwise fall through to avoid cancellation in the above (RC(x,y) -> 1/x^0.5 as x -> y)78 }79 }80 if(x == y)81 {82 if(BOOST_MATH_GPU_SAFE_MAX(x, z) / BOOST_MATH_GPU_SAFE_MIN(x, z) > T(1.3))83 return 3 * (ellint_rc_imp(z, x, pol) - 1 / sqrt(z)) / (z - x);84 // Otherwise fall through to avoid cancellation in the above (RC(x,y) -> 1/x^0.5 as x -> y)85 }86 if(y == 0)87 {88 BOOST_MATH_GPU_SAFE_SWAP(x, y);89 }90 if(x == 0)91 {92 //93 // Special handling for common case, from94 // Numerical Computation of Real or Complex Elliptic Integrals, eq.4795 //96 T xn = sqrt(y);97 T yn = sqrt(z);98 T x0 = xn;99 T y0 = yn;100 T sum = 0;101 T sum_pow = 0.25f;102 103 while(fabs(xn - yn) >= T(2.7) * tools::root_epsilon<T>() * fabs(xn))104 {105 T t = sqrt(xn * yn);106 xn = (xn + yn) / 2;107 yn = t;108 sum_pow *= 2;109 const auto temp = (xn - yn);110 sum += sum_pow * temp * temp;111 }112 T RF = constants::pi<T>() / (xn + yn);113 //114 // This following calculation suffers from serious cancellation when y ~ z115 // unless we combine terms. We have:116 //117 // ( ((x0 + y0)/2)^2 - z ) / (z(y-z))118 //119 // Substituting y = x0^2 and z = y0^2 and simplifying we get the following:120 //121 T pt = (x0 + 3 * y0) / (4 * z * (x0 + y0));122 //123 // Since we've moved the denominator from eq.47 inside the expression, we124 // need to also scale "sum" by the same value:125 //126 pt -= sum / (z * (y - z));127 return pt * RF * 3;128 }129 130 T xn = x;131 T yn = y;132 T zn = z;133 T An = (x + y + 3 * z) / 5;134 T A0 = An;135 // This has an extra 1.2 fudge factor which is really only needed when x, y and z are close in magnitude:136 T Q = pow(tools::epsilon<T>() / 4, -T(1) / 8) * BOOST_MATH_GPU_SAFE_MAX(BOOST_MATH_GPU_SAFE_MAX(An - x, An - y), An - z) * 1.2f;137 BOOST_MATH_INSTRUMENT_VARIABLE(Q);138 T lambda, rx, ry, rz;139 unsigned k = 0;140 T fn = 1;141 T RD_sum = 0;142 143 for(; k < policies::get_max_series_iterations<Policy>(); ++k)144 {145 rx = sqrt(xn);146 ry = sqrt(yn);147 rz = sqrt(zn);148 lambda = rx * ry + rx * rz + ry * rz;149 RD_sum += fn / (rz * (zn + lambda));150 An = (An + lambda) / 4;151 xn = (xn + lambda) / 4;152 yn = (yn + lambda) / 4;153 zn = (zn + lambda) / 4;154 fn /= 4;155 Q /= 4;156 BOOST_MATH_INSTRUMENT_VARIABLE(k);157 BOOST_MATH_INSTRUMENT_VARIABLE(RD_sum);158 BOOST_MATH_INSTRUMENT_VARIABLE(Q);159 if(Q < An)160 break;161 }162 163 policies::check_series_iterations<T, Policy>(function, k, pol);164 165 T X = fn * (A0 - x) / An;166 T Y = fn * (A0 - y) / An;167 T Z = -(X + Y) / 3;168 T E2 = X * Y - 6 * Z * Z;169 T E3 = (3 * X * Y - 8 * Z * Z) * Z;170 T E4 = 3 * (X * Y - Z * Z) * Z * Z;171 T E5 = X * Y * Z * Z * Z;172 173 T result = fn * pow(An, T(-3) / 2) *174 (1 - 3 * E2 / 14 + E3 / 6 + 9 * E2 * E2 / 88 - 3 * E4 / 22 - 9 * E2 * E3 / 52 + 3 * E5 / 26 - E2 * E2 * E2 / 16175 + 3 * E3 * E3 / 40 + 3 * E2 * E4 / 20 + 45 * E2 * E2 * E3 / 272 - 9 * (E3 * E4 + E2 * E5) / 68);176 BOOST_MATH_INSTRUMENT_VARIABLE(result);177 result += 3 * RD_sum;178 179 return result;180}181 182} // namespace detail183 184template <class T1, class T2, class T3, class Policy>185BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type 186 ellint_rd(T1 x, T2 y, T3 z, const Policy& pol)187{188 typedef typename tools::promote_args<T1, T2, T3>::type result_type;189 typedef typename policies::evaluation<result_type, Policy>::type value_type;190 return policies::checked_narrowing_cast<result_type, Policy>(191 detail::ellint_rd_imp(192 static_cast<value_type>(x),193 static_cast<value_type>(y),194 static_cast<value_type>(z), pol), "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)");195}196 197template <class T1, class T2, class T3>198BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type 199 ellint_rd(T1 x, T2 y, T3 z)200{201 return ellint_rd(x, y, z, policies::policy<>());202}203 204}} // namespaces205 206#endif // BOOST_MATH_ELLINT_RD_HPP207 208