1// Special functions -*- C++ -*- 
2 
3// Copyright (C) 2006-2019 Free Software Foundation, Inc. 
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 
8// Free Software Foundation; either version 3, or (at your option) 
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// 
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/>. 
24 
25/** @file tr1/ell_integral.tcc 
26 * This is an internal header file, included by other library headers. 
27 * Do not attempt to use it directly. @headername{tr1/cmath} 
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 
42#ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC 
43#define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1 
44 
45namespace std _GLIBCXX_VISIBILITY(default) 
46{ 
47_GLIBCXX_BEGIN_NAMESPACE_VERSION 
48 
49#if _GLIBCXX_USE_STD_SPEC_FUNCS 
50#elif defined(_GLIBCXX_TR1_CMATH) 
51namespace tr1 
52{ 
53#else 
54# error do not include this header directly, use <cmath> or <tr1/cmath> 
55#endif 
56 // [5.2] Special functions 
57 
58 // Implementation-space details. 
59 namespace __detail 
60 { 
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 
78 __ellint_rf(_Tp __x, _Tp __y, _Tp __z) 
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 
157 __comp_ellint_1_series(_Tp __k) 
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 
195 __comp_ellint_1(_Tp __k) 
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 
223 __ellint_1(_Tp __k, _Tp __phi) 
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 
270 __comp_ellint_2_series(_Tp __k) 
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) 
308 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 
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 
318 __ellint_rd(_Tp __x, _Tp __y, _Tp __z) 
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 
373 // Note: __ea is an SPU badname. 
374 _Tp __eaa = __xndev * __yndev; 
375 _Tp __eb = __zndev * __zndev; 
376 _Tp __ec = __eaa - __eb; 
377 _Tp __ed = __eaa - _Tp(6) * __eb; 
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 
384 + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa)); 
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 
406 __comp_ellint_2(_Tp __k) 
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 
440 __ellint_2(_Tp __k, _Tp __phi) 
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) 
490 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 
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 
499 __ellint_rc(_Tp __x, _Tp __y) 
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) 
559 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 
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 
570 __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p) 
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; 
624 const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot) 
625 + __xnroot * __ynroot * __znroot; 
626 const _Tp __alpha2 = __alpha1 * __alpha1; 
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 
637 // Note: __ea is an SPU badname. 
638 _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev; 
639 _Tp __eb = __xndev * __yndev * __zndev; 
640 _Tp __ec = __pndev * __pndev; 
641 _Tp __e2 = __eaa - _Tp(3) * __ec; 
642 _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec); 
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)); 
647 _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3) 
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 
674 __comp_ellint_3(_Tp __k, _Tp __nu) 
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 
714 __ellint_3(_Tp __k, _Tp __nu, _Tp __phi) 
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 } 
748 } // namespace __detail 
749#if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) 
750} // namespace tr1 
751#endif 
752 
753_GLIBCXX_END_NAMESPACE_VERSION 
754} 
755 
756#endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC 
757 
758