|
| 1 | +// transcendentals_highprecision.cpp: identity-driven high-precision validation of |
| 2 | +// the elreal transcendentals (exp/log/pow, sin/cos/tan, atan/asin/acos, |
| 3 | +// sinh/cosh/tanh) to >= 300 decimal digits (#1049). |
| 4 | +// |
| 5 | +// Why this exists |
| 6 | +// --------------- |
| 7 | +// The per-PR transcendental tests (exponent.cpp, hyperbolic.cpp, trigonometry.cpp) |
| 8 | +// only validate to host-double tolerance (~1e-10, ~16 digits). For an exact-lazy |
| 9 | +// real type that is almost no coverage: the argument-reduction + series + mul/div |
| 10 | +// machinery these compose from is untested past ~16 digits, so a high-precision |
| 11 | +// operator bug is invisible (#1044 already surfaced two under cancellation). 300-digit |
| 12 | +// external transcendental tables mostly do not exist, so the primary mechanism here |
| 13 | +// is SELF-CHECKING IDENTITIES validated exactly via the independent dyadic oracle |
| 14 | +// (elreal_reference_digits.hpp, per #987): an identity residual is converted to an |
| 15 | +// exact dyadic and compared with no shared code path through elreal arithmetic, so |
| 16 | +// an operator bug cannot hide inside a loose tolerance. A handful of absolute |
| 17 | +// spot-checks against mpmath 320-digit strings (reference_constants.hpp) catch a |
| 18 | +// wrong-but-self-consistent series that an identity alone would mask. |
| 19 | +// |
| 20 | +// Cost / placement |
| 21 | +// ---------------- |
| 22 | +// The 300-digit checks run the functions at depth ~20 (O(depth^4), ~8000x a depth-2 |
| 23 | +// run). They are gated to REGRESSION_LEVEL_4 (stress / manual) and are a no-op at the |
| 24 | +// sanity level CI / sanitizers / coverage use. A small always-on fast-tier identity |
| 25 | +// section (depth 6, double tolerance) provides per-PR coverage. |
| 26 | +// |
| 27 | +// Copyright (C) 2017 Stillwater Supercomputing, Inc. |
| 28 | +// SPDX-License-Identifier: MIT |
| 29 | +// |
| 30 | +// This file is part of the universal numbers project, which is released under an MIT Open Source license. |
| 31 | +#include <universal/utility/directives.hpp> |
| 32 | + |
| 33 | +// Regression testing guards: typically set by the cmake configuration, but |
| 34 | +// MANUAL_TESTING is an override. A bare manual compile (no -D) runs everything. |
| 35 | +#define MANUAL_TESTING 0 |
| 36 | +#ifndef REGRESSION_LEVEL_OVERRIDE |
| 37 | +#undef REGRESSION_LEVEL_1 |
| 38 | +#undef REGRESSION_LEVEL_2 |
| 39 | +#undef REGRESSION_LEVEL_3 |
| 40 | +#undef REGRESSION_LEVEL_4 |
| 41 | +#define REGRESSION_LEVEL_1 1 |
| 42 | +#define REGRESSION_LEVEL_2 1 |
| 43 | +#define REGRESSION_LEVEL_3 1 |
| 44 | +#define REGRESSION_LEVEL_4 1 |
| 45 | +#endif |
| 46 | + |
| 47 | +// asin(x)+acos(x) == pi/2 exercises add() over an EXACT leading-term cancellation |
| 48 | +// (acos = pi/2 - asin), which add() does not renormalize to 0-overlap -- it aborts |
| 49 | +// the ZBCL invariant in assert-enabled builds (#1057). The value is still correct; |
| 50 | +// only the canonical form is violated. Gated off until #1057 lands; flip to 1 then. |
| 51 | +#define ELREAL_ADD_EXACT_CANCEL_FIXED 0 |
| 52 | + |
| 53 | +// exp / sin / cos (and tan / sinh / cosh / tanh built on them) cap at host-double |
| 54 | +// precision because their Taylor/Maclaurin series scale each term by a host-double |
| 55 | +// reciprocal coefficient (#1058, same bug class as e_zbcl #1054). Their 300-digit |
| 56 | +// checks are gated off until #1058 lands; flip to 1 then. log/atan/asin/acos use the |
| 57 | +// full-precision odd_power_series and already reach 300+, so they stay active. |
| 58 | +#define ELREAL_EXP_SERIES_HIGH_PRECISION 0 |
| 59 | + |
| 60 | +#include <cmath> |
| 61 | +#include <iostream> |
| 62 | +#include <string> |
| 63 | +#include <string_view> |
| 64 | + |
| 65 | +#include <universal/number/elreal/elreal.hpp> |
| 66 | +#include <universal/number/elreal/math/constants.hpp> |
| 67 | +#include <universal/number/elreal/math/exponent.hpp> |
| 68 | +#include <universal/number/elreal/math/hyperbolic.hpp> |
| 69 | +#include <universal/number/elreal/math/trigonometry.hpp> |
| 70 | +#include <universal/verification/elreal_oracle.hpp> |
| 71 | +#include <universal/verification/elreal_reference_digits.hpp> |
| 72 | +#include <universal/verification/test_suite.hpp> |
| 73 | +#include <math/constants/reference_constants.hpp> |
| 74 | + |
| 75 | +namespace { |
| 76 | + |
| 77 | +using sw::universal::ZBCL; |
| 78 | +using ER = ZBCL<double>; |
| 79 | + |
| 80 | +// exact subtraction a - b via the elreal EFT add (no dedicated sub). |
| 81 | +ER sub(const ER& a, const ER& b) { return sw::universal::add(a, sw::universal::negate(b)); } |
| 82 | + |
| 83 | +// ---- fast tier (always on): identities to host-double tolerance ---------------- |
| 84 | +// Mirrors the per-PR coverage in exponent/hyperbolic/trigonometry.cpp so the new |
| 85 | +// file still exercises the identities when REGRESSION_LEVEL_4 is off (CI). Uses the |
| 86 | +// double approximation of each result (elreal_oracle::approx). |
| 87 | +int fast_identities(bool reportTestCases) { |
| 88 | + using namespace sw::universal; |
| 89 | + namespace est = sw::universal::elreal_oracle; |
| 90 | + constexpr std::size_t D = 6; |
| 91 | + const double tol = 1e-12; |
| 92 | + const double pi = std::acos(-1.0); |
| 93 | + int n = 0; |
| 94 | + |
| 95 | + auto close = [&](const char* name, double got, double want) { |
| 96 | + if (std::abs(got - want) > tol) { std::cout << " FAIL fast " << name << ": " << got << " != " << want << '\n'; ++n; } |
| 97 | + else if (reportTestCases) { std::cout << " ok fast " << name << '\n'; } |
| 98 | + }; |
| 99 | + |
| 100 | + ER half = from_native<double>(0.5); |
| 101 | + close("log(exp(0.5))==0.5", est::approx(log(exp(half, D), D)), 0.5); |
| 102 | + close("exp(log(2))==2", est::approx(exp(log(from_native<double>(2.0), D), D)), 2.0); |
| 103 | + |
| 104 | + ER s = sin(half, D), c = cos(half, D); |
| 105 | + close("sin^2+cos^2==1", est::approx(add(mul(s, s, D), mul(c, c, D))), 1.0); |
| 106 | + |
| 107 | + ER sh = sinh(half, D), ch = cosh(half, D); |
| 108 | + close("cosh^2-sinh^2==1", est::approx(sub(mul(ch, ch, D), mul(sh, sh, D))), 1.0); |
| 109 | + |
| 110 | +#if ELREAL_ADD_EXACT_CANCEL_FIXED |
| 111 | + close("asin+acos==pi/2", est::approx(add(asin(half, D), acos(half, D))), pi / 2.0); |
| 112 | +#endif |
| 113 | + close("4*atan(1)==pi", 4.0 * est::approx(atan(from_native<double>(1.0), D)), pi); |
| 114 | + |
| 115 | + // addition theorem exp(a+b) == exp(a)*exp(b) |
| 116 | + ER a = from_native<double>(0.5), b = from_native<double>(0.25); |
| 117 | + close("exp(a+b)==exp(a)exp(b)", |
| 118 | + est::approx(exp(add(a, b), D)) - est::approx(mul(exp(a, D), exp(b, D), D)) + 1.0, 1.0); |
| 119 | + |
| 120 | + // 0-overlap invariant on a representative result |
| 121 | + n += est::check_zero_overlap(sin(half, D), D, "fast sin(0.5)"); |
| 122 | + return n; |
| 123 | +} |
| 124 | + |
| 125 | +#if REGRESSION_LEVEL_4 |
| 126 | +using sw::universal::agreed_decimal_digits; |
| 127 | + |
| 128 | +constexpr std::size_t kDepth = 20; // function-evaluation depth (double-host ceiling) |
| 129 | +constexpr std::size_t kWork = 26; // working depth for identity mul/div (kDepth + guard) |
| 130 | +constexpr int kMinDigits = 300; |
| 131 | + |
| 132 | +// Absolute spot-check: z agrees with the decimal reference to >= kMinDigits digits. |
| 133 | +int check_value(const char* name, const ER& z, std::string_view ref, bool reportTestCases) { |
| 134 | + int got = agreed_decimal_digits(z, ref); |
| 135 | + if (got < kMinDigits) { |
| 136 | + std::cout << " FAIL " << name << ": agreed " << got << " digits, want >= " << kMinDigits << '\n'; |
| 137 | + return 1; |
| 138 | + } |
| 139 | + if (reportTestCases) std::cout << " ok " << name << ": " << got << " digits\n"; |
| 140 | + return 0; |
| 141 | +} |
| 142 | + |
| 143 | +// Identity check: lhs agrees with rhs (relative) to >= kMinDigits digits. Neither |
| 144 | +// side needs a closed-form decimal string; both are compared as exact dyadics. |
| 145 | +int check_identity(const char* name, const ER& lhs, const ER& rhs, bool reportTestCases) { |
| 146 | + int got = agreed_decimal_digits(lhs, rhs); |
| 147 | + if (got < kMinDigits) { |
| 148 | + std::cout << " FAIL " << name << ": identity holds to only " << got << " digits, want >= " << kMinDigits << '\n'; |
| 149 | + return 1; |
| 150 | + } |
| 151 | + if (reportTestCases) std::cout << " ok " << name << ": " << got << " digits\n"; |
| 152 | + return 0; |
| 153 | +} |
| 154 | +#endif |
| 155 | + |
| 156 | +} // anonymous |
| 157 | + |
| 158 | +int main() |
| 159 | +try { |
| 160 | + using namespace sw::universal; |
| 161 | + std::string test_suite = "elreal high-precision transcendental hardening (identity-driven, #1049)"; |
| 162 | + int nrOfFailedTestCases = 0; |
| 163 | + bool reportTestCases = true; |
| 164 | + ReportTestSuiteHeader(test_suite, reportTestCases); |
| 165 | + |
| 166 | + nrOfFailedTestCases += fast_identities(reportTestCases); |
| 167 | + |
| 168 | +#if REGRESSION_LEVEL_4 |
| 169 | + const std::size_t D = kDepth; |
| 170 | + const std::size_t W = kWork; |
| 171 | + ER half = from_native<double>(0.5); |
| 172 | + |
| 173 | + // === full-precision group (log / atan / asin / acos): reaches 300+ digits ===== |
| 174 | + // These compose div/mul/add/sum/odd_power_series, whose series divides by the |
| 175 | + // integer denominator at FULL precision, so they hit the double-host ceiling. |
| 176 | + // Absolute spot-checks against mpmath 320-digit references; arguments are exact |
| 177 | + // doubles (1, 1/2) so the reference equals what elreal sees. |
| 178 | + nrOfFailedTestCases += check_value("log(2)", log(from_native<double>(2.0), D), s_ln2, reportTestCases); |
| 179 | + nrOfFailedTestCases += check_value("atan(1)", atan(from_native<double>(1.0), D), s_pi_4, reportTestCases); |
| 180 | + nrOfFailedTestCases += check_value("asin(0.5)", asin(from_native<double>(0.5), D), s_pi_6, reportTestCases); |
| 181 | + nrOfFailedTestCases += check_value("acos(0.5)", acos(from_native<double>(0.5), D), s_pi_3, reportTestCases); |
| 182 | + // identities within the full-precision group |
| 183 | + nrOfFailedTestCases += check_value("atan(0.5)+atan(2)==pi/2", |
| 184 | + add(atan(half, D), atan(from_native<double>(2.0), D)), s_pi_2, reportTestCases); |
| 185 | + nrOfFailedTestCases += check_value("4*atan(1)==pi", |
| 186 | + mul(atan(from_native<double>(1.0), D), from_native<double>(4.0), W), s_pi, reportTestCases); |
| 187 | + ER x = from_native<double>(2.0), y = from_native<double>(3.0); |
| 188 | + nrOfFailedTestCases += check_identity("log(xy)==log(x)+log(y)", |
| 189 | + log(mul(x, y, W), D), add(log(x, D), log(y, D)), reportTestCases); |
| 190 | + |
| 191 | + // asin(x)+acos(x)==pi/2 exercises add() over an EXACT leading-term cancellation; |
| 192 | + // add() does not renormalize it to 0-overlap (aborts the ZBCL invariant). #1057. |
| 193 | +#if ELREAL_ADD_EXACT_CANCEL_FIXED |
| 194 | + nrOfFailedTestCases += check_value("asin(0.5)+acos(0.5)==pi/2", |
| 195 | + add(asin(half, D), acos(half, D)), s_pi_2, reportTestCases); |
| 196 | +#endif |
| 197 | + |
| 198 | + // === exp / sin / cos family: currently capped at ~17 digits (#1058) =========== |
| 199 | + // exp(), sin_series/cos_series multiply each Taylor term by a HOST-DOUBLE |
| 200 | + // reciprocal coefficient (1/n, 1/(2n+1)!), capping the sum at host-double |
| 201 | + // precision; sinh/cosh/tanh/tan compose them. Same bug class as e_zbcl (#1054) |
| 202 | + // but in the FUNCTIONS, not the constant. Gated off until #1058 lands; flip to 1 |
| 203 | + // then -- the references and identities below already pass at the fixed precision. |
| 204 | +#if ELREAL_EXP_SERIES_HIGH_PRECISION |
| 205 | + ER a = from_native<double>(0.5), b = from_native<double>(0.25); |
| 206 | + // absolute spot-checks |
| 207 | + nrOfFailedTestCases += check_value("exp(1)", exp(from_native<double>(1.0), D), s_e, reportTestCases); |
| 208 | + nrOfFailedTestCases += check_value("sin(0.5)", sin(from_native<double>(0.5), D), s_sin_half, reportTestCases); |
| 209 | + nrOfFailedTestCases += check_value("cos(0.5)", cos(from_native<double>(0.5), D), s_cos_half, reportTestCases); |
| 210 | + nrOfFailedTestCases += check_value("tan(0.5)", tan(from_native<double>(0.5), D), s_tan_half, reportTestCases); |
| 211 | + nrOfFailedTestCases += check_value("sinh(0.5)", sinh(from_native<double>(0.5), D), s_sinh_half, reportTestCases); |
| 212 | + nrOfFailedTestCases += check_value("cosh(0.5)", cosh(from_native<double>(0.5), D), s_cosh_half, reportTestCases); |
| 213 | + nrOfFailedTestCases += check_value("tanh(0.5)", tanh(from_native<double>(0.5), D), s_tanh_half, reportTestCases); |
| 214 | + // inverse round-trips |
| 215 | + nrOfFailedTestCases += check_value("log(exp(0.5))==0.5", log(exp(half, D), D), "0.5", reportTestCases); |
| 216 | + nrOfFailedTestCases += check_value("exp(log(2))==2", exp(log(from_native<double>(2.0), D), D), "2", reportTestCases); |
| 217 | + nrOfFailedTestCases += check_value("sin(asin(0.5))==0.5", sin(asin(half, D), D), "0.5", reportTestCases); |
| 218 | + nrOfFailedTestCases += check_value("cos(acos(0.5))==0.5", cos(acos(half, D), D), "0.5", reportTestCases); |
| 219 | + nrOfFailedTestCases += check_value("tan(atan(0.5))==0.5", tan(atan(half, D), D), "0.5", reportTestCases); |
| 220 | + // Pythagorean / hyperbolic |
| 221 | + ER s = sin(half, D), c = cos(half, D); |
| 222 | + nrOfFailedTestCases += check_value("sin^2+cos^2==1", add(mul(s, s, W), mul(c, c, W)), "1", reportTestCases); |
| 223 | + ER sh = sinh(half, D), ch = cosh(half, D); |
| 224 | + nrOfFailedTestCases += check_value("cosh^2-sinh^2==1", sub(mul(ch, ch, W), mul(sh, sh, W)), "1", reportTestCases); |
| 225 | + // addition theorems (compared side-to-side as exact dyadics) |
| 226 | + nrOfFailedTestCases += check_identity("exp(a+b)==exp(a)*exp(b)", |
| 227 | + exp(add(a, b), D), mul(exp(a, D), exp(b, D), W), reportTestCases); |
| 228 | + ER sa = sin(a, D), ca = cos(a, D), sb = sin(b, D), cb = cos(b, D); |
| 229 | + nrOfFailedTestCases += check_identity("sin(a+b)==sin a cos b+cos a sin b", |
| 230 | + sin(add(a, b), D), add(mul(sa, cb, W), mul(ca, sb, W)), reportTestCases); |
| 231 | + ER sha = sinh(a, D), cha = cosh(a, D), shb = sinh(b, D), chb = cosh(b, D); |
| 232 | + nrOfFailedTestCases += check_identity("cosh(a+b)==cosh a cosh b+sinh a sinh b", |
| 233 | + cosh(add(a, b), D), add(mul(cha, chb, W), mul(sha, shb, W)), reportTestCases); |
| 234 | +#endif |
| 235 | +#else |
| 236 | + std::cout << " high-precision transcendental validation runs at REGRESSION_LEVEL_4 (stress) only\n"; |
| 237 | +#endif |
| 238 | + |
| 239 | + ReportTestSuiteResults(test_suite, nrOfFailedTestCases); |
| 240 | + return (nrOfFailedTestCases > 0 ? EXIT_FAILURE : EXIT_SUCCESS); |
| 241 | +} |
| 242 | +catch (const std::exception& err) { |
| 243 | + std::cerr << "Caught unexpected exception: " << err.what() << std::endl; |
| 244 | + return EXIT_FAILURE; |
| 245 | +} |
0 commit comments