]>
Commit | Line | Data |
---|---|---|
7c62b943 BK |
1 | // Special functions -*- C++ -*- |
2 | ||
cbe34bb5 | 3 | // Copyright (C) 2006-2017 Free Software Foundation, Inc. |
7c62b943 BK |
4 | // |
5 | // This file is part of the GNU ISO C++ Library. This library is free | |
6 | // software; you can redistribute it and/or modify it under the | |
7 | // terms of the GNU General Public License as published by the | |
748086b7 | 8 | // Free Software Foundation; either version 3, or (at your option) |
7c62b943 BK |
9 | // any later version. |
10 | // | |
11 | // This library is distributed in the hope that it will be useful, | |
12 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
14 | // GNU General Public License for more details. | |
15 | // | |
748086b7 JJ |
16 | // Under Section 7 of GPL version 3, you are granted additional |
17 | // permissions described in the GCC Runtime Library Exception, version | |
18 | // 3.1, as published by the Free Software Foundation. | |
19 | ||
20 | // You should have received a copy of the GNU General Public License and | |
21 | // a copy of the GCC Runtime Library Exception along with this program; | |
22 | // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see | |
23 | // <http://www.gnu.org/licenses/>. | |
7c62b943 BK |
24 | |
25 | /** @file tr1/ell_integral.tcc | |
26 | * This is an internal header file, included by other library headers. | |
f910786b | 27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
7c62b943 BK |
28 | */ |
29 | ||
30 | // | |
31 | // ISO C++ 14882 TR1: 5.2 Special functions | |
32 | // | |
33 | ||
34 | // Written by Edward Smith-Rowland based on: | |
35 | // (1) B. C. Carlson Numer. Math. 33, 1 (1979) | |
36 | // (2) B. C. Carlson, Special Functions of Applied Mathematics (1977) | |
37 | // (3) The Gnu Scientific Library, http://www.gnu.org/software/gsl | |
38 | // (4) Numerical Recipes in C, 2nd ed, by W. H. Press, S. A. Teukolsky, | |
39 | // W. T. Vetterling, B. P. Flannery, Cambridge University Press | |
40 | // (1992), pp. 261-269 | |
41 | ||
e133ace8 PC |
42 | #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC |
43 | #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1 | |
7c62b943 | 44 | |
12ffa228 | 45 | namespace std _GLIBCXX_VISIBILITY(default) |
7c62b943 | 46 | { |
4a15d842 FD |
47 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
48 | ||
f8571e51 | 49 | #if _GLIBCXX_USE_STD_SPEC_FUNCS |
2be75957 | 50 | #elif defined(_GLIBCXX_TR1_CMATH) |
e133ace8 PC |
51 | namespace tr1 |
52 | { | |
2be75957 ESR |
53 | #else |
54 | # error do not include this header directly, use <cmath> or <tr1/cmath> | |
55 | #endif | |
7c62b943 BK |
56 | // [5.2] Special functions |
57 | ||
7c62b943 | 58 | // Implementation-space details. |
7c62b943 BK |
59 | namespace __detail |
60 | { | |
7c62b943 BK |
61 | /** |
62 | * @brief Return the Carlson elliptic function @f$ R_F(x,y,z) @f$ | |
63 | * of the first kind. | |
64 | * | |
65 | * The Carlson elliptic function of the first kind is defined by: | |
66 | * @f[ | |
67 | * R_F(x,y,z) = \frac{1}{2} \int_0^\infty | |
68 | * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}} | |
69 | * @f] | |
70 | * | |
71 | * @param __x The first of three symmetric arguments. | |
72 | * @param __y The second of three symmetric arguments. | |
73 | * @param __z The third of three symmetric arguments. | |
74 | * @return The Carlson elliptic function of the first kind. | |
75 | */ | |
76 | template<typename _Tp> | |
77 | _Tp | |
be59c932 | 78 | __ellint_rf(_Tp __x, _Tp __y, _Tp __z) |
7c62b943 BK |
79 | { |
80 | const _Tp __min = std::numeric_limits<_Tp>::min(); | |
81 | const _Tp __max = std::numeric_limits<_Tp>::max(); | |
82 | const _Tp __lolim = _Tp(5) * __min; | |
83 | const _Tp __uplim = __max / _Tp(5); | |
84 | ||
85 | if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) | |
86 | std::__throw_domain_error(__N("Argument less than zero " | |
87 | "in __ellint_rf.")); | |
88 | else if (__x + __y < __lolim || __x + __z < __lolim | |
89 | || __y + __z < __lolim) | |
90 | std::__throw_domain_error(__N("Argument too small in __ellint_rf")); | |
91 | else | |
92 | { | |
93 | const _Tp __c0 = _Tp(1) / _Tp(4); | |
94 | const _Tp __c1 = _Tp(1) / _Tp(24); | |
95 | const _Tp __c2 = _Tp(1) / _Tp(10); | |
96 | const _Tp __c3 = _Tp(3) / _Tp(44); | |
97 | const _Tp __c4 = _Tp(1) / _Tp(14); | |
98 | ||
99 | _Tp __xn = __x; | |
100 | _Tp __yn = __y; | |
101 | _Tp __zn = __z; | |
102 | ||
103 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
104 | const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6)); | |
105 | _Tp __mu; | |
106 | _Tp __xndev, __yndev, __zndev; | |
107 | ||
108 | const unsigned int __max_iter = 100; | |
109 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) | |
110 | { | |
111 | __mu = (__xn + __yn + __zn) / _Tp(3); | |
112 | __xndev = 2 - (__mu + __xn) / __mu; | |
113 | __yndev = 2 - (__mu + __yn) / __mu; | |
114 | __zndev = 2 - (__mu + __zn) / __mu; | |
115 | _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); | |
116 | __epsilon = std::max(__epsilon, std::abs(__zndev)); | |
117 | if (__epsilon < __errtol) | |
118 | break; | |
119 | const _Tp __xnroot = std::sqrt(__xn); | |
120 | const _Tp __ynroot = std::sqrt(__yn); | |
121 | const _Tp __znroot = std::sqrt(__zn); | |
122 | const _Tp __lambda = __xnroot * (__ynroot + __znroot) | |
123 | + __ynroot * __znroot; | |
124 | __xn = __c0 * (__xn + __lambda); | |
125 | __yn = __c0 * (__yn + __lambda); | |
126 | __zn = __c0 * (__zn + __lambda); | |
127 | } | |
128 | ||
129 | const _Tp __e2 = __xndev * __yndev - __zndev * __zndev; | |
130 | const _Tp __e3 = __xndev * __yndev * __zndev; | |
131 | const _Tp __s = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2 | |
132 | + __c4 * __e3; | |
133 | ||
134 | return __s / std::sqrt(__mu); | |
135 | } | |
136 | } | |
137 | ||
138 | ||
139 | /** | |
140 | * @brief Return the complete elliptic integral of the first kind | |
141 | * @f$ K(k) @f$ by series expansion. | |
142 | * | |
143 | * The complete elliptic integral of the first kind is defined as | |
144 | * @f[ | |
145 | * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} | |
146 | * {\sqrt{1 - k^2sin^2\theta}} | |
147 | * @f] | |
148 | * | |
149 | * This routine is not bad as long as |k| is somewhat smaller than 1 | |
150 | * but is not is good as the Carlson elliptic integral formulation. | |
151 | * | |
152 | * @param __k The argument of the complete elliptic function. | |
153 | * @return The complete elliptic function of the first kind. | |
154 | */ | |
155 | template<typename _Tp> | |
156 | _Tp | |
be59c932 | 157 | __comp_ellint_1_series(_Tp __k) |
7c62b943 BK |
158 | { |
159 | ||
160 | const _Tp __kk = __k * __k; | |
161 | ||
162 | _Tp __term = __kk / _Tp(4); | |
163 | _Tp __sum = _Tp(1) + __term; | |
164 | ||
165 | const unsigned int __max_iter = 1000; | |
166 | for (unsigned int __i = 2; __i < __max_iter; ++__i) | |
167 | { | |
168 | __term *= (2 * __i - 1) * __kk / (2 * __i); | |
169 | if (__term < std::numeric_limits<_Tp>::epsilon()) | |
170 | break; | |
171 | __sum += __term; | |
172 | } | |
173 | ||
174 | return __numeric_constants<_Tp>::__pi_2() * __sum; | |
175 | } | |
176 | ||
177 | ||
178 | /** | |
179 | * @brief Return the complete elliptic integral of the first kind | |
180 | * @f$ K(k) @f$ using the Carlson formulation. | |
181 | * | |
182 | * The complete elliptic integral of the first kind is defined as | |
183 | * @f[ | |
184 | * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} | |
185 | * {\sqrt{1 - k^2 sin^2\theta}} | |
186 | * @f] | |
187 | * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the | |
188 | * first kind. | |
189 | * | |
190 | * @param __k The argument of the complete elliptic function. | |
191 | * @return The complete elliptic function of the first kind. | |
192 | */ | |
193 | template<typename _Tp> | |
194 | _Tp | |
be59c932 | 195 | __comp_ellint_1(_Tp __k) |
7c62b943 BK |
196 | { |
197 | ||
198 | if (__isnan(__k)) | |
199 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
200 | else if (std::abs(__k) >= _Tp(1)) | |
201 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
202 | else | |
203 | return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1)); | |
204 | } | |
205 | ||
206 | ||
207 | /** | |
208 | * @brief Return the incomplete elliptic integral of the first kind | |
209 | * @f$ F(k,\phi) @f$ using the Carlson formulation. | |
210 | * | |
211 | * The incomplete elliptic integral of the first kind is defined as | |
212 | * @f[ | |
213 | * F(k,\phi) = \int_0^{\phi}\frac{d\theta} | |
214 | * {\sqrt{1 - k^2 sin^2\theta}} | |
215 | * @f] | |
216 | * | |
217 | * @param __k The argument of the elliptic function. | |
218 | * @param __phi The integral limit argument of the elliptic function. | |
219 | * @return The elliptic function of the first kind. | |
220 | */ | |
221 | template<typename _Tp> | |
222 | _Tp | |
be59c932 | 223 | __ellint_1(_Tp __k, _Tp __phi) |
7c62b943 BK |
224 | { |
225 | ||
226 | if (__isnan(__k) || __isnan(__phi)) | |
227 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
228 | else if (std::abs(__k) > _Tp(1)) | |
229 | std::__throw_domain_error(__N("Bad argument in __ellint_1.")); | |
230 | else | |
231 | { | |
232 | // Reduce phi to -pi/2 < phi < +pi/2. | |
233 | const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() | |
234 | + _Tp(0.5L)); | |
235 | const _Tp __phi_red = __phi | |
236 | - __n * __numeric_constants<_Tp>::__pi(); | |
237 | ||
238 | const _Tp __s = std::sin(__phi_red); | |
239 | const _Tp __c = std::cos(__phi_red); | |
240 | ||
241 | const _Tp __F = __s | |
242 | * __ellint_rf(__c * __c, | |
243 | _Tp(1) - __k * __k * __s * __s, _Tp(1)); | |
244 | ||
245 | if (__n == 0) | |
246 | return __F; | |
247 | else | |
248 | return __F + _Tp(2) * __n * __comp_ellint_1(__k); | |
249 | } | |
250 | } | |
251 | ||
252 | ||
253 | /** | |
254 | * @brief Return the complete elliptic integral of the second kind | |
255 | * @f$ E(k) @f$ by series expansion. | |
256 | * | |
257 | * The complete elliptic integral of the second kind is defined as | |
258 | * @f[ | |
259 | * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} | |
260 | * @f] | |
261 | * | |
262 | * This routine is not bad as long as |k| is somewhat smaller than 1 | |
263 | * but is not is good as the Carlson elliptic integral formulation. | |
264 | * | |
265 | * @param __k The argument of the complete elliptic function. | |
266 | * @return The complete elliptic function of the second kind. | |
267 | */ | |
268 | template<typename _Tp> | |
269 | _Tp | |
be59c932 | 270 | __comp_ellint_2_series(_Tp __k) |
7c62b943 BK |
271 | { |
272 | ||
273 | const _Tp __kk = __k * __k; | |
274 | ||
275 | _Tp __term = __kk; | |
276 | _Tp __sum = __term; | |
277 | ||
278 | const unsigned int __max_iter = 1000; | |
279 | for (unsigned int __i = 2; __i < __max_iter; ++__i) | |
280 | { | |
281 | const _Tp __i2m = 2 * __i - 1; | |
282 | const _Tp __i2 = 2 * __i; | |
283 | __term *= __i2m * __i2m * __kk / (__i2 * __i2); | |
284 | if (__term < std::numeric_limits<_Tp>::epsilon()) | |
285 | break; | |
286 | __sum += __term / __i2m; | |
287 | } | |
288 | ||
289 | return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum); | |
290 | } | |
291 | ||
292 | ||
293 | /** | |
294 | * @brief Return the Carlson elliptic function of the second kind | |
295 | * @f$ R_D(x,y,z) = R_J(x,y,z,z) @f$ where | |
296 | * @f$ R_J(x,y,z,p) @f$ is the Carlson elliptic function | |
297 | * of the third kind. | |
298 | * | |
299 | * The Carlson elliptic function of the second kind is defined by: | |
300 | * @f[ | |
301 | * R_D(x,y,z) = \frac{3}{2} \int_0^\infty | |
302 | * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}} | |
303 | * @f] | |
304 | * | |
305 | * Based on Carlson's algorithms: | |
306 | * - B. C. Carlson Numer. Math. 33, 1 (1979) | |
307 | * - B. C. Carlson, Special Functions of Applied Mathematics (1977) | |
28dac70a | 308 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
7c62b943 BK |
309 | * by Press, Teukolsky, Vetterling, Flannery (1992) |
310 | * | |
311 | * @param __x The first of two symmetric arguments. | |
312 | * @param __y The second of two symmetric arguments. | |
313 | * @param __z The third argument. | |
314 | * @return The Carlson elliptic function of the second kind. | |
315 | */ | |
316 | template<typename _Tp> | |
317 | _Tp | |
be59c932 | 318 | __ellint_rd(_Tp __x, _Tp __y, _Tp __z) |
7c62b943 BK |
319 | { |
320 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
321 | const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); | |
322 | const _Tp __min = std::numeric_limits<_Tp>::min(); | |
323 | const _Tp __max = std::numeric_limits<_Tp>::max(); | |
324 | const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3)); | |
325 | const _Tp __uplim = std::pow(_Tp(0.1L) * __errtol / __min, _Tp(2) / _Tp(3)); | |
326 | ||
327 | if (__x < _Tp(0) || __y < _Tp(0)) | |
328 | std::__throw_domain_error(__N("Argument less than zero " | |
329 | "in __ellint_rd.")); | |
330 | else if (__x + __y < __lolim || __z < __lolim) | |
331 | std::__throw_domain_error(__N("Argument too small " | |
332 | "in __ellint_rd.")); | |
333 | else | |
334 | { | |
335 | const _Tp __c0 = _Tp(1) / _Tp(4); | |
336 | const _Tp __c1 = _Tp(3) / _Tp(14); | |
337 | const _Tp __c2 = _Tp(1) / _Tp(6); | |
338 | const _Tp __c3 = _Tp(9) / _Tp(22); | |
339 | const _Tp __c4 = _Tp(3) / _Tp(26); | |
340 | ||
341 | _Tp __xn = __x; | |
342 | _Tp __yn = __y; | |
343 | _Tp __zn = __z; | |
344 | _Tp __sigma = _Tp(0); | |
345 | _Tp __power4 = _Tp(1); | |
346 | ||
347 | _Tp __mu; | |
348 | _Tp __xndev, __yndev, __zndev; | |
349 | ||
350 | const unsigned int __max_iter = 100; | |
351 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) | |
352 | { | |
353 | __mu = (__xn + __yn + _Tp(3) * __zn) / _Tp(5); | |
354 | __xndev = (__mu - __xn) / __mu; | |
355 | __yndev = (__mu - __yn) / __mu; | |
356 | __zndev = (__mu - __zn) / __mu; | |
357 | _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); | |
358 | __epsilon = std::max(__epsilon, std::abs(__zndev)); | |
359 | if (__epsilon < __errtol) | |
360 | break; | |
361 | _Tp __xnroot = std::sqrt(__xn); | |
362 | _Tp __ynroot = std::sqrt(__yn); | |
363 | _Tp __znroot = std::sqrt(__zn); | |
364 | _Tp __lambda = __xnroot * (__ynroot + __znroot) | |
365 | + __ynroot * __znroot; | |
366 | __sigma += __power4 / (__znroot * (__zn + __lambda)); | |
367 | __power4 *= __c0; | |
368 | __xn = __c0 * (__xn + __lambda); | |
369 | __yn = __c0 * (__yn + __lambda); | |
370 | __zn = __c0 * (__zn + __lambda); | |
371 | } | |
372 | ||
22493a73 PC |
373 | // Note: __ea is an SPU badname. |
374 | _Tp __eaa = __xndev * __yndev; | |
7c62b943 | 375 | _Tp __eb = __zndev * __zndev; |
22493a73 PC |
376 | _Tp __ec = __eaa - __eb; |
377 | _Tp __ed = __eaa - _Tp(6) * __eb; | |
7c62b943 BK |
378 | _Tp __ef = __ed + __ec + __ec; |
379 | _Tp __s1 = __ed * (-__c1 + __c3 * __ed | |
380 | / _Tp(3) - _Tp(3) * __c4 * __zndev * __ef | |
381 | / _Tp(2)); | |
382 | _Tp __s2 = __zndev | |
383 | * (__c2 * __ef | |
22493a73 | 384 | + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa)); |
7c62b943 BK |
385 | |
386 | return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2) | |
387 | / (__mu * std::sqrt(__mu)); | |
388 | } | |
389 | } | |
390 | ||
391 | ||
392 | /** | |
393 | * @brief Return the complete elliptic integral of the second kind | |
394 | * @f$ E(k) @f$ using the Carlson formulation. | |
395 | * | |
396 | * The complete elliptic integral of the second kind is defined as | |
397 | * @f[ | |
398 | * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} | |
399 | * @f] | |
400 | * | |
401 | * @param __k The argument of the complete elliptic function. | |
402 | * @return The complete elliptic function of the second kind. | |
403 | */ | |
404 | template<typename _Tp> | |
405 | _Tp | |
be59c932 | 406 | __comp_ellint_2(_Tp __k) |
7c62b943 BK |
407 | { |
408 | ||
409 | if (__isnan(__k)) | |
410 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
411 | else if (std::abs(__k) == 1) | |
412 | return _Tp(1); | |
413 | else if (std::abs(__k) > _Tp(1)) | |
414 | std::__throw_domain_error(__N("Bad argument in __comp_ellint_2.")); | |
415 | else | |
416 | { | |
417 | const _Tp __kk = __k * __k; | |
418 | ||
419 | return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1)) | |
420 | - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3); | |
421 | } | |
422 | } | |
423 | ||
424 | ||
425 | /** | |
426 | * @brief Return the incomplete elliptic integral of the second kind | |
427 | * @f$ E(k,\phi) @f$ using the Carlson formulation. | |
428 | * | |
429 | * The incomplete elliptic integral of the second kind is defined as | |
430 | * @f[ | |
431 | * E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta} | |
432 | * @f] | |
433 | * | |
434 | * @param __k The argument of the elliptic function. | |
435 | * @param __phi The integral limit argument of the elliptic function. | |
436 | * @return The elliptic function of the second kind. | |
437 | */ | |
438 | template<typename _Tp> | |
439 | _Tp | |
be59c932 | 440 | __ellint_2(_Tp __k, _Tp __phi) |
7c62b943 BK |
441 | { |
442 | ||
443 | if (__isnan(__k) || __isnan(__phi)) | |
444 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
445 | else if (std::abs(__k) > _Tp(1)) | |
446 | std::__throw_domain_error(__N("Bad argument in __ellint_2.")); | |
447 | else | |
448 | { | |
449 | // Reduce phi to -pi/2 < phi < +pi/2. | |
450 | const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() | |
451 | + _Tp(0.5L)); | |
452 | const _Tp __phi_red = __phi | |
453 | - __n * __numeric_constants<_Tp>::__pi(); | |
454 | ||
455 | const _Tp __kk = __k * __k; | |
456 | const _Tp __s = std::sin(__phi_red); | |
457 | const _Tp __ss = __s * __s; | |
458 | const _Tp __sss = __ss * __s; | |
459 | const _Tp __c = std::cos(__phi_red); | |
460 | const _Tp __cc = __c * __c; | |
461 | ||
462 | const _Tp __E = __s | |
463 | * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1)) | |
464 | - __kk * __sss | |
465 | * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1)) | |
466 | / _Tp(3); | |
467 | ||
468 | if (__n == 0) | |
469 | return __E; | |
470 | else | |
471 | return __E + _Tp(2) * __n * __comp_ellint_2(__k); | |
472 | } | |
473 | } | |
474 | ||
475 | ||
476 | /** | |
477 | * @brief Return the Carlson elliptic function | |
478 | * @f$ R_C(x,y) = R_F(x,y,y) @f$ where @f$ R_F(x,y,z) @f$ | |
479 | * is the Carlson elliptic function of the first kind. | |
480 | * | |
481 | * The Carlson elliptic function is defined by: | |
482 | * @f[ | |
483 | * R_C(x,y) = \frac{1}{2} \int_0^\infty | |
484 | * \frac{dt}{(t + x)^{1/2}(t + y)} | |
485 | * @f] | |
486 | * | |
487 | * Based on Carlson's algorithms: | |
488 | * - B. C. Carlson Numer. Math. 33, 1 (1979) | |
489 | * - B. C. Carlson, Special Functions of Applied Mathematics (1977) | |
28dac70a | 490 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
7c62b943 BK |
491 | * by Press, Teukolsky, Vetterling, Flannery (1992) |
492 | * | |
493 | * @param __x The first argument. | |
494 | * @param __y The second argument. | |
495 | * @return The Carlson elliptic function. | |
496 | */ | |
497 | template<typename _Tp> | |
498 | _Tp | |
be59c932 | 499 | __ellint_rc(_Tp __x, _Tp __y) |
7c62b943 BK |
500 | { |
501 | const _Tp __min = std::numeric_limits<_Tp>::min(); | |
502 | const _Tp __max = std::numeric_limits<_Tp>::max(); | |
503 | const _Tp __lolim = _Tp(5) * __min; | |
504 | const _Tp __uplim = __max / _Tp(5); | |
505 | ||
506 | if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim) | |
507 | std::__throw_domain_error(__N("Argument less than zero " | |
508 | "in __ellint_rc.")); | |
509 | else | |
510 | { | |
511 | const _Tp __c0 = _Tp(1) / _Tp(4); | |
512 | const _Tp __c1 = _Tp(1) / _Tp(7); | |
513 | const _Tp __c2 = _Tp(9) / _Tp(22); | |
514 | const _Tp __c3 = _Tp(3) / _Tp(10); | |
515 | const _Tp __c4 = _Tp(3) / _Tp(8); | |
516 | ||
517 | _Tp __xn = __x; | |
518 | _Tp __yn = __y; | |
519 | ||
520 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
521 | const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6)); | |
522 | _Tp __mu; | |
523 | _Tp __sn; | |
524 | ||
525 | const unsigned int __max_iter = 100; | |
526 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) | |
527 | { | |
528 | __mu = (__xn + _Tp(2) * __yn) / _Tp(3); | |
529 | __sn = (__yn + __mu) / __mu - _Tp(2); | |
530 | if (std::abs(__sn) < __errtol) | |
531 | break; | |
532 | const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn) | |
533 | + __yn; | |
534 | __xn = __c0 * (__xn + __lambda); | |
535 | __yn = __c0 * (__yn + __lambda); | |
536 | } | |
537 | ||
538 | _Tp __s = __sn * __sn | |
539 | * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2))); | |
540 | ||
541 | return (_Tp(1) + __s) / std::sqrt(__mu); | |
542 | } | |
543 | } | |
544 | ||
545 | ||
546 | /** | |
547 | * @brief Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$ | |
548 | * of the third kind. | |
549 | * | |
550 | * The Carlson elliptic function of the third kind is defined by: | |
551 | * @f[ | |
552 | * R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty | |
553 | * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)} | |
554 | * @f] | |
555 | * | |
556 | * Based on Carlson's algorithms: | |
557 | * - B. C. Carlson Numer. Math. 33, 1 (1979) | |
558 | * - B. C. Carlson, Special Functions of Applied Mathematics (1977) | |
28dac70a | 559 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
7c62b943 BK |
560 | * by Press, Teukolsky, Vetterling, Flannery (1992) |
561 | * | |
562 | * @param __x The first of three symmetric arguments. | |
563 | * @param __y The second of three symmetric arguments. | |
564 | * @param __z The third of three symmetric arguments. | |
565 | * @param __p The fourth argument. | |
566 | * @return The Carlson elliptic function of the fourth kind. | |
567 | */ | |
568 | template<typename _Tp> | |
569 | _Tp | |
be59c932 | 570 | __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p) |
7c62b943 BK |
571 | { |
572 | const _Tp __min = std::numeric_limits<_Tp>::min(); | |
573 | const _Tp __max = std::numeric_limits<_Tp>::max(); | |
574 | const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3)); | |
575 | const _Tp __uplim = _Tp(0.3L) | |
576 | * std::pow(_Tp(0.2L) * __max, _Tp(1)/_Tp(3)); | |
577 | ||
578 | if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) | |
579 | std::__throw_domain_error(__N("Argument less than zero " | |
580 | "in __ellint_rj.")); | |
581 | else if (__x + __y < __lolim || __x + __z < __lolim | |
582 | || __y + __z < __lolim || __p < __lolim) | |
583 | std::__throw_domain_error(__N("Argument too small " | |
584 | "in __ellint_rj")); | |
585 | else | |
586 | { | |
587 | const _Tp __c0 = _Tp(1) / _Tp(4); | |
588 | const _Tp __c1 = _Tp(3) / _Tp(14); | |
589 | const _Tp __c2 = _Tp(1) / _Tp(3); | |
590 | const _Tp __c3 = _Tp(3) / _Tp(22); | |
591 | const _Tp __c4 = _Tp(3) / _Tp(26); | |
592 | ||
593 | _Tp __xn = __x; | |
594 | _Tp __yn = __y; | |
595 | _Tp __zn = __z; | |
596 | _Tp __pn = __p; | |
597 | _Tp __sigma = _Tp(0); | |
598 | _Tp __power4 = _Tp(1); | |
599 | ||
600 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
601 | const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); | |
602 | ||
603 | _Tp __lambda, __mu; | |
604 | _Tp __xndev, __yndev, __zndev, __pndev; | |
605 | ||
606 | const unsigned int __max_iter = 100; | |
607 | for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) | |
608 | { | |
609 | __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5); | |
610 | __xndev = (__mu - __xn) / __mu; | |
611 | __yndev = (__mu - __yn) / __mu; | |
612 | __zndev = (__mu - __zn) / __mu; | |
613 | __pndev = (__mu - __pn) / __mu; | |
614 | _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); | |
615 | __epsilon = std::max(__epsilon, std::abs(__zndev)); | |
616 | __epsilon = std::max(__epsilon, std::abs(__pndev)); | |
617 | if (__epsilon < __errtol) | |
618 | break; | |
619 | const _Tp __xnroot = std::sqrt(__xn); | |
620 | const _Tp __ynroot = std::sqrt(__yn); | |
621 | const _Tp __znroot = std::sqrt(__zn); | |
622 | const _Tp __lambda = __xnroot * (__ynroot + __znroot) | |
623 | + __ynroot * __znroot; | |
f070285a | 624 | const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot) |
7c62b943 | 625 | + __xnroot * __ynroot * __znroot; |
f070285a | 626 | const _Tp __alpha2 = __alpha1 * __alpha1; |
7c62b943 BK |
627 | const _Tp __beta = __pn * (__pn + __lambda) |
628 | * (__pn + __lambda); | |
629 | __sigma += __power4 * __ellint_rc(__alpha2, __beta); | |
630 | __power4 *= __c0; | |
631 | __xn = __c0 * (__xn + __lambda); | |
632 | __yn = __c0 * (__yn + __lambda); | |
633 | __zn = __c0 * (__zn + __lambda); | |
634 | __pn = __c0 * (__pn + __lambda); | |
635 | } | |
636 | ||
22493a73 PC |
637 | // Note: __ea is an SPU badname. |
638 | _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev; | |
7c62b943 BK |
639 | _Tp __eb = __xndev * __yndev * __zndev; |
640 | _Tp __ec = __pndev * __pndev; | |
22493a73 PC |
641 | _Tp __e2 = __eaa - _Tp(3) * __ec; |
642 | _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec); | |
7c62b943 BK |
643 | _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4) |
644 | - _Tp(3) * __c4 * __e3 / _Tp(2)); | |
645 | _Tp __s2 = __eb * (__c2 / _Tp(2) | |
646 | + __pndev * (-__c3 - __c3 + __pndev * __c4)); | |
22493a73 | 647 | _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3) |
7c62b943 BK |
648 | - __c2 * __pndev * __ec; |
649 | ||
650 | return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3) | |
651 | / (__mu * std::sqrt(__mu)); | |
652 | } | |
653 | } | |
654 | ||
655 | ||
656 | /** | |
657 | * @brief Return the complete elliptic integral of the third kind | |
658 | * @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ using the | |
659 | * Carlson formulation. | |
660 | * | |
661 | * The complete elliptic integral of the third kind is defined as | |
662 | * @f[ | |
663 | * \Pi(k,\nu) = \int_0^{\pi/2} | |
664 | * \frac{d\theta} | |
665 | * {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}} | |
666 | * @f] | |
667 | * | |
668 | * @param __k The argument of the elliptic function. | |
669 | * @param __nu The second argument of the elliptic function. | |
670 | * @return The complete elliptic function of the third kind. | |
671 | */ | |
672 | template<typename _Tp> | |
673 | _Tp | |
be59c932 | 674 | __comp_ellint_3(_Tp __k, _Tp __nu) |
7c62b943 BK |
675 | { |
676 | ||
677 | if (__isnan(__k) || __isnan(__nu)) | |
678 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
679 | else if (__nu == _Tp(1)) | |
680 | return std::numeric_limits<_Tp>::infinity(); | |
681 | else if (std::abs(__k) > _Tp(1)) | |
682 | std::__throw_domain_error(__N("Bad argument in __comp_ellint_3.")); | |
683 | else | |
684 | { | |
685 | const _Tp __kk = __k * __k; | |
686 | ||
687 | return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1)) | |
688 | - __nu | |
689 | * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu) | |
690 | / _Tp(3); | |
691 | } | |
692 | } | |
693 | ||
694 | ||
695 | /** | |
696 | * @brief Return the incomplete elliptic integral of the third kind | |
697 | * @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation. | |
698 | * | |
699 | * The incomplete elliptic integral of the third kind is defined as | |
700 | * @f[ | |
701 | * \Pi(k,\nu,\phi) = \int_0^{\phi} | |
702 | * \frac{d\theta} | |
703 | * {(1 - \nu \sin^2\theta) | |
704 | * \sqrt{1 - k^2 \sin^2\theta}} | |
705 | * @f] | |
706 | * | |
707 | * @param __k The argument of the elliptic function. | |
708 | * @param __nu The second argument of the elliptic function. | |
709 | * @param __phi The integral limit argument of the elliptic function. | |
710 | * @return The elliptic function of the third kind. | |
711 | */ | |
712 | template<typename _Tp> | |
713 | _Tp | |
be59c932 | 714 | __ellint_3(_Tp __k, _Tp __nu, _Tp __phi) |
7c62b943 BK |
715 | { |
716 | ||
717 | if (__isnan(__k) || __isnan(__nu) || __isnan(__phi)) | |
718 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
719 | else if (std::abs(__k) > _Tp(1)) | |
720 | std::__throw_domain_error(__N("Bad argument in __ellint_3.")); | |
721 | else | |
722 | { | |
723 | // Reduce phi to -pi/2 < phi < +pi/2. | |
724 | const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() | |
725 | + _Tp(0.5L)); | |
726 | const _Tp __phi_red = __phi | |
727 | - __n * __numeric_constants<_Tp>::__pi(); | |
728 | ||
729 | const _Tp __kk = __k * __k; | |
730 | const _Tp __s = std::sin(__phi_red); | |
731 | const _Tp __ss = __s * __s; | |
732 | const _Tp __sss = __ss * __s; | |
733 | const _Tp __c = std::cos(__phi_red); | |
734 | const _Tp __cc = __c * __c; | |
735 | ||
736 | const _Tp __Pi = __s | |
737 | * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1)) | |
738 | - __nu * __sss | |
739 | * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1), | |
740 | _Tp(1) + __nu * __ss) / _Tp(3); | |
741 | ||
742 | if (__n == 0) | |
743 | return __Pi; | |
744 | else | |
745 | return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu); | |
746 | } | |
747 | } | |
2be75957 | 748 | } // namespace __detail |
f8571e51 | 749 | #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) |
2be75957 ESR |
750 | } // namespace tr1 |
751 | #endif | |
4a15d842 FD |
752 | |
753 | _GLIBCXX_END_NAMESPACE_VERSION | |
7c62b943 BK |
754 | } |
755 | ||
e133ace8 | 756 | #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC |
7c62b943 | 757 |