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 hermite(unsigned n, double x); 14 // float hermite(unsigned n, float x); 15 // long double hermite(unsigned n, long double x); 16 // float hermitef(unsigned n, float x); 17 // long double hermitel(unsigned n, long double x); 18 // template <class Integer> 19 // double hermite(unsigned n, Integer x); 20 21 #include <array> 22 #include <cassert> 23 #include <cmath> 24 #include <limits> 25 #include <vector> 26 27 #include "type_algorithms.h" 28 29 template <class Real> 30 constexpr unsigned get_maximal_order() { 31 if constexpr (std::numeric_limits<Real>::is_iec559) 32 return 128; 33 else { // Workaround for z/OS HexFloat. 34 // Note |H_n(x)| < 10^75 for n < 39 and x in sample_points(). 35 static_assert(std::numeric_limits<Real>::max_exponent10 == 75); 36 return 39; 37 } 38 } 39 40 template <class T> 41 std::array<T, 11> sample_points() { 42 return {-12.34, -7.42, -1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0, 5.67, 15.67}; 43 } 44 45 template <class Real> 46 class CompareFloatingValues { 47 private: 48 Real abs_tol; 49 Real rel_tol; 50 51 public: 52 CompareFloatingValues() { 53 abs_tol = []() -> Real { 54 if (std::is_same_v<Real, float>) 55 return 1e-5f; 56 else if (std::is_same_v<Real, double>) 57 return 1e-11; 58 else 59 return 1e-12l; 60 }(); 61 62 rel_tol = abs_tol; 63 } 64 65 bool operator()(Real result, Real expected) const { 66 if (std::isinf(expected) && std::isinf(result)) 67 return result == expected; 68 69 if (std::isnan(expected) || std::isnan(result)) 70 return false; 71 72 Real tol = abs_tol + std::abs(expected) * rel_tol; 73 return std::abs(result - expected) < tol; 74 } 75 }; 76 77 // Roots are taken from 78 // Salzer, Herbert E., Ruth Zucker, and Ruth Capuano. 79 // Table of the zeros and weight factors of the first twenty Hermite 80 // polynomials. US Government Printing Office, 1952. 81 template <class T> 82 std::vector<T> get_roots(unsigned n) { 83 switch (n) { 84 case 0: 85 return {}; 86 case 1: 87 return {T(0)}; 88 case 2: 89 return {T(0.707106781186548)}; 90 case 3: 91 return {T(0), T(1.224744871391589)}; 92 case 4: 93 return {T(0.524647623275290), T(1.650680123885785)}; 94 case 5: 95 return {T(0), T(0.958572464613819), T(2.020182870456086)}; 96 case 6: 97 return {T(0.436077411927617), T(1.335849074013697), T(2.350604973674492)}; 98 case 7: 99 return {T(0), T(0.816287882858965), T(1.673551628767471), T(2.651961356835233)}; 100 case 8: 101 return {T(0.381186990207322), T(1.157193712446780), T(1.981656756695843), T(2.930637420257244)}; 102 case 9: 103 return {T(0), T(0.723551018752838), T(1.468553289216668), T(2.266580584531843), T(3.190993201781528)}; 104 case 10: 105 return { 106 T(0.342901327223705), T(1.036610829789514), T(1.756683649299882), T(2.532731674232790), T(3.436159118837738)}; 107 case 11: 108 return {T(0), 109 T(0.65680956682100), 110 T(1.326557084494933), 111 T(2.025948015825755), 112 T(2.783290099781652), 113 T(3.668470846559583)}; 114 115 case 12: 116 return {T(0.314240376254359), 117 T(0.947788391240164), 118 T(1.597682635152605), 119 T(2.279507080501060), 120 T(3.020637025120890), 121 T(3.889724897869782)}; 122 123 case 13: 124 return {T(0), 125 T(0.605763879171060), 126 T(1.220055036590748), 127 T(1.853107651601512), 128 T(2.519735685678238), 129 T(3.246608978372410), 130 T(4.101337596178640)}; 131 132 case 14: 133 return {T(0.29174551067256), 134 T(0.87871378732940), 135 T(1.47668273114114), 136 T(2.09518325850772), 137 T(2.74847072498540), 138 T(3.46265693360227), 139 T(4.30444857047363)}; 140 141 case 15: 142 return {T(0.00000000000000), 143 T(0.56506958325558), 144 T(1.13611558521092), 145 T(1.71999257518649), 146 T(2.32573248617386), 147 T(2.96716692790560), 148 T(3.66995037340445), 149 T(4.49999070730939)}; 150 151 case 16: 152 return {T(0.27348104613815), 153 T(0.82295144914466), 154 T(1.38025853919888), 155 T(1.95178799091625), 156 T(2.54620215784748), 157 T(3.17699916197996), 158 T(3.86944790486012), 159 T(4.68873893930582)}; 160 161 case 17: 162 return {T(0), 163 T(0.5316330013427), 164 T(1.0676487257435), 165 T(1.6129243142212), 166 T(2.1735028266666), 167 T(2.7577629157039), 168 T(3.3789320911415), 169 T(4.0619466758755), 170 T(4.8713451936744)}; 171 172 case 18: 173 return {T(0.2582677505191), 174 T(0.7766829192674), 175 T(1.3009208583896), 176 T(1.8355316042616), 177 T(2.3862990891667), 178 T(2.9613775055316), 179 T(3.5737690684863), 180 T(4.2481178735681), 181 T(5.0483640088745)}; 182 183 case 19: 184 return {T(0), 185 T(0.5035201634239), 186 T(1.0103683871343), 187 T(1.5241706193935), 188 T(2.0492317098506), 189 T(2.5911337897945), 190 T(3.1578488183476), 191 T(3.7621873519640), 192 T(4.4285328066038), 193 T(5.2202716905375)}; 194 195 case 20: 196 return {T(0.2453407083009), 197 T(0.7374737285454), 198 T(1.2340762153953), 199 T(1.7385377121166), 200 T(2.2549740020893), 201 T(2.7888060584281), 202 T(3.347854567332), 203 T(3.9447640401156), 204 T(4.6036824495507), 205 T(5.3874808900112)}; 206 207 default: // polynom degree n>20 is unsupported 208 assert(false); 209 return {T(-42)}; 210 } 211 } 212 213 template <class Real> 214 void test() { 215 if constexpr ( 216 std::numeric_limits<Real>::has_quiet_NaN && 217 std::numeric_limits< 218 Real>::has_signaling_NaN) { // checks if NaNs are reported correctly (i.e. output == input for input == NaN) 219 using nl = std::numeric_limits<Real>; 220 for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()}) 221 for (unsigned n = 0; n < get_maximal_order<Real>(); ++n) 222 assert(std::isnan(std::hermite(n, NaN))); 223 } 224 225 if constexpr (std::numeric_limits<Real>::has_quiet_NaN && 226 std::numeric_limits< 227 Real>::has_signaling_NaN) { // simple sample points for n=0..127 should not produce NaNs. 228 for (Real x : sample_points<Real>()) 229 for (unsigned n = 0; n < get_maximal_order<Real>(); ++n) 230 assert(!std::isnan(std::hermite(n, x))); 231 } 232 233 { // checks std::hermite(n, x) for n=0..5 against analytic polynoms 234 const auto h0 = [](Real) -> Real { return 1; }; 235 const auto h1 = [](Real y) -> Real { return 2 * y; }; 236 const auto h2 = [](Real y) -> Real { return 4 * y * y - 2; }; 237 const auto h3 = [](Real y) -> Real { return y * (8 * y * y - 12); }; 238 const auto h4 = [](Real y) -> Real { return (16 * std::pow(y, 4) - 48 * y * y + 12); }; 239 const auto h5 = [](Real y) -> Real { return y * (32 * std::pow(y, 4) - 160 * y * y + 120); }; 240 241 for (Real x : sample_points<Real>()) { 242 const CompareFloatingValues<Real> compare; 243 assert(compare(std::hermite(0, x), h0(x))); 244 assert(compare(std::hermite(1, x), h1(x))); 245 assert(compare(std::hermite(2, x), h2(x))); 246 assert(compare(std::hermite(3, x), h3(x))); 247 assert(compare(std::hermite(4, x), h4(x))); 248 assert(compare(std::hermite(5, x), h5(x))); 249 } 250 } 251 252 { // checks std::hermitef for bitwise equality with std::hermite(unsigned, float) 253 if constexpr (std::is_same_v<Real, float>) 254 for (unsigned n = 0; n < get_maximal_order<Real>(); ++n) 255 for (float x : sample_points<float>()) 256 assert(std::hermite(n, x) == std::hermitef(n, x)); 257 } 258 259 { // checks std::hermitel for bitwise equality with std::hermite(unsigned, long double) 260 if constexpr (std::is_same_v<Real, long double>) 261 for (unsigned n = 0; n < get_maximal_order<Real>(); ++n) 262 for (long double x : sample_points<long double>()) 263 assert(std::hermite(n, x) == std::hermitel(n, x)); 264 } 265 266 { // Checks if the characteristic recurrence relation holds: H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x) 267 for (Real x : sample_points<Real>()) { 268 for (unsigned n = 1; n < get_maximal_order<Real>() - 1; ++n) { 269 Real H_next = std::hermite(n + 1, x); 270 Real H_next_recurrence = 2 * (x * std::hermite(n, x) - n * std::hermite(n - 1, x)); 271 272 if (std::isinf(H_next)) 273 break; 274 const CompareFloatingValues<Real> compare; 275 assert(compare(H_next, H_next_recurrence)); 276 } 277 } 278 } 279 280 { // sanity checks: hermite polynoms need to change signs at (simple) roots. checked upto order n<=20. 281 282 // root tolerance: must be smaller than the smallest difference between adjacent roots 283 Real tol = []() -> Real { 284 if (std::is_same_v<Real, float>) 285 return 1e-5f; 286 else if (std::is_same_v<Real, double>) 287 return 1e-9; 288 else 289 return 1e-10l; 290 }(); 291 292 const auto is_sign_change = [tol](unsigned n, Real x) -> bool { 293 return std::hermite(n, x - tol) * std::hermite(n, x + tol) < 0; 294 }; 295 296 for (unsigned n = 0; n <= 20u; ++n) { 297 for (Real x : get_roots<Real>(n)) { 298 // the roots are symmetric: if x is a root, so is -x 299 if (x > 0) 300 assert(is_sign_change(n, -x)); 301 assert(is_sign_change(n, x)); 302 } 303 } 304 } 305 306 if constexpr (std::numeric_limits<Real>::has_infinity) { // check input infinity is handled correctly 307 Real inf = std::numeric_limits<Real>::infinity(); 308 for (unsigned n = 1; n < get_maximal_order<Real>(); ++n) { 309 assert(std::hermite(n, +inf) == inf); 310 assert(std::hermite(n, -inf) == ((n & 1) ? -inf : inf)); 311 } 312 } 313 314 if constexpr (std::numeric_limits< 315 Real>::has_infinity) { // check: if overflow occurs that it is mapped to the correct infinity 316 if constexpr (std::is_same_v<Real, double>) { 317 // Q: Why only double? 318 // A: The numeric values (e.g. overflow threshold `n`) below are different for other types. 319 static_assert(sizeof(double) == 8); 320 for (unsigned n = 0; n < get_maximal_order<Real>(); ++n) { 321 // Q: Why n=111 and x=300? 322 // A: Both are chosen s.t. the first overlow occurs for some `n<get_maximal_order<Real>()`. 323 if (n < 111) { 324 assert(std::isfinite(std::hermite(n, +300.0))); 325 assert(std::isfinite(std::hermite(n, -300.0))); 326 } else { 327 double inf = std::numeric_limits<double>::infinity(); 328 assert(std::hermite(n, +300.0) == inf); 329 assert(std::hermite(n, -300.0) == ((n & 1) ? -inf : inf)); 330 } 331 } 332 } 333 } 334 } 335 336 struct TestFloat { 337 template <class Real> 338 void operator()() { 339 test<Real>(); 340 } 341 }; 342 343 struct TestInt { 344 template <class Integer> 345 void operator()() { 346 // checks that std::hermite(unsigned, Integer) actually wraps std::hermite(unsigned, double) 347 for (unsigned n = 0; n < get_maximal_order<double>(); ++n) 348 for (Integer x : {-42, -7, -5, -1, 0, 1, 5, 7, 42}) 349 assert(std::hermite(n, x) == std::hermite(n, static_cast<double>(x))); 350 } 351 }; 352 353 int main() { 354 types::for_each(types::floating_point_types(), TestFloat()); 355 types::for_each(types::type_list<short, int, long, long long>(), TestInt()); 356 } 357