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

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