* [libstc++] Don't throw in std::assoc_legendre for m > l
@ 2019-03-03 21:09 André Brand
0 siblings, 0 replies; 3+ messages in thread
From: André Brand @ 2019-03-03 21:09 UTC (permalink / raw)
To: gcc-patches, libstdc++
[-- Attachment #1: Type: text/plain, Size: 567 bytes --]
The return value specified in "8.1.2 associated Legendre polynomials"
of ISO/IEC JTC 1/SC 22/WG 21 N3060 (which is identical to the
expression in the doxygen comment of the patched function) is well-
defined for m>l: it is always zero because $ P_l(x) $ is a polynomial
of degree l.
The standard does not enforce an exception in this case because none of
the requirements in 8.1 (5) on page 11 of ISO/IEC JTC 1/SC 22/WG 21
N3060 are met.
Note: the implementation of st::assoc_legendre in Visual Studio 2017
(tested with Visual Studio 15.9.7) silently returns zero.
[-- Attachment #2: Type: text/x-patch, Size: 3096 bytes --]
Index: libstdc++-v3/include/tr1/legendre_function.tcc
===================================================================
--- libstdc++-v3/include/tr1/legendre_function.tcc (revision 269352)
+++ libstdc++-v3/include/tr1/legendre_function.tcc (working copy)
@@ -67,13 +67,13 @@
/**
* @brief Return the Legendre polynomial by recursion on degree
* @f$ l @f$.
- *
+ *
* The Legendre function of @f$ l @f$ and @f$ x @f$,
* @f$ P_l(x) @f$, is defined by:
* @f[
* P_l(x) = \frac{1}{2^l l!}\frac{d^l}{dx^l}(x^2 - 1)^{l}
* @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$.
*/
@@ -120,17 +120,17 @@
/**
* @brief Return the associated Legendre function by recursion
* on @f$ l @f$.
- *
+ *
* The associated Legendre function is derived from the Legendre function
* @f$ P_l(x) @f$ by the Rodrigues formula:
* @f[
* P_l^m(x) = (1 - x^2)^{m/2}\frac{d^m}{dx^m}P_l(x)
* @f]
- *
+ *
* @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$.
+ * @f$ m >= 0 @f$.
* @param x The argument of the associated Legendre function.
* @f$ |x| <= 1 @f$.
* @param phase The phase of the associated Legendre function.
@@ -146,8 +146,7 @@
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."));
+ return _Tp(0);
else if (__isnan(__x))
return std::numeric_limits<_Tp>::quiet_NaN();
else if (__m == 0)
@@ -192,7 +191,7 @@
/**
* @brief Return the spherical associated Legendre function.
- *
+ *
* The spherical associated Legendre function of @f$ l @f$, @f$ m @f$,
* and @f$ \theta @f$ is defined as @f$ Y_l^m(\theta,0) @f$ where
* @f[
@@ -202,7 +201,7 @@
* @f]
* is the spherical harmonic function and @f$ P_l^m(x) @f$ is the
* associated Legendre function.
- *
+ *
* This function differs from the associated Legendre function by
* argument (@f$x = \cos(\theta)@f$) and by a normalization factor
* but this factor is rather large for large @f$ l @f$ and @f$ m @f$
@@ -210,7 +209,7 @@
* 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 degree of the spherical associated Legendre function.
* @f$ l >= 0 @f$.
* @param m The order of the spherical associated Legendre function.
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: [libstc++] Don't throw in std::assoc_legendre for m > l
2019-03-04 18:57 Ed Smith-Rowland
@ 2019-03-06 1:47 ` Jonathan Wakely
0 siblings, 0 replies; 3+ messages in thread
From: Jonathan Wakely @ 2019-03-06 1:47 UTC (permalink / raw)
To: Ed Smith-Rowland; +Cc: libstdc++, gcc-patches, andre.j.brand
On 04/03/19 13:57 -0500, Ed Smith-Rowland wrote:
>This is actually PR libstdc++/86655.
>
>Thank you for reminding me Andre.
>
>I remove the throw for m > l and just return 0. This is also done for
>sph_legendre.
>
>This build and tests clean on x86_64-linux.
>
>OK?
OK for trunk, thanks.
^ permalink raw reply [flat|nested] 3+ messages in thread
* Re: [libstc++] Don't throw in std::assoc_legendre for m > l
@ 2019-03-04 18:57 Ed Smith-Rowland
2019-03-06 1:47 ` Jonathan Wakely
0 siblings, 1 reply; 3+ messages in thread
From: Ed Smith-Rowland @ 2019-03-04 18:57 UTC (permalink / raw)
To: libstdc++, gcc-patches, andre.j.brand, Jonathan Wakely
[-- Attachment #1: Type: text/plain, Size: 218 bytes --]
This is actually PR libstdc++/86655.
Thank you for reminding me Andre.
I remove the throw for m > l and just return 0. This is also done for
sph_legendre.
This build and tests clean on x86_64-linux.
OK?
Ed
[-- Attachment #2: CL_pr86655 --]
[-- Type: text/plain, Size: 662 bytes --]
2018-03-04 Edward Smith-Rowland <3dw4rd@verizon.net>
PR libstdc++/86655 - std::assoc_legendre should not constrain
the value of m (or x).
* include/tr1/legendre_function.tcc (__assoc_legendre_p,
__sph_legendre): If degree > order Don't throw, return 0.
(__legendre_p, __assoc_legendre_p): Don't constrain x either.
* testsuite/special_functions/02_assoc_legendre/pr86655.cc: New test.
* testsuite/special_functions/20_sph_legendre/pr86655.cc: New test.
* testsuite/tr1/5_numerical_facilities/special_functions/
02_assoc_legendre/pr86655.cc: New test.
* testsuite/tr1/5_numerical_facilities/special_functions/
22_sph_legendre/pr86655.cc: New test.
[-- Attachment #3: patch_pr86655 --]
[-- Type: text/plain, Size: 11187 bytes --]
Index: include/tr1/legendre_function.tcc
===================================================================
--- include/tr1/legendre_function.tcc (revision 269347)
+++ include/tr1/legendre_function.tcc (working copy)
@@ -82,10 +82,7 @@
__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,11 +123,11 @@
* @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 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.
@@ -142,12 +139,8 @@
_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)
@@ -209,12 +202,12 @@
* 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.
+ * 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 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.
*/
@@ -227,11 +220,8 @@
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);
@@ -284,7 +274,7 @@
_Tp __y_lm = _Tp(0);
// Compute Y_l^m, l > m+1, upward recursion on l.
- for (unsigned 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);
Index: testsuite/special_functions/02_assoc_legendre/pr86655.cc
===================================================================
--- testsuite/special_functions/02_assoc_legendre/pr86655.cc (nonexistent)
+++ testsuite/special_functions/02_assoc_legendre/pr86655.cc (working copy)
@@ -0,0 +1,56 @@
+// { dg-do run { target c++11 } }
+// { dg-options "-D__STDCPP_WANT_MATH_SPEC_FUNCS__ -ffp-contract=off" }
+
+// Copyright (C) 2019 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
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3. If not see
+// <http://www.gnu.org/licenses/>.
+
+#include <initializer_list>
+#include <cmath>
+#if defined(__TEST_DEBUG)
+# include <iostream>
+# define VERIFY(A) \
+ if (!(A)) \
+ { \
+ std::cout << "line " << __LINE__ \
+ << " std::assoc_legendre(l, m, x) == 0: " << (A) \
+ << '\n'; \
+ }
+#else
+# include <testsuite_hooks.h>
+#endif
+
+template<typename _Tp>
+ void
+ test_m_gt_l()
+ {
+ bool test __attribute__((unused)) = true;
+ for (auto l : {0u, 1u, 2u, 5u})
+ for (auto m : {l + 1u, l + 2u})
+ for (auto i : {-2, -1, 0, 1, 2})
+ {
+ auto x = _Tp(i * 0.5L);
+ VERIFY(std::assoc_legendre(l, m, x) == _Tp(0));
+ }
+ }
+
+int
+main()
+{
+ test_m_gt_l<float>();
+ test_m_gt_l<double>();
+ test_m_gt_l<long double>();
+}
Index: testsuite/special_functions/20_sph_legendre/pr86655.cc
===================================================================
--- testsuite/special_functions/20_sph_legendre/pr86655.cc (nonexistent)
+++ testsuite/special_functions/20_sph_legendre/pr86655.cc (working copy)
@@ -0,0 +1,56 @@
+// { dg-do run { target c++11 } }
+// { dg-options "-D__STDCPP_WANT_MATH_SPEC_FUNCS__ -ffp-contract=off" }
+
+// Copyright (C) 2019 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
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3. If not see
+// <http://www.gnu.org/licenses/>.
+
+#include <initializer_list>
+#include <cmath>
+#if defined(__TEST_DEBUG)
+# include <iostream>
+# define VERIFY(A) \
+ if (!(A)) \
+ { \
+ std::cout << "line " << __LINE__ \
+ << " std::sph_legendre(l, m, x) == 0: " << (A) \
+ << '\n'; \
+ }
+#else
+# include <testsuite_hooks.h>
+#endif
+
+template<typename _Tp>
+ void
+ test_m_gt_l()
+ {
+ bool test __attribute__((unused)) = true;
+ for (auto l : {0u, 1u, 2u, 5u})
+ for (auto m : {l + 1u, l + 2u})
+ for (auto i : {-2, -1, 0, 1, 2})
+ {
+ auto theta = std::acos(_Tp(i * 0.5L));
+ VERIFY(std::sph_legendre(l, m, theta) == _Tp(0));
+ }
+ }
+
+int
+main()
+{
+ test_m_gt_l<float>();
+ test_m_gt_l<double>();
+ test_m_gt_l<long double>();
+}
Index: testsuite/tr1/5_numerical_facilities/special_functions/02_assoc_legendre/pr86655.cc
===================================================================
--- testsuite/tr1/5_numerical_facilities/special_functions/02_assoc_legendre/pr86655.cc (nonexistent)
+++ testsuite/tr1/5_numerical_facilities/special_functions/02_assoc_legendre/pr86655.cc (working copy)
@@ -0,0 +1,56 @@
+// { dg-do run }
+// { dg-options "-std=c++98 -ffp-contract=off" }
+
+// Copyright (C) 2019 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
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3. If not see
+// <http://www.gnu.org/licenses/>.
+
+#include <tr1/cmath>
+#if defined(__TEST_DEBUG)
+# include <iostream>
+# define VERIFY(A) \
+ if (!(A)) \
+ { \
+ std::cout << "line " << __LINE__ \
+ << " std::tr1::assoc_legendre(l, m, x) == 0: " << (A) \
+ << '\n'; \
+ }
+#else
+# include <testsuite_hooks.h>
+#endif
+
+template<typename _Tp>
+ void
+ test_m_gt_l()
+ {
+ bool test __attribute__((unused)) = true;
+ unsigned int larr[4] = {0u, 1u, 2u, 5u};
+ for (unsigned int l = 0; l < 4; ++l)
+ for (unsigned int m = larr[l] + 1u; m <= larr[l] + 2u; ++m)
+ for (int i = -2; i <= +2; ++i)
+ {
+ _Tp x = _Tp(i * 0.5L);
+ VERIFY(std::tr1::assoc_legendre(larr[l], m, x) == _Tp(0));
+ }
+ }
+
+int
+main()
+{
+ test_m_gt_l<float>();
+ test_m_gt_l<double>();
+ test_m_gt_l<long double>();
+}
Index: testsuite/tr1/5_numerical_facilities/special_functions/22_sph_legendre/pr86655.cc
===================================================================
--- testsuite/tr1/5_numerical_facilities/special_functions/22_sph_legendre/pr86655.cc (nonexistent)
+++ testsuite/tr1/5_numerical_facilities/special_functions/22_sph_legendre/pr86655.cc (working copy)
@@ -0,0 +1,56 @@
+// { dg-do run }
+// { dg-options "-std=c++98 -ffp-contract=off" }
+
+// Copyright (C) 2019 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
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3. If not see
+// <http://www.gnu.org/licenses/>.
+
+#include <tr1/cmath>
+#if defined(__TEST_DEBUG)
+# include <iostream>
+# define VERIFY(A) \
+ if (!(A)) \
+ { \
+ std::cout << "line " << __LINE__ \
+ << " std::sph_legendre(l, m, x) == 0: " << (A) \
+ << '\n'; \
+ }
+#else
+# include <testsuite_hooks.h>
+#endif
+
+template<typename _Tp>
+ void
+ test_m_gt_l()
+ {
+ bool test __attribute__((unused)) = true;
+ unsigned int larr[4] = {0u, 1u, 2u, 5u};
+ for (unsigned int l = 0; l < 4; ++l)
+ for (unsigned int m = larr[l] + 1u; m <= larr[l] + 2u; ++m)
+ for (int i = -2; i <= +2; ++i)
+ {
+ _Tp theta = std::acos(_Tp(i * 0.5L));
+ VERIFY(std::tr1::sph_legendre(larr[l], m, theta) == _Tp(0));
+ }
+ }
+
+int
+main()
+{
+ test_m_gt_l<float>();
+ test_m_gt_l<double>();
+ test_m_gt_l<long double>();
+}
^ permalink raw reply [flat|nested] 3+ messages in thread
end of thread, other threads:[~2019-03-05 23:26 UTC | newest]
Thread overview: 3+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2019-03-03 21:09 [libstc++] Don't throw in std::assoc_legendre for m > l André Brand
2019-03-04 18:57 Ed Smith-Rowland
2019-03-06 1:47 ` Jonathan Wakely
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).