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