libcxx/test/std/numerics/c.math/hermite.pass.cpp

   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 }