| b34f60ac |
1 | // Special functions -*- C++ -*- |
| 2 | |
| f1717362 |
3 | // Copyright (C) 2006-2016 Free Software Foundation, Inc. |
| b34f60ac |
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 |
| 6bc9506f |
8 | // Free Software Foundation; either version 3, or (at your option) |
| b34f60ac |
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 | // |
| 6bc9506f |
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/>. |
| b34f60ac |
24 | |
| 25 | /** @file tr1/ell_integral.tcc |
| 26 | * This is an internal header file, included by other library headers. |
| 5846aeac |
27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
| b34f60ac |
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 | |
| c17b0a1c |
42 | #ifndef _GLIBCXX_TR1_ELL_INTEGRAL_TCC |
| 43 | #define _GLIBCXX_TR1_ELL_INTEGRAL_TCC 1 |
| b34f60ac |
44 | |
| 2948dd21 |
45 | namespace std _GLIBCXX_VISIBILITY(default) |
| b34f60ac |
46 | { |
| c17b0a1c |
47 | namespace tr1 |
| 48 | { |
| b34f60ac |
49 | // [5.2] Special functions |
| 50 | |
| b34f60ac |
51 | // Implementation-space details. |
| b34f60ac |
52 | namespace __detail |
| 53 | { |
| 2948dd21 |
54 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
| b34f60ac |
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 |
| cd7f5f45 |
73 | __ellint_rf(_Tp __x, _Tp __y, _Tp __z) |
| b34f60ac |
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 |
| cd7f5f45 |
152 | __comp_ellint_1_series(_Tp __k) |
| b34f60ac |
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 |
| cd7f5f45 |
190 | __comp_ellint_1(_Tp __k) |
| b34f60ac |
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 |
| cd7f5f45 |
218 | __ellint_1(_Tp __k, _Tp __phi) |
| b34f60ac |
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 |
| cd7f5f45 |
265 | __comp_ellint_2_series(_Tp __k) |
| b34f60ac |
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) |
| 9fc1117c |
303 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
| b34f60ac |
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 |
| cd7f5f45 |
313 | __ellint_rd(_Tp __x, _Tp __y, _Tp __z) |
| b34f60ac |
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 | |
| 12bd5954 |
368 | // Note: __ea is an SPU badname. |
| 369 | _Tp __eaa = __xndev * __yndev; |
| b34f60ac |
370 | _Tp __eb = __zndev * __zndev; |
| 12bd5954 |
371 | _Tp __ec = __eaa - __eb; |
| 372 | _Tp __ed = __eaa - _Tp(6) * __eb; |
| b34f60ac |
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 |
| 12bd5954 |
379 | + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa)); |
| b34f60ac |
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 |
| cd7f5f45 |
401 | __comp_ellint_2(_Tp __k) |
| b34f60ac |
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 |
| cd7f5f45 |
435 | __ellint_2(_Tp __k, _Tp __phi) |
| b34f60ac |
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) |
| 9fc1117c |
485 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
| b34f60ac |
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 |
| cd7f5f45 |
494 | __ellint_rc(_Tp __x, _Tp __y) |
| b34f60ac |
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) |
| 9fc1117c |
554 | * - Numerical Recipes in C, 2nd ed, pp. 261-269, |
| b34f60ac |
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 |
| cd7f5f45 |
565 | __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p) |
| b34f60ac |
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; |
| 41456429 |
619 | const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot) |
| b34f60ac |
620 | + __xnroot * __ynroot * __znroot; |
| 41456429 |
621 | const _Tp __alpha2 = __alpha1 * __alpha1; |
| b34f60ac |
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 | |
| 12bd5954 |
632 | // Note: __ea is an SPU badname. |
| 633 | _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev; |
| b34f60ac |
634 | _Tp __eb = __xndev * __yndev * __zndev; |
| 635 | _Tp __ec = __pndev * __pndev; |
| 12bd5954 |
636 | _Tp __e2 = __eaa - _Tp(3) * __ec; |
| 637 | _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec); |
| b34f60ac |
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)); |
| 12bd5954 |
642 | _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3) |
| b34f60ac |
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 |
| cd7f5f45 |
669 | __comp_ellint_3(_Tp __k, _Tp __nu) |
| b34f60ac |
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 |
| cd7f5f45 |
709 | __ellint_3(_Tp __k, _Tp __nu, _Tp __phi) |
| b34f60ac |
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 | |
| 2948dd21 |
744 | _GLIBCXX_END_NAMESPACE_VERSION |
| b34f60ac |
745 | } // namespace std::tr1::__detail |
| c17b0a1c |
746 | } |
| b34f60ac |
747 | } |
| 748 | |
| c17b0a1c |
749 | #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC |
| b34f60ac |
750 | |
projects / thirdparty / gcc.git / blame