1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006-2018 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 _GLIBCXX_BEGIN_NAMESPACE_VERSION 48 49 #if _GLIBCXX_USE_STD_SPEC_FUNCS 50 #elif defined(_GLIBCXX_TR1_CMATH) 51 namespace 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 __ellint_rf(_Tp __x,_Tp __y,_Tp __z)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 __comp_ellint_1_series(_Tp __k)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 __comp_ellint_1(_Tp __k)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 __ellint_1(_Tp __k,_Tp __phi)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 __comp_ellint_2_series(_Tp __k)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 __ellint_rd(_Tp __x,_Tp __y,_Tp __z)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 __comp_ellint_2(_Tp __k)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 __ellint_2(_Tp __k,_Tp __phi)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 __ellint_rc(_Tp __x,_Tp __y)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 __ellint_rj(_Tp __x,_Tp __y,_Tp __z,_Tp __p)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 __comp_ellint_3(_Tp __k,_Tp __nu)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 __ellint_3(_Tp __k,_Tp __nu,_Tp __phi)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