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