PR libstdc++/83140 - assoc_legendre returns negated value when m is odd

2018-05-10  Edward Smith-Rowland  <3dw4rd@verizon.net>

	PR libstdc++/83140 - assoc_legendre returns negated value when m is odd
	* include/tr1/legendre_function.tcc (__assoc_legendre_p): Add __phase
	argument defaulted to +1.  Doxy comments on same.
	* testsuite/special_functions/02_assoc_legendre/
	check_assoc_legendre.cc: Regen.
	* testsuite/tr1/5_numerical_facilities/special_functions/
	02_assoc_legendre/check_tr1_assoc_legendre.cc: Regen.

From-SVN: r260115

libstdc++-v3/ChangeLog

@@ -1,3 +1,13 @@
+2018-05-10  Edward Smith-Rowland  <3dw4rd@verizon.net>
+
+	PR libstdc++/83140 - assoc_legendre returns negated value when m is odd
+	* include/tr1/legendre_function.tcc (__assoc_legendre_p): Add __phase
+	argument defaulted to +1.  Doxy comments on same.
+	* testsuite/special_functions/02_assoc_legendre/
+	check_assoc_legendre.cc: Regen.
+	* testsuite/tr1/5_numerical_facilities/special_functions/
+	02_assoc_legendre/check_tr1_assoc_legendre.cc: Regen.
+
 2018-05-10  Jonathan Wakely  <jwakely@redhat.com>
 
 	PR libstdc++/85729

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

@@ -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>
@@ -127,16 +127,19 @@ namespace tr1
      *     P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x)
      *   @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))
@@ -160,7 +163,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,8 +208,10 @@ 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.
      * 
-     *   @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$.
@@ -265,19 +270,15 @@ 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)
+          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;
-            }
           else if (__l == __m + 1)
-            {
             return __y_mp1m;
-            }
           else
             {
               _Tp __y_lm = _Tp(0);