Update copyright years.

[thirdparty/gcc.git] / libstdc++-v3 / include / tr1 / legendre_function.tcc

diff --git a/libstdc++-v3/include/tr1/legendre_function.tcc b/libstdc++-v3/include/tr1/legendre_function.tcc
index e2d1880ecb003b4d9ddc6a8b4fbbd5904f98a16b..126355beb31bc4378fec63452a0685de05e21f3f 100644 (file)
--- a/libstdc++-v3/include/tr1/legendre_function.tcc
+++ b/libstdc++-v3/include/tr1/legendre_function.tcc
@@ -1,6 +1,6 @@
 // Special functions -*- C++ -*-
 
-// Copyright (C) 2006-2018 Free Software Foundation, Inc.
+// Copyright (C) 2006-2020 Free Software Foundation, Inc.
 //
 // This file is part of the GNU ISO C++ Library.  This library is free
 // software; you can redistribute it and/or modify it under the
@@ -44,7 +44,7 @@
 #ifndef _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC
 #define _GLIBCXX_TR1_LEGENDRE_FUNCTION_TCC 1
 
-#include "special_function_util.h"
+#include <tr1/special_function_util.h>
 
 namespace std _GLIBCXX_VISIBILITY(default)
 {
@@ -65,7 +65,7 @@ namespace tr1
   namespace __detail
   {
     /**
-     *   @brief  Return the Legendre polynomial by recursion on order
+     *   @brief  Return the Legendre polynomial by recursion on degree
      *           @f$ l @f$.
      * 
      *   The Legendre function of @f$ l @f$ and @f$ x @f$,
@@ -74,7 +74,7 @@ namespace tr1
      *     P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l}
      *   @f]
      * 
-     *   @param  l  The order of the Legendre polynomial.  @f$l >= 0@f$.
+     *   @param  l  The degree of the Legendre polynomial.  @f$l >= 0@f$.
      *   @param  x  The argument of the Legendre polynomial.  @f$|x| <= 1@f$.
      */
     template<typename _Tp>
@@ -82,10 +82,7 @@ namespace tr1
     __poly_legendre_p(unsigned int __l, _Tp __x)
     {
 
-      if ((__x < _Tp(-1)) || (__x > _Tp(+1)))
-        std::__throw_domain_error(__N("Argument out of range"
-                                      " in __poly_legendre_p."));
-      else if (__isnan(__x))
+      if (__isnan(__x))
         return std::numeric_limits<_Tp>::quiet_NaN();
       else if (__x == +_Tp(1))
         return +_Tp(1);
@@ -126,25 +123,24 @@ namespace tr1
      *   @f[
      *     P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x)
      *   @f]
+     *   @note @f$ P_l^m(x) = 0 @f$ if @f$ m > l @f$.
      * 
-     *   @param  l  The order of the associated Legendre function.
+     *   @param  l  The degree of the associated Legendre function.
      *              @f$ l >= 0 @f$.
      *   @param  m  The order of the associated Legendre function.
-     *              @f$ m <= l @f$.
      *   @param  x  The argument of the associated Legendre function.
      *              @f$ |x| <= 1 @f$.
+     *   @param  phase  The phase of the associated Legendre function.
+     *                  Use -1 for the Condon-Shortley phase convention.
      */
     template<typename _Tp>
     _Tp
-    __assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x)
+    __assoc_legendre_p(unsigned int __l, unsigned int __m, _Tp __x,
+                      _Tp __phase = _Tp(+1))
     {
 
-      if (__x < _Tp(-1) || __x > _Tp(+1))
-        std::__throw_domain_error(__N("Argument out of range"
-                                      " in __assoc_legendre_p."));
-      else if (__m > __l)
-        std::__throw_domain_error(__N("Degree out of range"
-                                      " in __assoc_legendre_p."));
+      if (__m > __l)
+        return _Tp(0);
       else if (__isnan(__x))
         return std::numeric_limits<_Tp>::quiet_NaN();
       else if (__m == 0)
@@ -160,7 +156,7 @@ namespace tr1
               _Tp __fact = _Tp(1);
               for (unsigned int __i = 1; __i <= __m; ++__i)
                 {
-                  __p_mm *= -__fact * __root;
+                  __p_mm *= __phase * __fact * __root;
                   __fact += _Tp(2);
                 }
             }
@@ -205,11 +201,13 @@ namespace tr1
      *   but this factor is rather large for large @f$ l @f$ and @f$ m @f$
      *   and so this function is stable for larger differences of @f$ l @f$
      *   and @f$ m @f$.
+     *   @note Unlike the case for __assoc_legendre_p the Condon-Shortley
+     *         phase factor @f$ (-1)^m @f$ is present here.
+     *   @note @f$ Y_l^m(\theta) = 0 @f$ if @f$ m > l @f$.
      * 
-     *   @param  l  The order of the spherical associated Legendre function.
+     *   @param  l  The degree of the spherical associated Legendre function.
      *              @f$ l >= 0 @f$.
      *   @param  m  The order of the spherical associated Legendre function.
-     *              @f$ m <= l @f$.
      *   @param  theta  The radian angle argument of the spherical associated
      *                  Legendre function.
      */
@@ -222,11 +220,8 @@ namespace tr1
 
       const _Tp __x = std::cos(__theta);
 
-      if (__l < __m)
-        {
-          std::__throw_domain_error(__N("Bad argument "
-                                        "in __sph_legendre."));
-        }
+      if (__m > __l)
+        return _Tp(0);
       else if (__m == 0)
         {
           _Tp __P = __poly_legendre_p(__l, __x);
@@ -265,25 +260,21 @@ namespace tr1
           const _Tp __lnpre_val =
                     -_Tp(0.25L) * __numeric_constants<_Tp>::__lnpi()
                     + _Tp(0.5L) * (__lnpoch + __m * __lncirc);
-          _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m)
-                   / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
+          const _Tp __sr = std::sqrt((_Tp(2) + _Tp(1) / __m)
+                         / (_Tp(4) * __numeric_constants<_Tp>::__pi()));
           _Tp __y_mm = __sgn * __sr * std::exp(__lnpre_val);
           _Tp __y_mp1m = __y_mp1m_factor * __y_mm;
 
           if (__l == __m)
-            {
-              return __y_mm;
-            }
+            return __y_mm;
           else if (__l == __m + 1)
-            {
-              return __y_mp1m;
-            }
+            return __y_mp1m;
           else
             {
               _Tp __y_lm = _Tp(0);
 
               // Compute Y_l^m, l > m+1, upward recursion on l.
-              for ( int __ll = __m + 2; __ll <= __l; ++__ll)
+              for (int __ll = __m + 2; __ll <= __l; ++__ll)
                 {
                   const _Tp __rat1 = _Tp(__ll - __m) / _Tp(__ll + __m);
                   const _Tp __rat2 = _Tp(__ll - __m - 1) / _Tp(__ll + __m - 1);