|
18 | 18 | // template <class Integer>
|
19 | 19 | // double legendre(unsigned n, Integer x);
|
20 | 20 |
|
| 21 | +#include <array> |
21 | 22 | #include <cassert>
|
22 | 23 | #include <cmath>
|
23 | 24 | #include <limits>
|
24 | 25 |
|
| 26 | +#include "type_algorithms.h" |
| 27 | + |
| 28 | +inline constexpr unsigned g_max_n = 128; |
| 29 | + |
25 | 30 | 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 |
| - } |
| 31 | +std::array<T, 7> sample_points() { |
| 32 | + return {-1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0}; |
32 | 33 | }
|
33 | 34 |
|
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))); |
| 35 | +template <class Real> |
| 36 | +class CompareFloatingValues { |
| 37 | +private: |
| 38 | + Real tol; |
| 39 | + |
| 40 | +public: |
| 41 | + CompareFloatingValues() { |
| 42 | + if (std::is_same_v<Real, float>) |
| 43 | + tol = 1e-6f; |
| 44 | + else if (std::is_same_v<Real, double>) |
| 45 | + tol = 1e-9; |
| 46 | + else |
| 47 | + tol = 1e-12l; |
40 | 48 | }
|
41 |
| -} |
42 | 49 |
|
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); |
| 50 | + bool operator()(Real result, Real expected) const { |
| 51 | + if (std::isinf(expected) && std::isinf(result)) |
| 52 | + return result == expected; |
| 53 | + |
| 54 | + if (std::isnan(expected) || std::isnan(result)) |
| 55 | + return false; |
| 56 | + |
| 57 | + return std::abs(result - expected) < (tol + std::abs(expected) * tol); |
| 58 | + } |
| 59 | +}; |
| 60 | + |
| 61 | +template <class Real> |
| 62 | +void test() { |
| 63 | + { // checks if NaNs are reported correctly (i.e. output == input for input == NaN) |
| 64 | + using nl = std::numeric_limits<Real>; |
| 65 | + for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()}) |
| 66 | + for (unsigned n = 0; n < g_max_n; ++n) |
| 67 | + assert(std::isnan(std::legendre(n, NaN))); |
55 | 68 | }
|
| 69 | + |
| 70 | + { // simple sample points for n=0..127 should not produce NaNs. |
| 71 | + for (Real x : sample_points<Real>()) |
| 72 | + for (unsigned n = 0; n < g_max_n; ++n) |
| 73 | + assert(!std::isnan(std::legendre(n, x))); |
| 74 | + } |
| 75 | + |
| 76 | + { // For x with abs(x) > 1 an domain_error exception should be thrown. |
| 77 | +#ifndef _LIBCPP_NO_EXCEPTIONS |
| 78 | + for (double absX : {2.0, 7.77, 42.42, std::numeric_limits<double>::infinity()}) |
| 79 | + for (Real x : {-absX, absX}) |
| 80 | + for (unsigned n = 0; n < g_max_n; ++n) { |
| 81 | + bool throws = false; |
| 82 | + try { |
| 83 | + std::legendre(n, x); |
| 84 | + } catch (const std::domain_error&) { |
| 85 | + throws = true; |
| 86 | + } |
| 87 | + assert(throws); |
| 88 | + } |
56 | 89 | #endif // _LIBCPP_NO_EXCEPTIONS
|
57 |
| -} |
| 90 | + } |
58 | 91 |
|
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; |
| 92 | + { // check against analytic polynoms for order n=0..6 |
| 93 | + for (Real x : sample_points<Real>()) { |
| 94 | + const auto l0 = [](Real) -> Real { return 1; }; |
| 95 | + const auto l1 = [](Real y) -> Real { return y; }; |
| 96 | + const auto l2 = [](Real y) -> Real { return (3 * y * y - 1) / 2; }; |
| 97 | + const auto l3 = [](Real y) -> Real { return (5 * y * y - 3) * y / 2; }; |
| 98 | + const auto l4 = [](Real y) -> Real { return (35 * std::pow(y, 4) - 30 * y * y + 3) / 8; }; |
| 99 | + const auto l5 = [](Real y) -> Real { return (63 * std::pow(y, 4) - 70 * y * y + 15) * y / 8; }; |
| 100 | + const auto l6 = [](Real y) -> Real { |
| 101 | + return (231 * std::pow(y, 6) - 315 * std::pow(y, 4) + 105 * y * y - 5) / 16; |
| 102 | + }; |
| 103 | + |
| 104 | + const CompareFloatingValues<Real> compare; |
| 105 | + assert(compare(std::legendre(0, x), l0(x))); |
| 106 | + assert(compare(std::legendre(1, x), l1(x))); |
| 107 | + assert(compare(std::legendre(2, x), l2(x))); |
| 108 | + assert(compare(std::legendre(3, x), l3(x))); |
| 109 | + assert(compare(std::legendre(4, x), l4(x))); |
| 110 | + assert(compare(std::legendre(5, x), l5(x))); |
| 111 | + assert(compare(std::legendre(6, x), l6(x))); |
| 112 | + } |
| 113 | + } |
65 | 114 |
|
66 |
| - if (std::isnan(ExpectedResult) || std::isnan(Result)) |
67 |
| - return false; |
| 115 | + { // checks std::legendref for bitwise equality with std::legendre(unsigned, float) |
| 116 | + if constexpr (std::is_same_v<Real, float>) |
| 117 | + for (unsigned n = 0; n < g_max_n; ++n) |
| 118 | + for (float x : sample_points<float>()) |
| 119 | + assert(std::legendre(n, x) == std::legendref(n, x)); |
| 120 | + } |
68 | 121 |
|
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 |
| -} |
| 122 | + { // checks std::legendrel for bitwise equality with std::legendre(unsigned, long double) |
| 123 | + if constexpr (std::is_same_v<Real, long double>) |
| 124 | + for (unsigned n = 0; n < g_max_n; ++n) |
| 125 | + for (long double x : sample_points<long double>()) { |
| 126 | + assert(std::legendre(n, x) == std::legendrel(n, x)); |
| 127 | + assert(std::legendre(n, x) == std::legendrel(n, x)); |
| 128 | + } |
| 129 | + } |
92 | 130 |
|
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)); |
| 131 | + { // evaluation at x=1: P_n(1) = 1 |
| 132 | + const CompareFloatingValues<Real> compare; |
| 133 | + for (unsigned n = 0; n < g_max_n; ++n) |
| 134 | + assert(compare(std::legendre(n, Real{1}), 1)); |
| 135 | + } |
98 | 136 |
|
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)}; |
| 137 | + { // evaluation at x=-1: P_n(-1) = (-1)^n |
| 138 | + const CompareFloatingValues<Real> compare; |
| 139 | + for (unsigned n = 0; n < g_max_n; ++n) |
| 140 | + assert(compare(std::legendre(n, Real{-1}), std::pow(-1, n))); |
| 141 | + } |
100 | 142 |
|
101 |
| - for (T x : Samples) { |
102 |
| - testLegendreNotNaN(x); |
103 |
| - testLegendreAnalytic(x, AbsTolerance, RelTolerance); |
| 143 | + { // evaluation at x=0: |
| 144 | + // P_{2n }(0) = (-1)^n (2n-1)!! / (2n)!! |
| 145 | + // P_{2n+1}(0) = 0 |
| 146 | + const CompareFloatingValues<Real> compare; |
| 147 | + Real doubleFactorials{1}; |
| 148 | + for (unsigned n = 0; n < g_max_n; ++n) { |
| 149 | + if (n & 1) // odd |
| 150 | + doubleFactorials *= n; |
| 151 | + else if (n != 0) // even and not zero |
| 152 | + doubleFactorials /= n; |
| 153 | + |
| 154 | + assert(compare(std::legendre(n, Real{0}), (n & 1) ? Real{0} : std::pow(-1, n / 2) * doubleFactorials)); |
| 155 | + } |
104 | 156 | }
|
105 | 157 | }
|
106 | 158 |
|
| 159 | +struct TestFloat { |
| 160 | + template <class Real> |
| 161 | + void operator()() { |
| 162 | + test<Real>(); |
| 163 | + } |
| 164 | +}; |
| 165 | + |
| 166 | +struct TestInt { |
| 167 | + template <class Integer> |
| 168 | + void operator()() { |
| 169 | + // checks that std::legendre(unsigned, Integer) actually wraps std::legendre(unsigned, double) |
| 170 | + for (unsigned n = 0; n < g_max_n; ++n) |
| 171 | + for (Integer x : {-1, 0, 1}) |
| 172 | + assert(std::legendre(n, x) == std::legendre(n, static_cast<double>(x))); |
| 173 | + } |
| 174 | +}; |
| 175 | + |
107 | 176 | int main() {
|
108 |
| - testLegendre<float>(1e-6f, 1e-6f); |
109 |
| - testLegendre<double>(1e-9, 1e-9); |
110 |
| - testLegendre<long double>(1e-12, 1e-12); |
| 177 | + types::for_each(types::floating_point_types(), TestFloat()); |
| 178 | + types::for_each(types::type_list<short, int, long, long long>(), TestInt()); |
111 | 179 | }
|
0 commit comments