359 lines · 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-exception6//7//===----------------------------------------------------------------------===//8 9// UNSUPPORTED: c++03, c++11, c++1410 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 29template <class Real>30constexpr 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 40template <class T>41std::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 45template <class Real>46class CompareFloatingValues {47private:48 Real abs_tol;49 Real rel_tol;50 51public: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 else59 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 from78// Salzer, Herbert E., Ruth Zucker, and Ruth Capuano.79// Table of the zeros and weight factors of the first twenty Hermite80// polynomials. US Government Printing Office, 1952.81template <class T>82std::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 unsupported208 assert(false);209 return {T(-42)};210 }211}212 213template <class Real>214void 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 polynoms234 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 roots283 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 else289 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 -x299 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 correctly307 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 infinity316 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 overflow 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 336struct TestFloat {337 template <class Real>338 void operator()() {339 test<Real>();340 }341};342 343struct 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 353int main(int, char**) {354 types::for_each(types::floating_point_types(), TestFloat());355 types::for_each(types::type_list<short, int, long, long long>(), TestInt());356 357 return 0;358}359