170 lines · plain
1// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock2// 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 to fit into the10// Boost.Math conceptual framework better, and to handle11// types longer than 80-bit reals.12// Updated 2015 to use Carlson's latest methods.13//14#ifndef BOOST_MATH_ELLINT_RF_HPP15#define BOOST_MATH_ELLINT_RF_HPP16 17#ifdef _MSC_VER18#pragma once19#endif20 21#include <boost/math/tools/config.hpp>22#include <boost/math/tools/numeric_limits.hpp>23#include <boost/math/special_functions/math_fwd.hpp>24#include <boost/math/constants/constants.hpp>25#include <boost/math/policies/error_handling.hpp>26#include <boost/math/special_functions/ellint_rc.hpp>27 28// Carlson's elliptic integral of the first kind29// R_F(x, y, z) = 0.5 * \int_{0}^{\infty} [(t+x)(t+y)(t+z)]^{-1/2} dt30// Carlson, Numerische Mathematik, vol 33, 1 (1979)31 32namespace boost { namespace math { namespace detail{33 34 template <typename T, typename Policy>35 BOOST_MATH_GPU_ENABLED T ellint_rf_imp(T x, T y, T z, const Policy& pol)36 {37 BOOST_MATH_STD_USING38 using namespace boost::math;39 40 constexpr auto function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)";41 42 if(x < 0 || y < 0 || z < 0)43 {44 return policies::raise_domain_error<T>(function, "domain error, all arguments must be non-negative, only sensible result is %1%.", boost::math::numeric_limits<T>::quiet_NaN(), pol);45 }46 if(x + y == 0 || y + z == 0 || z + x == 0)47 {48 return policies::raise_domain_error<T>(function, "domain error, at most one argument can be zero, only sensible result is %1%.", boost::math::numeric_limits<T>::quiet_NaN(), pol);49 }50 //51 // Special cases from http://dlmf.nist.gov/19.20#i52 //53 if(x == y)54 {55 if(x == z)56 {57 // x, y, z equal:58 return 1 / sqrt(x);59 }60 else61 {62 // 2 equal, x and y:63 if(z == 0)64 return constants::pi<T>() / (2 * sqrt(x));65 else66 return ellint_rc_imp(z, x, pol);67 }68 }69 if(x == z)70 {71 if(y == 0)72 return constants::pi<T>() / (2 * sqrt(x));73 else74 return ellint_rc_imp(y, x, pol);75 }76 if(y == z)77 {78 if(x == 0)79 return constants::pi<T>() / (2 * sqrt(y));80 else81 return ellint_rc_imp(x, y, pol);82 }83 if(x == 0)84 BOOST_MATH_GPU_SAFE_SWAP(x, z);85 else if(y == 0)86 BOOST_MATH_GPU_SAFE_SWAP(y, z);87 if(z == 0)88 {89 //90 // Special case for one value zero:91 //92 T xn = sqrt(x);93 T yn = sqrt(y);94 95 while(fabs(xn - yn) >= T(2.7) * tools::root_epsilon<T>() * fabs(xn))96 {97 T t = sqrt(xn * yn);98 xn = (xn + yn) / 2;99 yn = t;100 }101 return constants::pi<T>() / (xn + yn);102 }103 104 T xn = x;105 T yn = y;106 T zn = z;107 T An = (x + y + z) / 3;108 T A0 = An;109 T Q = pow(3 * boost::math::tools::epsilon<T>(), T(-1) / 8) * BOOST_MATH_GPU_SAFE_MAX(BOOST_MATH_GPU_SAFE_MAX(fabs(An - xn), fabs(An - yn)), fabs(An - zn));110 T fn = 1;111 112 113 // duplication114 unsigned k = 1;115 for(; k < boost::math::policies::get_max_series_iterations<Policy>(); ++k)116 {117 T root_x = sqrt(xn);118 T root_y = sqrt(yn);119 T root_z = sqrt(zn);120 T lambda = root_x * root_y + root_x * root_z + root_y * root_z;121 An = (An + lambda) / 4;122 xn = (xn + lambda) / 4;123 yn = (yn + lambda) / 4;124 zn = (zn + lambda) / 4;125 Q /= 4;126 fn *= 4;127 if(Q < fabs(An))128 break;129 }130 // Check to see if we gave up too soon:131 policies::check_series_iterations<T>(function, k, pol);132 BOOST_MATH_INSTRUMENT_VARIABLE(k);133 134 T X = (A0 - x) / (An * fn);135 T Y = (A0 - y) / (An * fn);136 T Z = -X - Y;137 138 // Taylor series expansion to the 7th order139 T E2 = X * Y - Z * Z;140 T E3 = X * Y * Z;141 return (1 + E3 * (T(1) / 14 + 3 * E3 / 104) + E2 * (T(-1) / 10 + E2 / 24 - (3 * E3) / 44 - 5 * E2 * E2 / 208 + E2 * E3 / 16)) / sqrt(An);142 }143 144} // namespace detail145 146template <class T1, class T2, class T3, class Policy>147BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type 148 ellint_rf(T1 x, T2 y, T3 z, const Policy& pol)149{150 typedef typename tools::promote_args<T1, T2, T3>::type result_type;151 typedef typename policies::evaluation<result_type, Policy>::type value_type;152 return policies::checked_narrowing_cast<result_type, Policy>(153 detail::ellint_rf_imp(154 static_cast<value_type>(x),155 static_cast<value_type>(y),156 static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)");157}158 159template <class T1, class T2, class T3>160BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2, T3>::type 161 ellint_rf(T1 x, T2 y, T3 z)162{163 return ellint_rf(x, y, z, policies::policy<>());164}165 166}} // namespaces167 168#endif // BOOST_MATH_ELLINT_RF_HPP169 170