1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006-2024 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 78 __ellint_rf(_Tp __x, _Tp __y, _Tp __z) 79 { 80 const _Tp __min = std::numeric_limits<_Tp>::min(); 81 const _Tp __lolim = _Tp(5) * __min; 82 83 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) 84 std::__throw_domain_error(__N("Argument less than zero " 85 "in __ellint_rf.")); 86 else if (__x + __y < __lolim || __x + __z < __lolim 87 || __y + __z < __lolim) 88 std::__throw_domain_error(__N("Argument too small in __ellint_rf")); 89 else 90 { 91 const _Tp __c0 = _Tp(1) / _Tp(4); 92 const _Tp __c1 = _Tp(1) / _Tp(24); 93 const _Tp __c2 = _Tp(1) / _Tp(10); 94 const _Tp __c3 = _Tp(3) / _Tp(44); 95 const _Tp __c4 = _Tp(1) / _Tp(14); 96 97 _Tp __xn = __x; 98 _Tp __yn = __y; 99 _Tp __zn = __z; 100 101 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 102 const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6)); 103 _Tp __mu; 104 _Tp __xndev, __yndev, __zndev; 105 106 const unsigned int __max_iter = 100; 107 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) 108 { 109 __mu = (__xn + __yn + __zn) / _Tp(3); 110 __xndev = 2 - (__mu + __xn) / __mu; 111 __yndev = 2 - (__mu + __yn) / __mu; 112 __zndev = 2 - (__mu + __zn) / __mu; 113 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); 114 __epsilon = std::max(__epsilon, std::abs(__zndev)); 115 if (__epsilon < __errtol) 116 break; 117 const _Tp __xnroot = std::sqrt(__xn); 118 const _Tp __ynroot = std::sqrt(__yn); 119 const _Tp __znroot = std::sqrt(__zn); 120 const _Tp __lambda = __xnroot * (__ynroot + __znroot) 121 + __ynroot * __znroot; 122 __xn = __c0 * (__xn + __lambda); 123 __yn = __c0 * (__yn + __lambda); 124 __zn = __c0 * (__zn + __lambda); 125 } 126 127 const _Tp __e2 = __xndev * __yndev - __zndev * __zndev; 128 const _Tp __e3 = __xndev * __yndev * __zndev; 129 const _Tp __s = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2 130 + __c4 * __e3; 131 132 return __s / std::sqrt(__mu); 133 } 134 } 135 136 137 /** 138 * @brief Return the complete elliptic integral of the first kind 139 * @f$ K(k) @f$ by series expansion. 140 * 141 * The complete elliptic integral of the first kind is defined as 142 * @f[ 143 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} 144 * {\sqrt{1 - k^2sin^2\theta}} 145 * @f] 146 * 147 * This routine is not bad as long as |k| is somewhat smaller than 1 148 * but is not is good as the Carlson elliptic integral formulation. 149 * 150 * @param __k The argument of the complete elliptic function. 151 * @return The complete elliptic function of the first kind. 152 */ 153 template<typename _Tp> 154 _Tp 155 __comp_ellint_1_series(_Tp __k) 156 { 157 158 const _Tp __kk = __k * __k; 159 160 _Tp __term = __kk / _Tp(4); 161 _Tp __sum = _Tp(1) + __term; 162 163 const unsigned int __max_iter = 1000; 164 for (unsigned int __i = 2; __i < __max_iter; ++__i) 165 { 166 __term *= (2 * __i - 1) * __kk / (2 * __i); 167 if (__term < std::numeric_limits<_Tp>::epsilon()) 168 break; 169 __sum += __term; 170 } 171 172 return __numeric_constants<_Tp>::__pi_2() * __sum; 173 } 174 175 176 /** 177 * @brief Return the complete elliptic integral of the first kind 178 * @f$ K(k) @f$ using the Carlson formulation. 179 * 180 * The complete elliptic integral of the first kind is defined as 181 * @f[ 182 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} 183 * {\sqrt{1 - k^2 sin^2\theta}} 184 * @f] 185 * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the 186 * first kind. 187 * 188 * @param __k The argument of the complete elliptic function. 189 * @return The complete elliptic function of the first kind. 190 */ 191 template<typename _Tp> 192 _Tp 193 __comp_ellint_1(_Tp __k) 194 { 195 196 if (__isnan(__k)) 197 return std::numeric_limits<_Tp>::quiet_NaN(); 198 else if (std::abs(__k) >= _Tp(1)) 199 return std::numeric_limits<_Tp>::quiet_NaN(); 200 else 201 return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1)); 202 } 203 204 205 /** 206 * @brief Return the incomplete elliptic integral of the first kind 207 * @f$ F(k,\phi) @f$ using the Carlson formulation. 208 * 209 * The incomplete elliptic integral of the first kind is defined as 210 * @f[ 211 * F(k,\phi) = \int_0^{\phi}\frac{d\theta} 212 * {\sqrt{1 - k^2 sin^2\theta}} 213 * @f] 214 * 215 * @param __k The argument of the elliptic function. 216 * @param __phi The integral limit argument of the elliptic function. 217 * @return The elliptic function of the first kind. 218 */ 219 template<typename _Tp> 220 _Tp 221 __ellint_1(_Tp __k, _Tp __phi) 222 { 223 224 if (__isnan(__k) || __isnan(__phi)) 225 return std::numeric_limits<_Tp>::quiet_NaN(); 226 else if (std::abs(__k) > _Tp(1)) 227 std::__throw_domain_error(__N("Bad argument in __ellint_1.")); 228 else 229 { 230 // Reduce phi to -pi/2 < phi < +pi/2. 231 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() 232 + _Tp(0.5L)); 233 const _Tp __phi_red = __phi 234 - __n * __numeric_constants<_Tp>::__pi(); 235 236 const _Tp __s = std::sin(__phi_red); 237 const _Tp __c = std::cos(__phi_red); 238 239 const _Tp __F = __s 240 * __ellint_rf(__c * __c, 241 _Tp(1) - __k * __k * __s * __s, _Tp(1)); 242 243 if (__n == 0) 244 return __F; 245 else 246 return __F + _Tp(2) * __n * __comp_ellint_1(__k); 247 } 248 } 249 250 251 /** 252 * @brief Return the complete elliptic integral of the second kind 253 * @f$ E(k) @f$ by series expansion. 254 * 255 * The complete elliptic integral of the second kind is defined as 256 * @f[ 257 * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} 258 * @f] 259 * 260 * This routine is not bad as long as |k| is somewhat smaller than 1 261 * but is not is good as the Carlson elliptic integral formulation. 262 * 263 * @param __k The argument of the complete elliptic function. 264 * @return The complete elliptic function of the second kind. 265 */ 266 template<typename _Tp> 267 _Tp 268 __comp_ellint_2_series(_Tp __k) 269 { 270 271 const _Tp __kk = __k * __k; 272 273 _Tp __term = __kk; 274 _Tp __sum = __term; 275 276 const unsigned int __max_iter = 1000; 277 for (unsigned int __i = 2; __i < __max_iter; ++__i) 278 { 279 const _Tp __i2m = 2 * __i - 1; 280 const _Tp __i2 = 2 * __i; 281 __term *= __i2m * __i2m * __kk / (__i2 * __i2); 282 if (__term < std::numeric_limits<_Tp>::epsilon()) 283 break; 284 __sum += __term / __i2m; 285 } 286 287 return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum); 288 } 289 290 291 /** 292 * @brief Return the Carlson elliptic function of the second kind 293 * @f$ R_D(x,y,z) = R_J(x,y,z,z) @f$ where 294 * @f$ R_J(x,y,z,p) @f$ is the Carlson elliptic function 295 * of the third kind. 296 * 297 * The Carlson elliptic function of the second kind is defined by: 298 * @f[ 299 * R_D(x,y,z) = \frac{3}{2} \int_0^\infty 300 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}} 301 * @f] 302 * 303 * Based on Carlson's algorithms: 304 * - B. C. Carlson Numer. Math. 33, 1 (1979) 305 * - B. C. Carlson, Special Functions of Applied Mathematics (1977) 306 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 307 * by Press, Teukolsky, Vetterling, Flannery (1992) 308 * 309 * @param __x The first of two symmetric arguments. 310 * @param __y The second of two symmetric arguments. 311 * @param __z The third argument. 312 * @return The Carlson elliptic function of the second kind. 313 */ 314 template<typename _Tp> 315 _Tp 316 __ellint_rd(_Tp __x, _Tp __y, _Tp __z) 317 { 318 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 319 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); 320 const _Tp __max = std::numeric_limits<_Tp>::max(); 321 const _Tp __lolim = _Tp(2) / std::pow(__max, _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 _Tp __ea = __xndev * __yndev; 370 _Tp __eb = __zndev * __zndev; 371 _Tp __ec = __ea - __eb; 372 _Tp __ed = __ea - _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 - __ea)); 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 __lolim = _Tp(5) * __min; 498 499 if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim) 500 std::__throw_domain_error(__N("Argument less than zero " 501 "in __ellint_rc.")); 502 else 503 { 504 const _Tp __c0 = _Tp(1) / _Tp(4); 505 const _Tp __c1 = _Tp(1) / _Tp(7); 506 const _Tp __c2 = _Tp(9) / _Tp(22); 507 const _Tp __c3 = _Tp(3) / _Tp(10); 508 const _Tp __c4 = _Tp(3) / _Tp(8); 509 510 _Tp __xn = __x; 511 _Tp __yn = __y; 512 513 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 514 const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6)); 515 _Tp __mu; 516 _Tp __sn; 517 518 const unsigned int __max_iter = 100; 519 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) 520 { 521 __mu = (__xn + _Tp(2) * __yn) / _Tp(3); 522 __sn = (__yn + __mu) / __mu - _Tp(2); 523 if (std::abs(__sn) < __errtol) 524 break; 525 const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn) 526 + __yn; 527 __xn = __c0 * (__xn + __lambda); 528 __yn = __c0 * (__yn + __lambda); 529 } 530 531 _Tp __s = __sn * __sn 532 * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2))); 533 534 return (_Tp(1) + __s) / std::sqrt(__mu); 535 } 536 } 537 538 539 /** 540 * @brief Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$ 541 * of the third kind. 542 * 543 * The Carlson elliptic function of the third kind is defined by: 544 * @f[ 545 * R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty 546 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)} 547 * @f] 548 * 549 * Based on Carlson's algorithms: 550 * - B. C. Carlson Numer. Math. 33, 1 (1979) 551 * - B. C. Carlson, Special Functions of Applied Mathematics (1977) 552 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 553 * by Press, Teukolsky, Vetterling, Flannery (1992) 554 * 555 * @param __x The first of three symmetric arguments. 556 * @param __y The second of three symmetric arguments. 557 * @param __z The third of three symmetric arguments. 558 * @param __p The fourth argument. 559 * @return The Carlson elliptic function of the fourth kind. 560 */ 561 template<typename _Tp> 562 _Tp 563 __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p) 564 { 565 const _Tp __min = std::numeric_limits<_Tp>::min(); 566 const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3)); 567 568 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) 569 std::__throw_domain_error(__N("Argument less than zero " 570 "in __ellint_rj.")); 571 else if (__x + __y < __lolim || __x + __z < __lolim 572 || __y + __z < __lolim || __p < __lolim) 573 std::__throw_domain_error(__N("Argument too small " 574 "in __ellint_rj")); 575 else 576 { 577 const _Tp __c0 = _Tp(1) / _Tp(4); 578 const _Tp __c1 = _Tp(3) / _Tp(14); 579 const _Tp __c2 = _Tp(1) / _Tp(3); 580 const _Tp __c3 = _Tp(3) / _Tp(22); 581 const _Tp __c4 = _Tp(3) / _Tp(26); 582 583 _Tp __xn = __x; 584 _Tp __yn = __y; 585 _Tp __zn = __z; 586 _Tp __pn = __p; 587 _Tp __sigma = _Tp(0); 588 _Tp __power4 = _Tp(1); 589 590 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 591 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); 592 593 _Tp __mu; 594 _Tp __xndev, __yndev, __zndev, __pndev; 595 596 const unsigned int __max_iter = 100; 597 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) 598 { 599 __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5); 600 __xndev = (__mu - __xn) / __mu; 601 __yndev = (__mu - __yn) / __mu; 602 __zndev = (__mu - __zn) / __mu; 603 __pndev = (__mu - __pn) / __mu; 604 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); 605 __epsilon = std::max(__epsilon, std::abs(__zndev)); 606 __epsilon = std::max(__epsilon, std::abs(__pndev)); 607 if (__epsilon < __errtol) 608 break; 609 const _Tp __xnroot = std::sqrt(__xn); 610 const _Tp __ynroot = std::sqrt(__yn); 611 const _Tp __znroot = std::sqrt(__zn); 612 const _Tp __lambda = __xnroot * (__ynroot + __znroot) 613 + __ynroot * __znroot; 614 const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot) 615 + __xnroot * __ynroot * __znroot; 616 const _Tp __alpha2 = __alpha1 * __alpha1; 617 const _Tp __beta = __pn * (__pn + __lambda) 618 * (__pn + __lambda); 619 __sigma += __power4 * __ellint_rc(__alpha2, __beta); 620 __power4 *= __c0; 621 __xn = __c0 * (__xn + __lambda); 622 __yn = __c0 * (__yn + __lambda); 623 __zn = __c0 * (__zn + __lambda); 624 __pn = __c0 * (__pn + __lambda); 625 } 626 627 _Tp __ea = __xndev * (__yndev + __zndev) + __yndev * __zndev; 628 _Tp __eb = __xndev * __yndev * __zndev; 629 _Tp __ec = __pndev * __pndev; 630 _Tp __e2 = __ea - _Tp(3) * __ec; 631 _Tp __e3 = __eb + _Tp(2) * __pndev * (__ea - __ec); 632 _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4) 633 - _Tp(3) * __c4 * __e3 / _Tp(2)); 634 _Tp __s2 = __eb * (__c2 / _Tp(2) 635 + __pndev * (-__c3 - __c3 + __pndev * __c4)); 636 _Tp __s3 = __pndev * __ea * (__c2 - __pndev * __c3) 637 - __c2 * __pndev * __ec; 638 639 return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3) 640 / (__mu * std::sqrt(__mu)); 641 } 642 } 643 644 645 /** 646 * @brief Return the complete elliptic integral of the third kind 647 * @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ using the 648 * Carlson formulation. 649 * 650 * The complete elliptic integral of the third kind is defined as 651 * @f[ 652 * \Pi(k,\nu) = \int_0^{\pi/2} 653 * \frac{d\theta} 654 * {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}} 655 * @f] 656 * 657 * @param __k The argument of the elliptic function. 658 * @param __nu The second argument of the elliptic function. 659 * @return The complete elliptic function of the third kind. 660 */ 661 template<typename _Tp> 662 _Tp 663 __comp_ellint_3(_Tp __k, _Tp __nu) 664 { 665 666 if (__isnan(__k) || __isnan(__nu)) 667 return std::numeric_limits<_Tp>::quiet_NaN(); 668 else if (__nu == _Tp(1)) 669 return std::numeric_limits<_Tp>::infinity(); 670 else if (std::abs(__k) > _Tp(1)) 671 std::__throw_domain_error(__N("Bad argument in __comp_ellint_3.")); 672 else 673 { 674 const _Tp __kk = __k * __k; 675 676 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1)) 677 + __nu 678 * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) - __nu) 679 / _Tp(3); 680 } 681 } 682 683 684 /** 685 * @brief Return the incomplete elliptic integral of the third kind 686 * @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation. 687 * 688 * The incomplete elliptic integral of the third kind is defined as 689 * @f[ 690 * \Pi(k,\nu,\phi) = \int_0^{\phi} 691 * \frac{d\theta} 692 * {(1 - \nu \sin^2\theta) 693 * \sqrt{1 - k^2 \sin^2\theta}} 694 * @f] 695 * 696 * @param __k The argument of the elliptic function. 697 * @param __nu The second argument of the elliptic function. 698 * @param __phi The integral limit argument of the elliptic function. 699 * @return The elliptic function of the third kind. 700 */ 701 template<typename _Tp> 702 _Tp 703 __ellint_3(_Tp __k, _Tp __nu, _Tp __phi) 704 { 705 706 if (__isnan(__k) || __isnan(__nu) || __isnan(__phi)) 707 return std::numeric_limits<_Tp>::quiet_NaN(); 708 else if (std::abs(__k) > _Tp(1)) 709 std::__throw_domain_error(__N("Bad argument in __ellint_3.")); 710 else 711 { 712 // Reduce phi to -pi/2 < phi < +pi/2. 713 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() 714 + _Tp(0.5L)); 715 const _Tp __phi_red = __phi 716 - __n * __numeric_constants<_Tp>::__pi(); 717 718 const _Tp __kk = __k * __k; 719 const _Tp __s = std::sin(__phi_red); 720 const _Tp __ss = __s * __s; 721 const _Tp __sss = __ss * __s; 722 const _Tp __c = std::cos(__phi_red); 723 const _Tp __cc = __c * __c; 724 725 const _Tp __Pi = __s 726 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1)) 727 + __nu * __sss 728 * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1), 729 _Tp(1) - __nu * __ss) / _Tp(3); 730 731 if (__n == 0) 732 return __Pi; 733 else 734 return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu); 735 } 736 } 737 } // namespace __detail 738 #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) 739 } // namespace tr1 740 #endif 741 742 _GLIBCXX_END_NAMESPACE_VERSION 743 } 744 745 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC 746 747
Indexes created Mon Mar 02 05:31:46 UTC 2026