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