libstdc++-v3/include/tr1/ell_integral.tcc

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