Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. // Special functions -*- C++ -*-
  2.  
  3. // Copyright (C) 2006-2015 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.  
  751.