From 0718545ae0ba92431780fb985a5f6478a1ae9a3d Mon Sep 17 00:00:00 2001 From: Ravenwater Date: Tue, 27 Aug 2024 18:56:39 -0400 Subject: [PATCH] adding a verification suite to pow() function for double-double --- include/universal/number/dd/dd_impl.hpp | 47 ++- include/universal/number/dd/math/classify.hpp | 10 +- static/dd/math/pow.cpp | 321 ++++++++++++++++-- 3 files changed, 355 insertions(+), 23 deletions(-) diff --git a/include/universal/number/dd/dd_impl.hpp b/include/universal/number/dd/dd_impl.hpp index 4ee30e73a..16cb7d93a 100644 --- a/include/universal/number/dd/dd_impl.hpp +++ b/include/universal/number/dd/dd_impl.hpp @@ -45,8 +45,10 @@ dd operator-(const dd&, const dd&); dd operator*(const dd&, const dd&); dd operator/(const dd&, const dd&); std::ostream& operator<<(std::ostream&, const dd&); +bool signbit(const dd&); dd pown(const dd&, int); dd frexp(const dd&, int*); +dd ldexp(const dd&, int); // dd is an unevaluated pair of IEEE-754 doubles that provides a (1,11,106) floating-point triple class dd { @@ -850,6 +852,17 @@ class dd { //////////////////////// precomputed constants of note ///////////////////////////////// +constexpr int dd_max_precision = 106; // in bits + +// simple constants +constexpr dd dd_thousand(1000.0, 0.0); +constexpr dd dd_hundred(100.0, 0.0); +constexpr dd dd_ten(10.0, 0.0); +constexpr dd dd_one(1.0, 0.0); +constexpr dd dd_half(0.5, 0.0); +constexpr dd dd_third(0.33333333333333331, 1.8503717077085941e-17); +constexpr dd dd_fourth(0.25, 0); + // precomputed double-double constants courtesy Scibuilders, Jack Poulson constexpr dd dd_2pi (6.283185307179586232e+00, 2.449293598294706414e-16); @@ -1231,7 +1244,7 @@ inline dd sqr(const dd& a) { p2 += 2.0 * a.high() * a.low(); p2 += a.low() * a.low(); - double s2, s1 = quick_two_sum(p1, p2, s2); + double s2{ 0 }, s1 = quick_two_sum(p1, p2, s2); return dd(s1, s2); } @@ -1256,6 +1269,35 @@ inline dd reciprocal(const dd& a) { } } +///////////////////////////////////////////////////////////////////////////// +// power functions + +dd cbrt(const dd& a) { + if (!isfinite(a) || a.iszero()) + return a; // NaN returns NaN; +/-Inf returns +/-Inf, +/-0.0 returns +/-0.0 + + bool signA = signbit(a); + int e; // 0.5 <= r < 1.0 + dd r = frexp(abs(a), &e); + while (e % 3 != 0) { + ++e; + r = ldexp(r, -1); + } + + // at this point, 0.125 <= r < 1.0 + dd x = pow(r.high(), -dd_third.high()); + + // refine estimate using Newton's iteration + x += x * (1.0 - r * sqr(x) * x) * dd_third; + x += x * (1.0 - r * sqr(x) * x) * dd_third; + x = reciprocal(x); + + if (signA) + x = -x; + + return ldexp(x, e / 3); +} + inline dd pown(const dd& a, int n) { if (a.isnan()) return a; @@ -1269,7 +1311,7 @@ inline dd pown(const dd& a, int n) { errno = EDOM; return dd(SpecificValue::qnan); } - return 1.0; + return dd(1.0); case 1: s = a; @@ -1277,6 +1319,7 @@ inline dd pown(const dd& a, int n) { case 2: s = sqr(a); + break; default: // Use binary exponentiation { diff --git a/include/universal/number/dd/math/classify.hpp b/include/universal/number/dd/math/classify.hpp index e1a6dd268..3ce43fdcf 100644 --- a/include/universal/number/dd/math/classify.hpp +++ b/include/universal/number/dd/math/classify.hpp @@ -47,9 +47,17 @@ inline bool iszero(const dd& a) { return (std::fpclassify(a.high()) == FP_ZERO); } -bool signbit(dd const& a) { +inline bool signbit(const dd& a) { auto signA = std::copysign(1.0, a.high()); return signA < 0.0; } +/* +inline dd copysign(const dd& a, const dd& b) +{ + auto signA = std::copysign(1.0, a.high()); + auto signB = std::copysign(1.0, b.high()); + return signA != signB ? -a : a; +} +*/ }} // namespace sw::universal diff --git a/static/dd/math/pow.cpp b/static/dd/math/pow.cpp index 5d56ecb2c..c5695bb21 100644 --- a/static/dd/math/pow.cpp +++ b/static/dd/math/pow.cpp @@ -5,32 +5,273 @@ // // This file is part of the universal numbers project, which is released under an MIT Open Source license. #include +#include #include #include -// generate specific test case -template -void GenerateTestCase(Ty fa, Ty fb) { - unsigned precision = 25; - unsigned width = 30; - Ty fref; - sw::universal::dd a, b, ref, v; - a = fa; - b = fb; - fref = std::pow(fa, fb); - ref = fref; - v = sw::universal::pow(a, b); - auto oldPrec = std::cout.precision(); - std::cout << std::setprecision(precision); - std::cout << " -> pow(" << fa << "," << fb << ") = " << std::setw(width) << fref << std::endl; - std::cout << " -> pow( " << a << "," << b << ") = " << v << '\n' << to_binary(v) << '\n'; - std::cout << to_binary(ref) << "\n -> reference\n"; - std::cout << (ref == v ? "PASS" : "FAIL") << std::endl << std::endl; - std::cout << std::setprecision(oldPrec); +// overload the uniform distribution for double-double + +namespace std +{ + + template <> + class uniform_real_distribution< sw::universal::dd > + { + public: + uniform_real_distribution(const sw::universal::dd& low, const sw::universal::dd& high) + : _mDist(low.high(), high.high()) + { + } + + template < typename engine > + sw::universal::dd operator()(engine& eng) { + return sw::universal::dd(_mDist(eng), 0.5 * std::numeric_limits< double >::epsilon() * _mDist(eng)); + } + + private: + uniform_real_distribution< double > _mDist; + }; + } +namespace sw { + namespace universal { + + // generate specific test case + template + void GenerateTestCase(Ty fa, Ty fb) { + unsigned precision = 25; + unsigned width = 30; + Ty fref; + sw::universal::dd a, b, ref, v; + a = fa; + b = fb; + fref = std::pow(fa, fb); + ref = fref; + v = sw::universal::pow(a, b); + auto oldPrec = std::cout.precision(); + std::cout << std::setprecision(precision); + std::cout << " -> pow(" << fa << "," << fb << ") = " << std::setw(width) << fref << std::endl; + std::cout << " -> pow( " << a << "," << b << ") = " << v << '\n' << to_binary(v) << '\n'; + std::cout << to_binary(ref) << "\n -> reference\n"; + std::cout << (ref == v ? "PASS" : "FAIL") << std::endl << std::endl; + std::cout << std::setprecision(oldPrec); + } + + template + int calculateNrOfValidBits(const Real& computed, const Real& expected) { + constexpr double LOG2E = 1.44269504088896340736; + + dd delta = computed - expected; + if (delta == 0.0) { + return dd_max_precision; + } + else { + if (expected == 0.0) { + return static_cast(-std::log(std::fabs(double(computed))) * LOG2E); + } + else { + delta /= expected; + double logOfDelta = std::log(std::fabs(double(delta))) * LOG2E; + return static_cast(-logOfDelta); + } + } + } + + static constexpr int NR_RANDOMS = 500; +#ifdef _DEBUG + static constexpr int PRECISION_THRESHOLD = 85; // in bits: 85 bits is ~ 25.5 digits out of 32 digits +#else + static constexpr int PRECISION_THRESHOLD = 75; // in bits: 85 bits is ~ 25.5 digits out of 32 digits +#endif + + int comparePowWithSqrt(bool reportTestCases, int precisionThreshold = PRECISION_THRESHOLD, int nrOfRandoms = NR_RANDOMS) { + std::default_random_engine generator; + std::uniform_real_distribution< dd > distribution(1.0, 1048576.0); + int maxValidBits, minValidBits; + int nrOfFailedTestCases{ 0 }; + + std::cerr << "smallest number of valid bits of pow(x, 0.5) = "; + if (reportTestCases) std::cerr << '\n'; + maxValidBits = 0; + minValidBits = dd_max_precision; + for (int i = 0; i < nrOfRandoms; ++i) { + dd x = distribution(generator); + dd expected = sqrt(x); + dd computed = pow(x, dd(0.5)); + + int nrOfValidBits = calculateNrOfValidBits(computed, expected); + if (nrOfValidBits < 0) { // something very wrong has happened, provide feedback + ReportValue(computed, "computed"); + ReportValue(expected, "expected"); + } + minValidBits = std::min(minValidBits, nrOfValidBits); + maxValidBits = std::max(maxValidBits, nrOfValidBits); + if (nrOfValidBits < precisionThreshold) { + ++nrOfFailedTestCases; + } + if (reportTestCases) std::cerr << "valid bits pow( " << x << ", 0.5) : " << nrOfValidBits << "\n"; + } + if (minValidBits == dd_max_precision) + std::cerr << "EXACT "; + else + std::cerr << "[ " << minValidBits << ", " << maxValidBits << "] "; + + std::cerr << (nrOfFailedTestCases ? "FAIL\n" : "PASS\n"); + return nrOfFailedTestCases; + } + + int comparePowWithCubeRoot(bool reportTestCases, int precisionThreshold = PRECISION_THRESHOLD, int nrOfRandoms = NR_RANDOMS) { + std::default_random_engine generator; + std::uniform_real_distribution< dd > distribution(1.0, 1048576.0); + int maxValidBits, minValidBits; + int nrOfFailedTestCases{ 0 }; + + std::cerr << "smallest number of valid bits of pow(x, 0.33333...) = "; + if (reportTestCases) std::cerr << '\n'; + maxValidBits = 0; + minValidBits = dd_max_precision; + for (int i = 0; i < nrOfRandoms; ++i) { + dd x = distribution(generator); + dd expected = cbrt(x); + dd computed = pow(x, dd_third); + + int nrOfValidBits = calculateNrOfValidBits(computed, expected); + if (nrOfValidBits < 0) { + ReportValue(computed, "computed"); + ReportValue(expected, "expected"); + } + minValidBits = std::min(minValidBits, nrOfValidBits); + maxValidBits = std::max(maxValidBits, nrOfValidBits); + if (nrOfValidBits < precisionThreshold) { + ++nrOfFailedTestCases; + } + if (reportTestCases) std::cerr << "valid bits pow( " << x << ", 0.3333...) : " << nrOfValidBits << "\n"; + } + if (minValidBits == dd_max_precision) + std::cerr << "EXACT "; + else + std::cerr << "[ " << minValidBits << ", " << maxValidBits << "] "; + + std::cerr << (nrOfFailedTestCases ? "FAIL\n" : "PASS\n"); + return nrOfFailedTestCases; + } + + int comparePowWithSquare(bool reportTestCases, int precisionThreshold = PRECISION_THRESHOLD, int nrOfRandoms = NR_RANDOMS) { + std::default_random_engine generator; + std::uniform_real_distribution< dd > distribution(1.0, 1048576.0); + int maxValidBits, minValidBits; + int nrOfFailedTestCases{ 0 }; + + std::cerr << "smallest number of valid bits of pow(x, 2.0) = "; + if (reportTestCases) std::cerr << '\n'; + maxValidBits = 0; + minValidBits = dd_max_precision; + for (int i = 0; i < nrOfRandoms; ++i) { + dd x = distribution(generator); + dd expected = x * x; + dd computed = pow(x, dd(2.0)); + + int nrOfValidBits = calculateNrOfValidBits(computed, expected); + if (nrOfValidBits < 0) { + ReportValue(computed, "computed"); + ReportValue(expected, "expected"); + } + minValidBits = std::min(minValidBits, nrOfValidBits); + maxValidBits = std::max(maxValidBits, nrOfValidBits); + if (nrOfValidBits < precisionThreshold) { + ++nrOfFailedTestCases; + } + if (reportTestCases) std::cerr << "valid bits pow( " << x << ", 2.0) : " << nrOfValidBits << "\n"; + } + if (minValidBits == dd_max_precision) + std::cerr << "EXACT "; + else + std::cerr << "[ " << minValidBits << ", " << maxValidBits << "] "; + + std::cerr << (nrOfFailedTestCases ? "FAIL\n" : "PASS\n"); + return nrOfFailedTestCases; + } + + int comparePowWithCube(bool reportTestCases, int precisionThreshold = PRECISION_THRESHOLD, int nrOfRandoms = NR_RANDOMS) { + std::default_random_engine generator; + std::uniform_real_distribution< dd > distribution(1.0, 1048576.0); + int maxValidBits, minValidBits; + int nrOfFailedTestCases{ 0 }; + + std::cerr << "smallest number of valid bits of pow(x, 3.0) = "; + if (reportTestCases) std::cerr << '\n'; + maxValidBits = 0; + minValidBits = dd_max_precision; + for (int i = 0; i < nrOfRandoms; ++i) { + dd x = distribution(generator); + dd expected = x * x * x; + dd computed = pow(x, dd(3.0)); + + int nrOfValidBits = calculateNrOfValidBits(computed, expected); + if (nrOfValidBits < 0) { + ReportValue(computed, "computed"); + ReportValue(expected, "expected"); + } + minValidBits = std::min(minValidBits, nrOfValidBits); + maxValidBits = std::max(maxValidBits, nrOfValidBits); + if (nrOfValidBits < precisionThreshold) { + ++nrOfFailedTestCases; + } + if (reportTestCases) std::cerr << "valid bits pow( " << x << ", 3.0) : " << nrOfValidBits << "\n"; + } + if (minValidBits == dd_max_precision) + std::cerr << "EXACT "; + else + std::cerr << "[ " << minValidBits << ", " << maxValidBits << "] "; + + std::cerr << (nrOfFailedTestCases ? "FAIL\n" : "PASS\n"); + return nrOfFailedTestCases; + } + + + int comparePowWithQuadratic(bool reportTestCases, int precisionThreshold = PRECISION_THRESHOLD, int nrOfRandoms = NR_RANDOMS) { + std::default_random_engine generator; + std::uniform_real_distribution< dd > distribution(1.0, 1048576.0); + int maxValidBits, minValidBits; + int nrOfFailedTestCases{ 0 }; + + std::cerr << "smallest number of valid bits of pow(x, 4.0) = "; + if (reportTestCases) std::cerr << '\n'; + maxValidBits = 0; + minValidBits = dd_max_precision; + for (int i = 0; i < nrOfRandoms; ++i) { + dd x = distribution(generator); + dd square = x * x; + dd expected = square * square; + dd computed = pow(x, dd(4.0)); + + int nrOfValidBits = calculateNrOfValidBits(computed, expected); + if (nrOfValidBits < 0) { + ReportValue(computed, "computed"); + ReportValue(expected, "expected"); + } + minValidBits = std::min(minValidBits, nrOfValidBits); + maxValidBits = std::max(maxValidBits, nrOfValidBits); + if (nrOfValidBits < precisionThreshold) { + ++nrOfFailedTestCases; + } + if (reportTestCases) std::cerr << "valid bits pow( " << x << ", 3.0) : " << nrOfValidBits << "\n"; + } + if (minValidBits == dd_max_precision) + std::cerr << "EXACT "; + else + std::cerr << "[ " << minValidBits << ", " << maxValidBits << "] "; + + std::cerr << (nrOfFailedTestCases ? "FAIL\n" : "PASS\n"); + return nrOfFailedTestCases; + } + } +} + + // Regression testing guards: typically set by the cmake configuration, but MANUAL_TESTING is an override -#define MANUAL_TESTING 1 +#define MANUAL_TESTING 0 // REGRESSION_LEVEL_OVERRIDE is set by the cmake file to drive a specific regression intensity // It is the responsibility of the regression test to organize the tests in a quartile progression. //#undef REGRESSION_LEVEL_OVERRIDE @@ -79,6 +320,46 @@ try { return EXIT_SUCCESS; // ignore errors #else + // this is very sad: we are loosing 10 bits of precision when running in Release as compared to Debug: +/* + * Release mode + * double-double mathlib power function validation: results only + * PRECISION_THRESHOLD set to 85 bits, which is approximate 25.7635 digits: out of a total of 32 digits + * smallest number of valid bits of pow(x, 0.5) = [ 81, 92] FAIL + * smallest number of valid bits of pow(x, 0.33333...) = [ 77, 86] FAIL + * smallest number of valid bits of pow(x, 2.0) = [ 78, 91] FAIL + * smallest number of valid bits of pow(x, 3.0) = [ 77, 90] FAIL + * smallest number of valid bits of pow(x, 4.0) = [ 77, 89] FAIL + * double-double mathlib power function validation: FAIL + * + * Debug mode + * double-double mathlib power function validation: results only + * PRECISION_THRESHOLD set to 85 bits, which is approximate 25.7635 digits: out of a total of 32 digits + * smallest number of valid bits of pow(x, 0.5) = [ 92, 110] PASS + * smallest number of valid bits of pow(x, 0.33333...) = [ 91, 110] PASS + * smallest number of valid bits of pow(x, 2.0) = [ 89, 110] PASS + * smallest number of valid bits of pow(x, 3.0) = [ 88, 108] PASS + * smallest number of valid bits of pow(x, 4.0) = [ 88, 108] PASS + * double-double mathlib power function validation: PASS + * + * Setting lower precision threshold to pass regressions + * double-double mathlib power function validation: results only + * PRECISION_THRESHOLD set to 75 bits, which is approximate 22.7325 digits: out of a total of 32 digits + * smallest number of valid bits of pow(x, 0.5) = [ 81, 92] PASS + * smallest number of valid bits of pow(x, 0.33333...) = [ 77, 86] PASS + * smallest number of valid bits of pow(x, 2.0) = [ 78, 91] PASS + * smallest number of valid bits of pow(x, 3.0) = [ 77, 90] PASS + * smallest number of valid bits of pow(x, 4.0) = [ 77, 89] PASS + * double-double mathlib power function validation: PASS + */ + std::cerr << "PRECISION_THRESHOLD set to " << PRECISION_THRESHOLD + << " bits, which is approximate " << (0.3031 * PRECISION_THRESHOLD) << " digits: out of a total of 32 digits\n"; + + nrOfFailedTestCases += comparePowWithSqrt(reportTestCases); + nrOfFailedTestCases += comparePowWithCubeRoot(reportTestCases); + nrOfFailedTestCases += comparePowWithSquare(reportTestCases); + nrOfFailedTestCases += comparePowWithCube(reportTestCases); + nrOfFailedTestCases += comparePowWithQuadratic(reportTestCases); ReportTestSuiteResults(test_suite, nrOfFailedTestCases); return (nrOfFailedTestCases > 0 ? EXIT_FAILURE : EXIT_SUCCESS);