]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | // | |
6 | // History: | |
7 | // XZ wrote the original of this file as part of the Google | |
8 | // Summer of Code 2006. JM modified it to fit into the | |
9 | // Boost.Math conceptual framework better, and to handle | |
10 | // types longer than 80-bit reals. | |
11 | // Updated 2015 to use Carlson's latest methods. | |
12 | // | |
13 | #ifndef BOOST_MATH_ELLINT_RF_HPP | |
14 | #define BOOST_MATH_ELLINT_RF_HPP | |
15 | ||
16 | #ifdef _MSC_VER | |
17 | #pragma once | |
18 | #endif | |
19 | ||
20 | #include <boost/math/special_functions/math_fwd.hpp> | |
21 | #include <boost/math/tools/config.hpp> | |
22 | #include <boost/math/constants/constants.hpp> | |
23 | #include <boost/math/policies/error_handling.hpp> | |
24 | #include <boost/math/special_functions/ellint_rc.hpp> | |
25 | ||
26 | // Carlson's elliptic integral of the first kind | |
27 | // R_F(x, y, z) = 0.5 * \int_{0}^{\infty} [(t+x)(t+y)(t+z)]^{-1/2} dt | |
28 | // Carlson, Numerische Mathematik, vol 33, 1 (1979) | |
29 | ||
30 | namespace boost { namespace math { namespace detail{ | |
31 | ||
32 | template <typename T, typename Policy> | |
33 | T ellint_rf_imp(T x, T y, T z, const Policy& pol) | |
34 | { | |
35 | BOOST_MATH_STD_USING | |
36 | using namespace boost::math; | |
37 | using std::swap; | |
38 | ||
39 | static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; | |
40 | ||
41 | if(x < 0 || y < 0 || z < 0) | |
42 | { | |
43 | return policies::raise_domain_error<T>(function, | |
44 | "domain error, all arguments must be non-negative, " | |
45 | "only sensible result is %1%.", | |
46 | std::numeric_limits<T>::quiet_NaN(), pol); | |
47 | } | |
48 | if(x + y == 0 || y + z == 0 || z + x == 0) | |
49 | { | |
50 | return policies::raise_domain_error<T>(function, | |
51 | "domain error, at most one argument can be zero, " | |
52 | "only sensible result is %1%.", | |
53 | std::numeric_limits<T>::quiet_NaN(), pol); | |
54 | } | |
55 | // | |
56 | // Special cases from http://dlmf.nist.gov/19.20#i | |
57 | // | |
58 | if(x == y) | |
59 | { | |
60 | if(x == z) | |
61 | { | |
62 | // x, y, z equal: | |
63 | return 1 / sqrt(x); | |
64 | } | |
65 | else | |
66 | { | |
67 | // 2 equal, x and y: | |
68 | if(z == 0) | |
69 | return constants::pi<T>() / (2 * sqrt(x)); | |
70 | else | |
71 | return ellint_rc_imp(z, x, pol); | |
72 | } | |
73 | } | |
74 | if(x == z) | |
75 | { | |
76 | if(y == 0) | |
77 | return constants::pi<T>() / (2 * sqrt(x)); | |
78 | else | |
79 | return ellint_rc_imp(y, x, pol); | |
80 | } | |
81 | if(y == z) | |
82 | { | |
83 | if(x == 0) | |
84 | return constants::pi<T>() / (2 * sqrt(y)); | |
85 | else | |
86 | return ellint_rc_imp(x, y, pol); | |
87 | } | |
88 | if(x == 0) | |
89 | swap(x, z); | |
90 | else if(y == 0) | |
91 | swap(y, z); | |
92 | if(z == 0) | |
93 | { | |
94 | // | |
95 | // Special case for one value zero: | |
96 | // | |
97 | T xn = sqrt(x); | |
98 | T yn = sqrt(y); | |
99 | ||
100 | while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn)) | |
101 | { | |
102 | T t = sqrt(xn * yn); | |
103 | xn = (xn + yn) / 2; | |
104 | yn = t; | |
105 | } | |
106 | return constants::pi<T>() / (xn + yn); | |
107 | } | |
108 | ||
109 | T xn = x; | |
110 | T yn = y; | |
111 | T zn = z; | |
112 | T An = (x + y + z) / 3; | |
113 | T A0 = An; | |
114 | T Q = pow(3 * boost::math::tools::epsilon<T>(), T(-1) / 8) * (std::max)((std::max)(fabs(An - xn), fabs(An - yn)), fabs(An - zn)); | |
115 | T fn = 1; | |
116 | ||
117 | ||
118 | // duplication | |
119 | unsigned k = 1; | |
120 | for(; k < boost::math::policies::get_max_series_iterations<Policy>(); ++k) | |
121 | { | |
122 | T root_x = sqrt(xn); | |
123 | T root_y = sqrt(yn); | |
124 | T root_z = sqrt(zn); | |
125 | T lambda = root_x * root_y + root_x * root_z + root_y * root_z; | |
126 | An = (An + lambda) / 4; | |
127 | xn = (xn + lambda) / 4; | |
128 | yn = (yn + lambda) / 4; | |
129 | zn = (zn + lambda) / 4; | |
130 | Q /= 4; | |
131 | fn *= 4; | |
132 | if(Q < fabs(An)) | |
133 | break; | |
134 | } | |
135 | // Check to see if we gave up too soon: | |
136 | policies::check_series_iterations<T>(function, k, pol); | |
137 | BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
138 | ||
139 | T X = (A0 - x) / (An * fn); | |
140 | T Y = (A0 - y) / (An * fn); | |
141 | T Z = -X - Y; | |
142 | ||
143 | // Taylor series expansion to the 7th order | |
144 | T E2 = X * Y - Z * Z; | |
145 | T E3 = X * Y * Z; | |
146 | 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); | |
147 | } | |
148 | ||
149 | } // namespace detail | |
150 | ||
151 | template <class T1, class T2, class T3, class Policy> | |
152 | inline typename tools::promote_args<T1, T2, T3>::type | |
153 | ellint_rf(T1 x, T2 y, T3 z, const Policy& pol) | |
154 | { | |
155 | typedef typename tools::promote_args<T1, T2, T3>::type result_type; | |
156 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
157 | return policies::checked_narrowing_cast<result_type, Policy>( | |
158 | detail::ellint_rf_imp( | |
159 | static_cast<value_type>(x), | |
160 | static_cast<value_type>(y), | |
161 | static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"); | |
162 | } | |
163 | ||
164 | template <class T1, class T2, class T3> | |
165 | inline typename tools::promote_args<T1, T2, T3>::type | |
166 | ellint_rf(T1 x, T2 y, T3 z) | |
167 | { | |
168 | return ellint_rf(x, y, z, policies::policy<>()); | |
169 | } | |
170 | ||
171 | }} // namespaces | |
172 | ||
173 | #endif // BOOST_MATH_ELLINT_RF_HPP | |
174 |