libstdc++
|
00001 // Special functions -*- C++ -*- 00002 00003 // Copyright (C) 2006, 2007, 2008, 2009 00004 // Free Software Foundation, Inc. 00005 // 00006 // This file is part of the GNU ISO C++ Library. This library is free 00007 // software; you can redistribute it and/or modify it under the 00008 // terms of the GNU General Public License as published by the 00009 // Free Software Foundation; either version 3, or (at your option) 00010 // any later version. 00011 // 00012 // This library is distributed in the hope that it will be useful, 00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 // GNU General Public License for more details. 00016 // 00017 // Under Section 7 of GPL version 3, you are granted additional 00018 // permissions described in the GCC Runtime Library Exception, version 00019 // 3.1, as published by the Free Software Foundation. 00020 00021 // You should have received a copy of the GNU General Public License and 00022 // a copy of the GCC Runtime Library Exception along with this program; 00023 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 00024 // <http://www.gnu.org/licenses/>. 00025 00026 /** @file tr1/ell_integral.tcc 00027 * This is an internal header file, included by other library headers. 00028 * You should not attempt to use it directly. 00029 */ 00030 00031 // 00032 // ISO C++ 14882 TR1: 5.2 Special functions 00033 // 00034 00035 // Written by Edward Smith-Rowland based on: 00036 // (1) B. C. Carlson Numer. Math. 33, 1 (1979) 00037 // (2) B. C. Carlson, Special Functions of Applied Mathematics (1977) 00038 // (3) The Gnu Scientific Library, http://www.gnu.org/software/gsl 00039 // (4) Numerical Recipes in C, 2nd ed, by W. H. Press, S. A. Teukolsky, 00040 // W. T. Vetterling, B. P. Flannery, Cambridge University Press 00041 // (1992), pp. 261-269 00042 00043 #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC 00044 #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1 00045 00046 namespace std 00047 { 00048 namespace tr1 00049 { 00050 00051 // [5.2] Special functions 00052 00053 // Implementation-space details. 00054 namespace __detail 00055 { 00056 00057 /** 00058 * @brief Return the Carlson elliptic function @f$ R_F(x,y,z) @f$ 00059 * of the first kind. 00060 * 00061 * The Carlson elliptic function of the first kind is defined by: 00062 * @f[ 00063 * R_F(x,y,z) = \frac{1}{2} \int_0^\infty 00064 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}} 00065 * @f] 00066 * 00067 * @param __x The first of three symmetric arguments. 00068 * @param __y The second of three symmetric arguments. 00069 * @param __z The third of three symmetric arguments. 00070 * @return The Carlson elliptic function of the first kind. 00071 */ 00072 template<typename _Tp> 00073 _Tp 00074 __ellint_rf(const _Tp __x, const _Tp __y, const _Tp __z) 00075 { 00076 const _Tp __min = std::numeric_limits<_Tp>::min(); 00077 const _Tp __max = std::numeric_limits<_Tp>::max(); 00078 const _Tp __lolim = _Tp(5) * __min; 00079 const _Tp __uplim = __max / _Tp(5); 00080 00081 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) 00082 std::__throw_domain_error(__N("Argument less than zero " 00083 "in __ellint_rf.")); 00084 else if (__x + __y < __lolim || __x + __z < __lolim 00085 || __y + __z < __lolim) 00086 std::__throw_domain_error(__N("Argument too small in __ellint_rf")); 00087 else 00088 { 00089 const _Tp __c0 = _Tp(1) / _Tp(4); 00090 const _Tp __c1 = _Tp(1) / _Tp(24); 00091 const _Tp __c2 = _Tp(1) / _Tp(10); 00092 const _Tp __c3 = _Tp(3) / _Tp(44); 00093 const _Tp __c4 = _Tp(1) / _Tp(14); 00094 00095 _Tp __xn = __x; 00096 _Tp __yn = __y; 00097 _Tp __zn = __z; 00098 00099 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 00100 const _Tp __errtol = std::pow(__eps, _Tp(1) / _Tp(6)); 00101 _Tp __mu; 00102 _Tp __xndev, __yndev, __zndev; 00103 00104 const unsigned int __max_iter = 100; 00105 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) 00106 { 00107 __mu = (__xn + __yn + __zn) / _Tp(3); 00108 __xndev = 2 - (__mu + __xn) / __mu; 00109 __yndev = 2 - (__mu + __yn) / __mu; 00110 __zndev = 2 - (__mu + __zn) / __mu; 00111 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); 00112 __epsilon = std::max(__epsilon, std::abs(__zndev)); 00113 if (__epsilon < __errtol) 00114 break; 00115 const _Tp __xnroot = std::sqrt(__xn); 00116 const _Tp __ynroot = std::sqrt(__yn); 00117 const _Tp __znroot = std::sqrt(__zn); 00118 const _Tp __lambda = __xnroot * (__ynroot + __znroot) 00119 + __ynroot * __znroot; 00120 __xn = __c0 * (__xn + __lambda); 00121 __yn = __c0 * (__yn + __lambda); 00122 __zn = __c0 * (__zn + __lambda); 00123 } 00124 00125 const _Tp __e2 = __xndev * __yndev - __zndev * __zndev; 00126 const _Tp __e3 = __xndev * __yndev * __zndev; 00127 const _Tp __s = _Tp(1) + (__c1 * __e2 - __c2 - __c3 * __e3) * __e2 00128 + __c4 * __e3; 00129 00130 return __s / std::sqrt(__mu); 00131 } 00132 } 00133 00134 00135 /** 00136 * @brief Return the complete elliptic integral of the first kind 00137 * @f$ K(k) @f$ by series expansion. 00138 * 00139 * The complete elliptic integral of the first kind is defined as 00140 * @f[ 00141 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} 00142 * {\sqrt{1 - k^2sin^2\theta}} 00143 * @f] 00144 * 00145 * This routine is not bad as long as |k| is somewhat smaller than 1 00146 * but is not is good as the Carlson elliptic integral formulation. 00147 * 00148 * @param __k The argument of the complete elliptic function. 00149 * @return The complete elliptic function of the first kind. 00150 */ 00151 template<typename _Tp> 00152 _Tp 00153 __comp_ellint_1_series(const _Tp __k) 00154 { 00155 00156 const _Tp __kk = __k * __k; 00157 00158 _Tp __term = __kk / _Tp(4); 00159 _Tp __sum = _Tp(1) + __term; 00160 00161 const unsigned int __max_iter = 1000; 00162 for (unsigned int __i = 2; __i < __max_iter; ++__i) 00163 { 00164 __term *= (2 * __i - 1) * __kk / (2 * __i); 00165 if (__term < std::numeric_limits<_Tp>::epsilon()) 00166 break; 00167 __sum += __term; 00168 } 00169 00170 return __numeric_constants<_Tp>::__pi_2() * __sum; 00171 } 00172 00173 00174 /** 00175 * @brief Return the complete elliptic integral of the first kind 00176 * @f$ K(k) @f$ using the Carlson formulation. 00177 * 00178 * The complete elliptic integral of the first kind is defined as 00179 * @f[ 00180 * K(k) = F(k,\pi/2) = \int_0^{\pi/2}\frac{d\theta} 00181 * {\sqrt{1 - k^2 sin^2\theta}} 00182 * @f] 00183 * where @f$ F(k,\phi) @f$ is the incomplete elliptic integral of the 00184 * first kind. 00185 * 00186 * @param __k The argument of the complete elliptic function. 00187 * @return The complete elliptic function of the first kind. 00188 */ 00189 template<typename _Tp> 00190 _Tp 00191 __comp_ellint_1(const _Tp __k) 00192 { 00193 00194 if (__isnan(__k)) 00195 return std::numeric_limits<_Tp>::quiet_NaN(); 00196 else if (std::abs(__k) >= _Tp(1)) 00197 return std::numeric_limits<_Tp>::quiet_NaN(); 00198 else 00199 return __ellint_rf(_Tp(0), _Tp(1) - __k * __k, _Tp(1)); 00200 } 00201 00202 00203 /** 00204 * @brief Return the incomplete elliptic integral of the first kind 00205 * @f$ F(k,\phi) @f$ using the Carlson formulation. 00206 * 00207 * The incomplete elliptic integral of the first kind is defined as 00208 * @f[ 00209 * F(k,\phi) = \int_0^{\phi}\frac{d\theta} 00210 * {\sqrt{1 - k^2 sin^2\theta}} 00211 * @f] 00212 * 00213 * @param __k The argument of the elliptic function. 00214 * @param __phi The integral limit argument of the elliptic function. 00215 * @return The elliptic function of the first kind. 00216 */ 00217 template<typename _Tp> 00218 _Tp 00219 __ellint_1(const _Tp __k, const _Tp __phi) 00220 { 00221 00222 if (__isnan(__k) || __isnan(__phi)) 00223 return std::numeric_limits<_Tp>::quiet_NaN(); 00224 else if (std::abs(__k) > _Tp(1)) 00225 std::__throw_domain_error(__N("Bad argument in __ellint_1.")); 00226 else 00227 { 00228 // Reduce phi to -pi/2 < phi < +pi/2. 00229 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() 00230 + _Tp(0.5L)); 00231 const _Tp __phi_red = __phi 00232 - __n * __numeric_constants<_Tp>::__pi(); 00233 00234 const _Tp __s = std::sin(__phi_red); 00235 const _Tp __c = std::cos(__phi_red); 00236 00237 const _Tp __F = __s 00238 * __ellint_rf(__c * __c, 00239 _Tp(1) - __k * __k * __s * __s, _Tp(1)); 00240 00241 if (__n == 0) 00242 return __F; 00243 else 00244 return __F + _Tp(2) * __n * __comp_ellint_1(__k); 00245 } 00246 } 00247 00248 00249 /** 00250 * @brief Return the complete elliptic integral of the second kind 00251 * @f$ E(k) @f$ by series expansion. 00252 * 00253 * The complete elliptic integral of the second kind is defined as 00254 * @f[ 00255 * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} 00256 * @f] 00257 * 00258 * This routine is not bad as long as |k| is somewhat smaller than 1 00259 * but is not is good as the Carlson elliptic integral formulation. 00260 * 00261 * @param __k The argument of the complete elliptic function. 00262 * @return The complete elliptic function of the second kind. 00263 */ 00264 template<typename _Tp> 00265 _Tp 00266 __comp_ellint_2_series(const _Tp __k) 00267 { 00268 00269 const _Tp __kk = __k * __k; 00270 00271 _Tp __term = __kk; 00272 _Tp __sum = __term; 00273 00274 const unsigned int __max_iter = 1000; 00275 for (unsigned int __i = 2; __i < __max_iter; ++__i) 00276 { 00277 const _Tp __i2m = 2 * __i - 1; 00278 const _Tp __i2 = 2 * __i; 00279 __term *= __i2m * __i2m * __kk / (__i2 * __i2); 00280 if (__term < std::numeric_limits<_Tp>::epsilon()) 00281 break; 00282 __sum += __term / __i2m; 00283 } 00284 00285 return __numeric_constants<_Tp>::__pi_2() * (_Tp(1) - __sum); 00286 } 00287 00288 00289 /** 00290 * @brief Return the Carlson elliptic function of the second kind 00291 * @f$ R_D(x,y,z) = R_J(x,y,z,z) @f$ where 00292 * @f$ R_J(x,y,z,p) @f$ is the Carlson elliptic function 00293 * of the third kind. 00294 * 00295 * The Carlson elliptic function of the second kind is defined by: 00296 * @f[ 00297 * R_D(x,y,z) = \frac{3}{2} \int_0^\infty 00298 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{3/2}} 00299 * @f] 00300 * 00301 * Based on Carlson's algorithms: 00302 * - B. C. Carlson Numer. Math. 33, 1 (1979) 00303 * - B. C. Carlson, Special Functions of Applied Mathematics (1977) 00304 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 00305 * by Press, Teukolsky, Vetterling, Flannery (1992) 00306 * 00307 * @param __x The first of two symmetric arguments. 00308 * @param __y The second of two symmetric arguments. 00309 * @param __z The third argument. 00310 * @return The Carlson elliptic function of the second kind. 00311 */ 00312 template<typename _Tp> 00313 _Tp 00314 __ellint_rd(const _Tp __x, const _Tp __y, const _Tp __z) 00315 { 00316 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 00317 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); 00318 const _Tp __min = std::numeric_limits<_Tp>::min(); 00319 const _Tp __max = std::numeric_limits<_Tp>::max(); 00320 const _Tp __lolim = _Tp(2) / std::pow(__max, _Tp(2) / _Tp(3)); 00321 const _Tp __uplim = std::pow(_Tp(0.1L) * __errtol / __min, _Tp(2) / _Tp(3)); 00322 00323 if (__x < _Tp(0) || __y < _Tp(0)) 00324 std::__throw_domain_error(__N("Argument less than zero " 00325 "in __ellint_rd.")); 00326 else if (__x + __y < __lolim || __z < __lolim) 00327 std::__throw_domain_error(__N("Argument too small " 00328 "in __ellint_rd.")); 00329 else 00330 { 00331 const _Tp __c0 = _Tp(1) / _Tp(4); 00332 const _Tp __c1 = _Tp(3) / _Tp(14); 00333 const _Tp __c2 = _Tp(1) / _Tp(6); 00334 const _Tp __c3 = _Tp(9) / _Tp(22); 00335 const _Tp __c4 = _Tp(3) / _Tp(26); 00336 00337 _Tp __xn = __x; 00338 _Tp __yn = __y; 00339 _Tp __zn = __z; 00340 _Tp __sigma = _Tp(0); 00341 _Tp __power4 = _Tp(1); 00342 00343 _Tp __mu; 00344 _Tp __xndev, __yndev, __zndev; 00345 00346 const unsigned int __max_iter = 100; 00347 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) 00348 { 00349 __mu = (__xn + __yn + _Tp(3) * __zn) / _Tp(5); 00350 __xndev = (__mu - __xn) / __mu; 00351 __yndev = (__mu - __yn) / __mu; 00352 __zndev = (__mu - __zn) / __mu; 00353 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); 00354 __epsilon = std::max(__epsilon, std::abs(__zndev)); 00355 if (__epsilon < __errtol) 00356 break; 00357 _Tp __xnroot = std::sqrt(__xn); 00358 _Tp __ynroot = std::sqrt(__yn); 00359 _Tp __znroot = std::sqrt(__zn); 00360 _Tp __lambda = __xnroot * (__ynroot + __znroot) 00361 + __ynroot * __znroot; 00362 __sigma += __power4 / (__znroot * (__zn + __lambda)); 00363 __power4 *= __c0; 00364 __xn = __c0 * (__xn + __lambda); 00365 __yn = __c0 * (__yn + __lambda); 00366 __zn = __c0 * (__zn + __lambda); 00367 } 00368 00369 // Note: __ea is an SPU badname. 00370 _Tp __eaa = __xndev * __yndev; 00371 _Tp __eb = __zndev * __zndev; 00372 _Tp __ec = __eaa - __eb; 00373 _Tp __ed = __eaa - _Tp(6) * __eb; 00374 _Tp __ef = __ed + __ec + __ec; 00375 _Tp __s1 = __ed * (-__c1 + __c3 * __ed 00376 / _Tp(3) - _Tp(3) * __c4 * __zndev * __ef 00377 / _Tp(2)); 00378 _Tp __s2 = __zndev 00379 * (__c2 * __ef 00380 + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa)); 00381 00382 return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2) 00383 / (__mu * std::sqrt(__mu)); 00384 } 00385 } 00386 00387 00388 /** 00389 * @brief Return the complete elliptic integral of the second kind 00390 * @f$ E(k) @f$ using the Carlson formulation. 00391 * 00392 * The complete elliptic integral of the second kind is defined as 00393 * @f[ 00394 * E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta} 00395 * @f] 00396 * 00397 * @param __k The argument of the complete elliptic function. 00398 * @return The complete elliptic function of the second kind. 00399 */ 00400 template<typename _Tp> 00401 _Tp 00402 __comp_ellint_2(const _Tp __k) 00403 { 00404 00405 if (__isnan(__k)) 00406 return std::numeric_limits<_Tp>::quiet_NaN(); 00407 else if (std::abs(__k) == 1) 00408 return _Tp(1); 00409 else if (std::abs(__k) > _Tp(1)) 00410 std::__throw_domain_error(__N("Bad argument in __comp_ellint_2.")); 00411 else 00412 { 00413 const _Tp __kk = __k * __k; 00414 00415 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1)) 00416 - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3); 00417 } 00418 } 00419 00420 00421 /** 00422 * @brief Return the incomplete elliptic integral of the second kind 00423 * @f$ E(k,\phi) @f$ using the Carlson formulation. 00424 * 00425 * The incomplete elliptic integral of the second kind is defined as 00426 * @f[ 00427 * E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta} 00428 * @f] 00429 * 00430 * @param __k The argument of the elliptic function. 00431 * @param __phi The integral limit argument of the elliptic function. 00432 * @return The elliptic function of the second kind. 00433 */ 00434 template<typename _Tp> 00435 _Tp 00436 __ellint_2(const _Tp __k, const _Tp __phi) 00437 { 00438 00439 if (__isnan(__k) || __isnan(__phi)) 00440 return std::numeric_limits<_Tp>::quiet_NaN(); 00441 else if (std::abs(__k) > _Tp(1)) 00442 std::__throw_domain_error(__N("Bad argument in __ellint_2.")); 00443 else 00444 { 00445 // Reduce phi to -pi/2 < phi < +pi/2. 00446 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() 00447 + _Tp(0.5L)); 00448 const _Tp __phi_red = __phi 00449 - __n * __numeric_constants<_Tp>::__pi(); 00450 00451 const _Tp __kk = __k * __k; 00452 const _Tp __s = std::sin(__phi_red); 00453 const _Tp __ss = __s * __s; 00454 const _Tp __sss = __ss * __s; 00455 const _Tp __c = std::cos(__phi_red); 00456 const _Tp __cc = __c * __c; 00457 00458 const _Tp __E = __s 00459 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1)) 00460 - __kk * __sss 00461 * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1)) 00462 / _Tp(3); 00463 00464 if (__n == 0) 00465 return __E; 00466 else 00467 return __E + _Tp(2) * __n * __comp_ellint_2(__k); 00468 } 00469 } 00470 00471 00472 /** 00473 * @brief Return the Carlson elliptic function 00474 * @f$ R_C(x,y) = R_F(x,y,y) @f$ where @f$ R_F(x,y,z) @f$ 00475 * is the Carlson elliptic function of the first kind. 00476 * 00477 * The Carlson elliptic function is defined by: 00478 * @f[ 00479 * R_C(x,y) = \frac{1}{2} \int_0^\infty 00480 * \frac{dt}{(t + x)^{1/2}(t + y)} 00481 * @f] 00482 * 00483 * Based on Carlson's algorithms: 00484 * - B. C. Carlson Numer. Math. 33, 1 (1979) 00485 * - B. C. Carlson, Special Functions of Applied Mathematics (1977) 00486 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 00487 * by Press, Teukolsky, Vetterling, Flannery (1992) 00488 * 00489 * @param __x The first argument. 00490 * @param __y The second argument. 00491 * @return The Carlson elliptic function. 00492 */ 00493 template<typename _Tp> 00494 _Tp 00495 __ellint_rc(const _Tp __x, const _Tp __y) 00496 { 00497 const _Tp __min = std::numeric_limits<_Tp>::min(); 00498 const _Tp __max = std::numeric_limits<_Tp>::max(); 00499 const _Tp __lolim = _Tp(5) * __min; 00500 const _Tp __uplim = __max / _Tp(5); 00501 00502 if (__x < _Tp(0) || __y < _Tp(0) || __x + __y < __lolim) 00503 std::__throw_domain_error(__N("Argument less than zero " 00504 "in __ellint_rc.")); 00505 else 00506 { 00507 const _Tp __c0 = _Tp(1) / _Tp(4); 00508 const _Tp __c1 = _Tp(1) / _Tp(7); 00509 const _Tp __c2 = _Tp(9) / _Tp(22); 00510 const _Tp __c3 = _Tp(3) / _Tp(10); 00511 const _Tp __c4 = _Tp(3) / _Tp(8); 00512 00513 _Tp __xn = __x; 00514 _Tp __yn = __y; 00515 00516 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 00517 const _Tp __errtol = std::pow(__eps / _Tp(30), _Tp(1) / _Tp(6)); 00518 _Tp __mu; 00519 _Tp __sn; 00520 00521 const unsigned int __max_iter = 100; 00522 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) 00523 { 00524 __mu = (__xn + _Tp(2) * __yn) / _Tp(3); 00525 __sn = (__yn + __mu) / __mu - _Tp(2); 00526 if (std::abs(__sn) < __errtol) 00527 break; 00528 const _Tp __lambda = _Tp(2) * std::sqrt(__xn) * std::sqrt(__yn) 00529 + __yn; 00530 __xn = __c0 * (__xn + __lambda); 00531 __yn = __c0 * (__yn + __lambda); 00532 } 00533 00534 _Tp __s = __sn * __sn 00535 * (__c3 + __sn*(__c1 + __sn * (__c4 + __sn * __c2))); 00536 00537 return (_Tp(1) + __s) / std::sqrt(__mu); 00538 } 00539 } 00540 00541 00542 /** 00543 * @brief Return the Carlson elliptic function @f$ R_J(x,y,z,p) @f$ 00544 * of the third kind. 00545 * 00546 * The Carlson elliptic function of the third kind is defined by: 00547 * @f[ 00548 * R_J(x,y,z,p) = \frac{3}{2} \int_0^\infty 00549 * \frac{dt}{(t + x)^{1/2}(t + y)^{1/2}(t + z)^{1/2}(t + p)} 00550 * @f] 00551 * 00552 * Based on Carlson's algorithms: 00553 * - B. C. Carlson Numer. Math. 33, 1 (1979) 00554 * - B. C. Carlson, Special Functions of Applied Mathematics (1977) 00555 * - Numerical Recipes in C, 2nd ed, pp. 261-269, 00556 * by Press, Teukolsky, Vetterling, Flannery (1992) 00557 * 00558 * @param __x The first of three symmetric arguments. 00559 * @param __y The second of three symmetric arguments. 00560 * @param __z The third of three symmetric arguments. 00561 * @param __p The fourth argument. 00562 * @return The Carlson elliptic function of the fourth kind. 00563 */ 00564 template<typename _Tp> 00565 _Tp 00566 __ellint_rj(const _Tp __x, const _Tp __y, const _Tp __z, const _Tp __p) 00567 { 00568 const _Tp __min = std::numeric_limits<_Tp>::min(); 00569 const _Tp __max = std::numeric_limits<_Tp>::max(); 00570 const _Tp __lolim = std::pow(_Tp(5) * __min, _Tp(1)/_Tp(3)); 00571 const _Tp __uplim = _Tp(0.3L) 00572 * std::pow(_Tp(0.2L) * __max, _Tp(1)/_Tp(3)); 00573 00574 if (__x < _Tp(0) || __y < _Tp(0) || __z < _Tp(0)) 00575 std::__throw_domain_error(__N("Argument less than zero " 00576 "in __ellint_rj.")); 00577 else if (__x + __y < __lolim || __x + __z < __lolim 00578 || __y + __z < __lolim || __p < __lolim) 00579 std::__throw_domain_error(__N("Argument too small " 00580 "in __ellint_rj")); 00581 else 00582 { 00583 const _Tp __c0 = _Tp(1) / _Tp(4); 00584 const _Tp __c1 = _Tp(3) / _Tp(14); 00585 const _Tp __c2 = _Tp(1) / _Tp(3); 00586 const _Tp __c3 = _Tp(3) / _Tp(22); 00587 const _Tp __c4 = _Tp(3) / _Tp(26); 00588 00589 _Tp __xn = __x; 00590 _Tp __yn = __y; 00591 _Tp __zn = __z; 00592 _Tp __pn = __p; 00593 _Tp __sigma = _Tp(0); 00594 _Tp __power4 = _Tp(1); 00595 00596 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 00597 const _Tp __errtol = std::pow(__eps / _Tp(8), _Tp(1) / _Tp(6)); 00598 00599 _Tp __lambda, __mu; 00600 _Tp __xndev, __yndev, __zndev, __pndev; 00601 00602 const unsigned int __max_iter = 100; 00603 for (unsigned int __iter = 0; __iter < __max_iter; ++__iter) 00604 { 00605 __mu = (__xn + __yn + __zn + _Tp(2) * __pn) / _Tp(5); 00606 __xndev = (__mu - __xn) / __mu; 00607 __yndev = (__mu - __yn) / __mu; 00608 __zndev = (__mu - __zn) / __mu; 00609 __pndev = (__mu - __pn) / __mu; 00610 _Tp __epsilon = std::max(std::abs(__xndev), std::abs(__yndev)); 00611 __epsilon = std::max(__epsilon, std::abs(__zndev)); 00612 __epsilon = std::max(__epsilon, std::abs(__pndev)); 00613 if (__epsilon < __errtol) 00614 break; 00615 const _Tp __xnroot = std::sqrt(__xn); 00616 const _Tp __ynroot = std::sqrt(__yn); 00617 const _Tp __znroot = std::sqrt(__zn); 00618 const _Tp __lambda = __xnroot * (__ynroot + __znroot) 00619 + __ynroot * __znroot; 00620 const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot) 00621 + __xnroot * __ynroot * __znroot; 00622 const _Tp __alpha2 = __alpha1 * __alpha1; 00623 const _Tp __beta = __pn * (__pn + __lambda) 00624 * (__pn + __lambda); 00625 __sigma += __power4 * __ellint_rc(__alpha2, __beta); 00626 __power4 *= __c0; 00627 __xn = __c0 * (__xn + __lambda); 00628 __yn = __c0 * (__yn + __lambda); 00629 __zn = __c0 * (__zn + __lambda); 00630 __pn = __c0 * (__pn + __lambda); 00631 } 00632 00633 // Note: __ea is an SPU badname. 00634 _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev; 00635 _Tp __eb = __xndev * __yndev * __zndev; 00636 _Tp __ec = __pndev * __pndev; 00637 _Tp __e2 = __eaa - _Tp(3) * __ec; 00638 _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec); 00639 _Tp __s1 = _Tp(1) + __e2 * (-__c1 + _Tp(3) * __c3 * __e2 / _Tp(4) 00640 - _Tp(3) * __c4 * __e3 / _Tp(2)); 00641 _Tp __s2 = __eb * (__c2 / _Tp(2) 00642 + __pndev * (-__c3 - __c3 + __pndev * __c4)); 00643 _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3) 00644 - __c2 * __pndev * __ec; 00645 00646 return _Tp(3) * __sigma + __power4 * (__s1 + __s2 + __s3) 00647 / (__mu * std::sqrt(__mu)); 00648 } 00649 } 00650 00651 00652 /** 00653 * @brief Return the complete elliptic integral of the third kind 00654 * @f$ \Pi(k,\nu) = \Pi(k,\nu,\pi/2) @f$ using the 00655 * Carlson formulation. 00656 * 00657 * The complete elliptic integral of the third kind is defined as 00658 * @f[ 00659 * \Pi(k,\nu) = \int_0^{\pi/2} 00660 * \frac{d\theta} 00661 * {(1 - \nu \sin^2\theta)\sqrt{1 - k^2 \sin^2\theta}} 00662 * @f] 00663 * 00664 * @param __k The argument of the elliptic function. 00665 * @param __nu The second argument of the elliptic function. 00666 * @return The complete elliptic function of the third kind. 00667 */ 00668 template<typename _Tp> 00669 _Tp 00670 __comp_ellint_3(const _Tp __k, const _Tp __nu) 00671 { 00672 00673 if (__isnan(__k) || __isnan(__nu)) 00674 return std::numeric_limits<_Tp>::quiet_NaN(); 00675 else if (__nu == _Tp(1)) 00676 return std::numeric_limits<_Tp>::infinity(); 00677 else if (std::abs(__k) > _Tp(1)) 00678 std::__throw_domain_error(__N("Bad argument in __comp_ellint_3.")); 00679 else 00680 { 00681 const _Tp __kk = __k * __k; 00682 00683 return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1)) 00684 - __nu 00685 * __ellint_rj(_Tp(0), _Tp(1) - __kk, _Tp(1), _Tp(1) + __nu) 00686 / _Tp(3); 00687 } 00688 } 00689 00690 00691 /** 00692 * @brief Return the incomplete elliptic integral of the third kind 00693 * @f$ \Pi(k,\nu,\phi) @f$ using the Carlson formulation. 00694 * 00695 * The incomplete elliptic integral of the third kind is defined as 00696 * @f[ 00697 * \Pi(k,\nu,\phi) = \int_0^{\phi} 00698 * \frac{d\theta} 00699 * {(1 - \nu \sin^2\theta) 00700 * \sqrt{1 - k^2 \sin^2\theta}} 00701 * @f] 00702 * 00703 * @param __k The argument of the elliptic function. 00704 * @param __nu The second argument of the elliptic function. 00705 * @param __phi The integral limit argument of the elliptic function. 00706 * @return The elliptic function of the third kind. 00707 */ 00708 template<typename _Tp> 00709 _Tp 00710 __ellint_3(const _Tp __k, const _Tp __nu, const _Tp __phi) 00711 { 00712 00713 if (__isnan(__k) || __isnan(__nu) || __isnan(__phi)) 00714 return std::numeric_limits<_Tp>::quiet_NaN(); 00715 else if (std::abs(__k) > _Tp(1)) 00716 std::__throw_domain_error(__N("Bad argument in __ellint_3.")); 00717 else 00718 { 00719 // Reduce phi to -pi/2 < phi < +pi/2. 00720 const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi() 00721 + _Tp(0.5L)); 00722 const _Tp __phi_red = __phi 00723 - __n * __numeric_constants<_Tp>::__pi(); 00724 00725 const _Tp __kk = __k * __k; 00726 const _Tp __s = std::sin(__phi_red); 00727 const _Tp __ss = __s * __s; 00728 const _Tp __sss = __ss * __s; 00729 const _Tp __c = std::cos(__phi_red); 00730 const _Tp __cc = __c * __c; 00731 00732 const _Tp __Pi = __s 00733 * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1)) 00734 - __nu * __sss 00735 * __ellint_rj(__cc, _Tp(1) - __kk * __ss, _Tp(1), 00736 _Tp(1) + __nu * __ss) / _Tp(3); 00737 00738 if (__n == 0) 00739 return __Pi; 00740 else 00741 return __Pi + _Tp(2) * __n * __comp_ellint_3(__k, __nu); 00742 } 00743 } 00744 00745 } // namespace std::tr1::__detail 00746 } 00747 } 00748 00749 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC 00750