l4/pkg/libstdc++-v3/contrib/libstdc++-v3-4.9/include/tr1/ell_integral.tcc

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006-2014 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 namespace tr1
  48 {
  49   // [5.2] Special functions
  50
  51   // Implementation-space details.
  52   namespace __detail
  53   {
  54   _GLIBCXX_BEGIN_NAMESPACE_VERSION
  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
  73     __ellint_rf(_Tp __x, _Tp __y, _Tp __z)
  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
 152     __comp_ellint_1_series(_Tp __k)
 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
 190     __comp_ellint_1(_Tp __k)
 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
 218     __ellint_1(_Tp __k, _Tp __phi)
 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
 265     __comp_ellint_2_series(_Tp __k)
 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)
 303      *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
 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
 313     __ellint_rd(_Tp __x, _Tp __y, _Tp __z)
 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
 368           // Note: __ea is an SPU badname.
 369           _Tp __eaa = __xndev * __yndev;
 370           _Tp __eb = __zndev * __zndev;
 371           _Tp __ec = __eaa - __eb;
 372           _Tp __ed = __eaa - _Tp(6) * __eb;
 373           _Tp __ef = __ed + __ec + __ec;
 374           _Tp __s1 = __ed * (-__c1 + __c3 * __ed
 375                                    / _Tp(3) - _Tp(3) * __c4 * __zndev * __ef
 376                                    / _Tp(2));
 377           _Tp __s2 = __zndev
 378                    * (__c2 * __ef
 379                     + __zndev * (-__c3 * __ec - __zndev * __c4 - __eaa));
 380
 381           return _Tp(3) * __sigma + __power4 * (_Tp(1) + __s1 + __s2)
 382                                         / (__mu * std::sqrt(__mu));
 383         }
 384     }
 385
 386
 387     /**
 388      *   @brief  Return the complete elliptic integral of the second kind
 389      *           @f$ E(k) @f$ using the Carlson formulation.
 390      *
 391      *   The complete elliptic integral of the second kind is defined as
 392      *   @f[
 393      *     E(k,\pi/2) = \int_0^{\pi/2}\sqrt{1 - k^2 sin^2\theta}
 394      *   @f]
 395      *
 396      *   @param  __k  The argument of the complete elliptic function.
 397      *   @return  The complete elliptic function of the second kind.
 398      */
 399     template<typename _Tp>
 400     _Tp
 401     __comp_ellint_2(_Tp __k)
 402     {
 403
 404       if (__isnan(__k))
 405         return std::numeric_limits<_Tp>::quiet_NaN();
 406       else if (std::abs(__k) == 1)
 407         return _Tp(1);
 408       else if (std::abs(__k) > _Tp(1))
 409         std::__throw_domain_error(__N("Bad argument in __comp_ellint_2."));
 410       else
 411         {
 412           const _Tp __kk = __k * __k;
 413
 414           return __ellint_rf(_Tp(0), _Tp(1) - __kk, _Tp(1))
 415                - __kk * __ellint_rd(_Tp(0), _Tp(1) - __kk, _Tp(1)) / _Tp(3);
 416         }
 417     }
 418
 419
 420     /**
 421      *   @brief  Return the incomplete elliptic integral of the second kind
 422      *           @f$ E(k,\phi) @f$ using the Carlson formulation.
 423      *
 424      *   The incomplete elliptic integral of the second kind is defined as
 425      *   @f[
 426      *     E(k,\phi) = \int_0^{\phi} \sqrt{1 - k^2 sin^2\theta}
 427      *   @f]
 428      *
 429      *   @param  __k  The argument of the elliptic function.
 430      *   @param  __phi  The integral limit argument of the elliptic function.
 431      *   @return  The elliptic function of the second kind.
 432      */
 433     template<typename _Tp>
 434     _Tp
 435     __ellint_2(_Tp __k, _Tp __phi)
 436     {
 437
 438       if (__isnan(__k) || __isnan(__phi))
 439         return std::numeric_limits<_Tp>::quiet_NaN();
 440       else if (std::abs(__k) > _Tp(1))
 441         std::__throw_domain_error(__N("Bad argument in __ellint_2."));
 442       else
 443         {
 444           //  Reduce phi to -pi/2 < phi < +pi/2.
 445           const int __n = std::floor(__phi / __numeric_constants<_Tp>::__pi()
 446                                    + _Tp(0.5L));
 447           const _Tp __phi_red = __phi
 448                               - __n * __numeric_constants<_Tp>::__pi();
 449
 450           const _Tp __kk = __k * __k;
 451           const _Tp __s = std::sin(__phi_red);
 452           const _Tp __ss = __s * __s;
 453           const _Tp __sss = __ss * __s;
 454           const _Tp __c = std::cos(__phi_red);
 455           const _Tp __cc = __c * __c;
 456
 457           const _Tp __E = __s
 458                         * __ellint_rf(__cc, _Tp(1) - __kk * __ss, _Tp(1))
 459                         - __kk * __sss
 460                         * __ellint_rd(__cc, _Tp(1) - __kk * __ss, _Tp(1))
 461                         / _Tp(3);
 462
 463           if (__n == 0)
 464             return __E;
 465           else
 466             return __E + _Tp(2) * __n * __comp_ellint_2(__k);
 467         }
 468     }
 469
 470
 471     /**
 472      *   @brief  Return the Carlson elliptic function
 473      *           @f$ R_C(x,y) = R_F(x,y,y) @f$ where @f$ R_F(x,y,z) @f$
 474      *           is the Carlson elliptic function of the first kind.
 475      *
 476      *   The Carlson elliptic function is defined by:
 477      *   @f[
 478      *       R_C(x,y) = \frac{1}{2} \int_0^\infty
 479      *                 \frac{dt}{(t + x)^{1/2}(t + y)}
 480      *   @f]
 481      *
 482      *   Based on Carlson's algorithms:
 483      *   -  B. C. Carlson Numer. Math. 33, 1 (1979)
 484      *   -  B. C. Carlson, Special Functions of Applied Mathematics (1977)
 485      *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
 486      *      by Press, Teukolsky, Vetterling, Flannery (1992)
 487      *
 488      *   @param  __x  The first argument.
 489      *   @param  __y  The second argument.
 490      *   @return  The Carlson elliptic function.
 491      */
 492     template<typename _Tp>
 493     _Tp
 494     __ellint_rc(_Tp __x, _Tp __y)
 495     {
 496       const _Tp __min = std::numeric_limits<_Tp>::min();
 497       const _Tp __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)
 554      *   -  Numerical Recipes in C, 2nd ed, pp. 261-269,
 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
 565     __ellint_rj(_Tp __x, _Tp __y, _Tp __z, _Tp __p)
 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;
 619               const _Tp __alpha1 = __pn * (__xnroot + __ynroot + __znroot)
 620                                 + __xnroot * __ynroot * __znroot;
 621               const _Tp __alpha2 = __alpha1 * __alpha1;
 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
 632           // Note: __ea is an SPU badname.
 633           _Tp __eaa = __xndev * (__yndev + __zndev) + __yndev * __zndev;
 634           _Tp __eb = __xndev * __yndev * __zndev;
 635           _Tp __ec = __pndev * __pndev;
 636           _Tp __e2 = __eaa - _Tp(3) * __ec;
 637           _Tp __e3 = __eb + _Tp(2) * __pndev * (__eaa - __ec);
 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));
 642           _Tp __s3 = __pndev * __eaa * (__c2 - __pndev * __c3)
 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
 669     __comp_ellint_3(_Tp __k, _Tp __nu)
 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
 709     __ellint_3(_Tp __k, _Tp __nu, _Tp __phi)
 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
 744   _GLIBCXX_END_NAMESPACE_VERSION
 745   } // namespace std::tr1::__detail
 746 }
 747 }
 748
 749 #endif // _GLIBCXX_TR1_ELL_INTEGRAL_TCC
 750