contrib/gcc-8.0/libstdc++-v3/include/tr1/ell_integral.tcc

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