Skip to content

Commit 9c6716e

Browse files
PaulPaulXiCao
Paul
authored and
PaulXiCao
committed
implementation + unit
1 parent cc2fb58 commit 9c6716e

File tree

4 files changed

+169
-1
lines changed

4 files changed

+169
-1
lines changed

‎libcxx/include/__math/special_functions.h

Lines changed: 48 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,11 +11,13 @@
1111
#define _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
1212

1313
#include <__config>
14+
#include <__math/abs.h>
1415
#include <__math/copysign.h>
1516
#include <__math/traits.h>
1617
#include <__type_traits/enable_if.h>
1718
#include <__type_traits/is_integral.h>
1819
#include <limits>
20+
#include <stdexcept>
1921

2022
#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
2123
# pragma GCC system_header
@@ -49,7 +51,7 @@ _LIBCPP_HIDE_FROM_ABI _Real __hermite(unsigned __n, _Real __x) {
4951
}
5052

5153
if (!__math::isfinite(__H_n)) {
52-
// Overflow occured. Two possible cases:
54+
// Overflow occurred. Two possible cases:
5355
// n is odd: return infinity of the same sign as x.
5456
// n is even: return +Inf
5557
_Real __inf = std::numeric_limits<_Real>::infinity();
@@ -77,6 +79,51 @@ _LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, _Integer __x) {
7779
return std::hermite(__n, static_cast<double>(__x));
7880
}
7981

82+
template <class _Real>
83+
_LIBCPP_HIDE_FROM_ABI _Real __legendre(unsigned __n, _Real __x) {
84+
// The Legendre polynomial P_n(x).
85+
// The implementation is based on the recurrence formula: (n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x).
86+
// Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing.
87+
// Cambridge university press, 2007, p. 183.
88+
89+
// NOLINTBEGIN(readability-identifier-naming)
90+
if (__math::isnan(__x))
91+
return __x;
92+
93+
if (__math::fabs(__x) > 1)
94+
std::__throw_domain_error("Argument `x` of Legendre function is out of range: `|x| <= 1`.");
95+
96+
_Real __P_0{1};
97+
if (__n == 0)
98+
return __P_0;
99+
100+
_Real __P_n_prev = __P_0;
101+
_Real __P_n = __x;
102+
for (unsigned __i = 1; __i < __n; ++__i) {
103+
_Real __P_n_next = ((2 * __i + 1) * __x * __P_n - __i * __P_n_prev) / static_cast<_Real>(__i + 1);
104+
__P_n_prev = __P_n;
105+
__P_n = __P_n_next;
106+
}
107+
108+
return __P_n;
109+
// NOLINTEND(readability-identifier-naming)
110+
}
111+
112+
inline _LIBCPP_HIDE_FROM_ABI double legendre(unsigned __n, double __x) { return std::__legendre(__n, __x); }
113+
114+
inline _LIBCPP_HIDE_FROM_ABI float legendre(unsigned __n, float __x) { return std::__legendre(__n, __x); }
115+
116+
inline _LIBCPP_HIDE_FROM_ABI long double legendre(unsigned __n, long double __x) { return std::__legendre(__n, __x); }
117+
118+
inline _LIBCPP_HIDE_FROM_ABI float legendref(unsigned __n, float __x) { return std::legendre(__n, __x); }
119+
120+
inline _LIBCPP_HIDE_FROM_ABI long double legendrel(unsigned __n, long double __x) { return std::legendre(__n, __x); }
121+
122+
template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
123+
_LIBCPP_HIDE_FROM_ABI double legendre(unsigned __n, _Integer __x) {
124+
return std::legendre(__n, static_cast<double>(__x));
125+
}
126+
80127
#endif // _LIBCPP_STD_VER >= 17
81128

82129
_LIBCPP_END_NAMESPACE_STD

‎libcxx/include/cmath

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,14 @@ int ilogb (arithmetic x);
224224
int ilogbf(float x);
225225
int ilogbl(long double x);
226226
227+
double legendre(unsigned n, double x); // C++17
228+
float legendre(unsigned n, float x); // C++17
229+
long double legendre(unsigned n, long double x); // C++17
230+
float legendref(unsigned n, float x); // C++17
231+
long double legendrel(unsigned n, long double x); // C++17
232+
template <class Integer>
233+
double legendre(unsigned n, Integer x); // C++17
234+
227235
floating_point lgamma (arithmetic x);
228236
float lgammaf(float x);
229237
long double lgammal(long double x);

‎libcxx/modules/std/cmath.inc

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -346,12 +346,14 @@ export namespace std {
346346
using std::laguerre;
347347
using std::laguerref;
348348
using std::laguerrel;
349+
#endif
349350

350351
// [sf.cmath.legendre], Legendre polynomials
351352
using std::legendre;
352353
using std::legendref;
353354
using std::legendrel;
354355

356+
#if 0
355357
// [sf.cmath.riemann.zeta], Riemann zeta function
356358
using std::riemann_zeta;
357359
using std::riemann_zetaf;

‎libcxx/test/std/numerics/c.math/legendre.pass.cpp

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
//===----------------------------------------------------------------------===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
// UNSUPPORTED: c++03, c++11, c++14
10+
11+
// <cmath>
12+
13+
// double legendre(unsigned n, double x);
14+
// float legendre(unsigned n, float x);
15+
// long double legendre(unsigned n, long double x);
16+
// float legendref(unsigned n, float x);
17+
// long double legendrel(unsigned n, long double x);
18+
// template <class Integer>
19+
// double legendre(unsigned n, Integer x);
20+
21+
#include <cassert>
22+
#include <cmath>
23+
#include <limits>
24+
25+
template <class T>
26+
void testLegendreNaNPropagation() {
27+
const unsigned MaxN = 127;
28+
const T x = std::numeric_limits<T>::quiet_NaN();
29+
for (unsigned n = 0; n <= MaxN; ++n) {
30+
assert(std::isnan(std::legendre(n, x)));
31+
}
32+
}
33+
34+
template <class T>
35+
void testLegendreNotNaN(const T x) {
36+
assert(!std::isnan(x));
37+
const unsigned MaxN = 127;
38+
for (unsigned n = 0; n <= MaxN; ++n) {
39+
assert(!std::isnan(std::legendre(n, x)));
40+
}
41+
}
42+
43+
template <class T>
44+
void testLegendreThrows(const T x) {
45+
#ifndef _LIBCPP_NO_EXCEPTIONS
46+
const unsigned MaxN = 127;
47+
for (unsigned n = 0; n <= MaxN; ++n) {
48+
bool Throws = false;
49+
try {
50+
std::legendre(n, x);
51+
} catch (const std::domain_error&) {
52+
Throws = true;
53+
}
54+
assert(Throws);
55+
}
56+
#endif // _LIBCPP_NO_EXCEPTIONS
57+
}
58+
59+
template <class T>
60+
void testLegendreAnalytic(const T x, const T AbsTolerance, const T RelTolerance) {
61+
assert(!std::isnan(x));
62+
const auto compareFloatingPoint = [AbsTolerance, RelTolerance](const T Result, const T ExpectedResult) {
63+
if (std::isinf(ExpectedResult) && std::isinf(Result))
64+
return true;
65+
66+
if (std::isnan(ExpectedResult) || std::isnan(Result))
67+
return false;
68+
69+
const T Tolerance = AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
70+
return std::abs(Result - ExpectedResult) < Tolerance;
71+
};
72+
73+
const auto l0 = [](T) { return T(1); };
74+
const auto l1 = [](T y) { return y; };
75+
const auto l2 = [](T y) { return (T(3) * y * y - T(1)) / T(2); };
76+
const auto l3 = [](T y) { return (T(5) * y * y - T(3)) * y / T(2); };
77+
const auto l4 = [](T y) { return (T(35) * y * y * y * y - T(30) * y * y + T(3)) / T(8); };
78+
const auto l5 = [](T y) { return (T(63) * y * y * y * y - T(70) * y * y + T(15)) * y / T(8); };
79+
const auto l6 = [](T y) {
80+
const T y2 = y * y;
81+
return (T(231) * y2 * y2 * y2 - T(315) * y2 * y2 + T(105) * y2 - T(5)) / T(16);
82+
};
83+
84+
assert(compareFloatingPoint(std::legendre(0, x), l0(x)));
85+
assert(compareFloatingPoint(std::legendre(1, x), l1(x)));
86+
assert(compareFloatingPoint(std::legendre(2, x), l2(x)));
87+
assert(compareFloatingPoint(std::legendre(3, x), l3(x)));
88+
assert(compareFloatingPoint(std::legendre(4, x), l4(x)));
89+
assert(compareFloatingPoint(std::legendre(5, x), l5(x)));
90+
assert(compareFloatingPoint(std::legendre(6, x), l6(x)));
91+
}
92+
93+
template <class T>
94+
void testLegendre(const T AbsTolerance, const T RelTolerance) {
95+
testLegendreNaNPropagation<T>();
96+
testLegendreThrows<T>(T(-5));
97+
testLegendreThrows<T>(T(5));
98+
99+
const T Samples[] = {T(-1.0), T(-0.5), T(-0.1), T(0.0), T(0.1), T(0.5), T(1.0)};
100+
101+
for (T x : Samples) {
102+
testLegendreNotNaN(x);
103+
testLegendreAnalytic(x, AbsTolerance, RelTolerance);
104+
}
105+
}
106+
107+
int main() {
108+
testLegendre<float>(1e-6f, 1e-6f);
109+
testLegendre<double>(1e-9, 1e-9);
110+
testLegendre<long double>(1e-12, 1e-12);
111+
}

0 commit comments

Comments
 (0)