1[/ 2Copyright (c) 2006 Xiaogang Zhang 3Copyright (c) 2006 John Maddock 4Use, modification and distribution are subject to the 5Boost Software License, Version 1.0. (See accompanying file 6LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 7] 8 9[section:ellint_carlson Elliptic Integrals - Carlson Form] 10 11[heading Synopsis] 12 13`` 14 #include <boost/math/special_functions/ellint_rf.hpp> 15`` 16 17 namespace boost { namespace math { 18 19 template <class T1, class T2, class T3> 20 ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z) 21 22 template <class T1, class T2, class T3, class ``__Policy``> 23 ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z, const ``__Policy``&) 24 25 }} // namespaces 26 27 28`` 29 #include <boost/math/special_functions/ellint_rd.hpp> 30`` 31 32 namespace boost { namespace math { 33 34 template <class T1, class T2, class T3> 35 ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z) 36 37 template <class T1, class T2, class T3, class ``__Policy``> 38 ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z, const ``__Policy``&) 39 40 }} // namespaces 41 42 43`` 44 #include <boost/math/special_functions/ellint_rj.hpp> 45`` 46 47 namespace boost { namespace math { 48 49 template <class T1, class T2, class T3, class T4> 50 ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p) 51 52 template <class T1, class T2, class T3, class T4, class ``__Policy``> 53 ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p, const ``__Policy``&) 54 55 }} // namespaces 56 57 58`` 59 #include <boost/math/special_functions/ellint_rc.hpp> 60`` 61 62 namespace boost { namespace math { 63 64 template <class T1, class T2> 65 ``__sf_result`` ellint_rc(T1 x, T2 y) 66 67 template <class T1, class T2, class ``__Policy``> 68 ``__sf_result`` ellint_rc(T1 x, T2 y, const ``__Policy``&) 69 70 }} // namespaces 71 72`` 73 #include <boost/math/special_functions/ellint_rg.hpp> 74`` 75 76 namespace boost { namespace math { 77 78 template <class T1, class T2, class T3> 79 ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z) 80 81 template <class T1, class T2, class T3, class ``__Policy``> 82 ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z, const ``__Policy``&) 83 84 }} // namespaces 85 86 87 88[heading Description] 89 90These functions return Carlson's symmetrical elliptic integrals, the functions 91have complicated behavior over all their possible domains, but the following 92graph gives an idea of their behavior: 93 94[graph ellint_carlson] 95 96The return type of these functions is computed using the __arg_promotion_rules 97when the arguments are of different types: otherwise the return is the same type 98as the arguments. 99 100 template <class T1, class T2, class T3> 101 ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z) 102 103 template <class T1, class T2, class T3, class ``__Policy``> 104 ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z, const ``__Policy``&) 105 106Returns Carlson's Elliptic Integral ['R[sub F]]: 107 108[equation ellint9] 109 110Requires that all of the arguments are non-negative, and at most 111one may be zero. Otherwise returns the result of __domain_error. 112 113[optional_policy] 114 115 template <class T1, class T2, class T3> 116 ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z) 117 118 template <class T1, class T2, class T3, class ``__Policy``> 119 ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z, const ``__Policy``&) 120 121Returns Carlson's elliptic integral R[sub D]: 122 123[equation ellint10] 124 125Requires that x and y are non-negative, with at most one of them 126zero, and that z >= 0. Otherwise returns the result of __domain_error. 127 128[optional_policy] 129 130 template <class T1, class T2, class T3, class T4> 131 ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p) 132 133 template <class T1, class T2, class T3, class T4, class ``__Policy``> 134 ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p, const ``__Policy``&) 135 136Returns Carlson's elliptic integral R[sub J]: 137 138[equation ellint11] 139 140Requires that x, y and z are non-negative, with at most one of them 141zero, and that ['p != 0]. Otherwise returns the result of __domain_error. 142 143[optional_policy] 144 145When ['p < 0] the function returns the 146[@http://en.wikipedia.org/wiki/Cauchy_principal_value Cauchy principal value] 147using the relation: 148 149[equation ellint17] 150 151 template <class T1, class T2> 152 ``__sf_result`` ellint_rc(T1 x, T2 y) 153 154 template <class T1, class T2, class ``__Policy``> 155 ``__sf_result`` ellint_rc(T1 x, T2 y, const ``__Policy``&) 156 157Returns Carlson's elliptic integral R[sub C]: 158 159[equation ellint12] 160 161Requires that ['x > 0] and that ['y != 0]. 162Otherwise returns the result of __domain_error. 163 164[optional_policy] 165 166When ['y < 0] the function returns the 167[@http://mathworld.wolfram.com/CauchyPrincipalValue.html Cauchy principal value] 168using the relation: 169 170[equation ellint18] 171 172 template <class T1, class T2, class T3> 173 ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z) 174 175 template <class T1, class T2, class T3, class ``__Policy``> 176 ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z, const ``__Policy``&) 177 178Returns Carlson's elliptic integral ['R[sub G]:] 179 180[equation ellint27] 181 182Requires that x and y are non-negative, otherwise returns the result of __domain_error. 183 184[optional_policy] 185 186[heading Testing] 187 188There are two sets of tests. 189 190Spot tests compare selected values with test data given in: 191 192[:B. C. Carlson, ['[@http://arxiv.org/abs/math.CA/9409227 193Numerical computation of real or complex elliptic integrals]]. Numerical Algorithms, 194Volume 10, Number 1 / March, 1995, pp 13-26.] 195 196Random test data generated using NTL::RR at 1000-bit precision and our 197implementation checks for rounding-errors and/or regressions. 198 199There are also sanity checks that use the inter-relations between the integrals 200to verify their correctness: see the above Carlson paper for details. 201 202[heading Accuracy] 203 204These functions are computed using only basic arithmetic operations, so 205there isn't much variation in accuracy over differing platforms. 206Note that only results for the widest floating-point type on the 207system are given as narrower types have __zero_error. All values 208are relative errors in units of epsilon. 209 210[table_ellint_rc] 211 212[table_ellint_rd] 213 214[table_ellint_rg] 215 216[table_ellint_rf] 217 218[table_ellint_rj] 219 220 221[heading Implementation] 222 223The key of Carlson's algorithm [[link ellint_ref_carlson79 Carlson79]] is the 224duplication theorem: 225 226[equation ellint13] 227 228By applying it repeatedly, ['x], ['y], ['z] get 229closer and closer. When they are nearly equal, the special case equation 230 231[equation ellint16] 232 233is used. More specifically, ['[R F]] is evaluated from a Taylor series 234expansion to the fifth order. The calculations of the other three integrals 235are analogous, except for R[sub C] which can be computed from elementary functions. 236 237For ['p < 0] in ['R[sub J](x, y, z, p)] and ['y < 0] in ['R[sub C](x, y)], 238the integrals are singular and their 239[@http://mathworld.wolfram.com/CauchyPrincipalValue.html Cauchy principal values] 240are returned via the relations: 241 242[equation ellint17] 243 244[equation ellint18] 245 246[endsect] [/section:ellint_carlson Elliptic Integrals - Carlson Form] 247 248