diff --git a/include/universal/native/ieee754_numeric.hpp b/include/universal/native/ieee754_numeric.hpp index 7e4d81ab6..b1b328f74 100644 --- a/include/universal/native/ieee754_numeric.hpp +++ b/include/universal/native/ieee754_numeric.hpp @@ -16,7 +16,7 @@ namespace sw { namespace universal { typename = typename ::std::enable_if< ::std::is_floating_point::value, Real >::type > inline Real ulp(const Real& a) { - return std::nextafter(a, a + a/2.0f) - a; + return std::nextafter(a, Real(INFINITY)) - a; } // check if the floating-point number is zero diff --git a/include/universal/number/dd/dd_impl.hpp b/include/universal/number/dd/dd_impl.hpp index 2f2c76c70..7df6af895 100644 --- a/include/universal/number/dd/dd_impl.hpp +++ b/include/universal/number/dd/dd_impl.hpp @@ -29,7 +29,8 @@ namespace sw { namespace universal { -// fwd references to free functions used in to_digits() +// fwd references to free functions +dd operator-(const dd& lhs, const dd&); dd operator*(const dd& lhs, const dd&); std::ostream& operator<<(std::ostream&, const dd&); dd pown(const dd&, int); @@ -865,6 +866,14 @@ inline std::string to_binary(const dd& number, bool bNibbleMarker = false) { //////////////////////// math functions ///////////////////////////////// +inline dd ulp(const dd& a) { + double hi{ a.high() }; + double nlo = std::nextafter(a.low(), INFINITY); + dd n(hi, nlo); + + return n - a; +} + inline dd abs(dd a) { double hi = a.high(); double lo = a.low(); diff --git a/include/universal/number/dd/math/next.hpp b/include/universal/number/dd/math/next.hpp index ca896806e..399194719 100644 --- a/include/universal/number/dd/math/next.hpp +++ b/include/universal/number/dd/math/next.hpp @@ -24,25 +24,16 @@ Return Value - And math_errhandling has MATH_ERRNO set: the global variable errno is set to ERANGE. - And math_errhandling has MATH_ERREXCEPT set: FE_OVERFLOW is raised. */ -dd nextafter(dd x, dd target) { - if (x == target) return target; - if (target.isnan()) { - if (x.isneg()) { - --x; - } - else { - ++x; - } +dd nextafter(const dd& x, const dd& target) { + double hi{ x.high() }; + double lo; + if (x < target) { + lo = std::nextafter(x.low(), +INFINITY); } else { - if (x > target) { - --x; - } - else { - ++x; - } + lo = std::nextafter(x.low(), -INFINITY); } - return x; + return dd(hi, lo); } /* TODO diff --git a/static/dd/api/api.cpp b/static/dd/api/api.cpp index 4d20afe61..3ba7ff6c2 100644 --- a/static/dd/api/api.cpp +++ b/static/dd/api/api.cpp @@ -61,7 +61,7 @@ namespace sw { ostr << str << '\n'; } - void construct_largest_doubledouble() { + void construct_largest_double_double() { using Scalar = dd; double firstLimb = std::numeric_limits::max(); @@ -89,6 +89,26 @@ namespace sw { std::cout << c << '\n'; } + template::value, Real >::type + > + Real emulateNextAfter(Real x, Real y) { + if (x == y) return y; + int direction = (x < y) ? 1 : -1; + Real eps = std::numeric_limits::epsilon(); + return x + direction * eps; + } + + void ulp_progression(const std::string& tag, const dd& start) { + std::cout << tag; + for (dd from = start, to, delta; + (delta = (to = nextafter(from, +INFINITY)) - from) < 10.0; + from *= 10.0) { + std::cout << "ulp(" << std::scientific << std::setprecision(0) << from + << ") gives " << to_binary(ulp(from)) << " : " + << std::fixed << std::setprecision(6) << ulp(from) << '\n'; + } + } } } @@ -183,22 +203,22 @@ try { } // report on the dynamic range of some standard configurations - std::cout << "+--------- Dynamic range doubledouble configurations ---------+\n"; + std::cout << "+--------- Dynamic range double-double configurations ---------+\n"; { dd a; // uninitialized a.maxpos(); - std::cout << "maxpos doubledouble : " << to_binary(a) << " : " << a << '\n'; + std::cout << "maxpos double-double : " << to_binary(a) << " : " << a << '\n'; a.setbits(0x0080); // positive min normal - std::cout << "minnorm doubledouble : " << to_binary(a) << " : " << a << '\n'; + std::cout << "minnorm double-double : " << to_binary(a) << " : " << a << '\n'; a.minpos(); - std::cout << "minpos doubledouble : " << to_binary(a) << " : " << a << '\n'; + std::cout << "minpos double-double : " << to_binary(a) << " : " << a << '\n'; a.zero(); std::cout << "zero : " << to_binary(a) << " : " << a << '\n'; a.minneg(); - std::cout << "minneg doubledouble : " << to_binary(a) << " : " << a << '\n'; + std::cout << "minneg double-double : " << to_binary(a) << " : " << a << '\n'; a.maxneg(); - std::cout << "maxneg doubledouble : " << to_binary(a) << " : " << a << '\n'; + std::cout << "maxneg double-double : " << to_binary(a) << " : " << a << '\n'; std::cout << "---\n"; } @@ -302,7 +322,7 @@ try { std::cout << dynamic_range
() << std::endl; } - std::cout << "+--------- doubledouble subnormal behavior --------+\n"; + std::cout << "+--------- double-double subnormal behavior --------+\n"; { constexpr double minpos = std::numeric_limits::min(); std::cout << to_binary(minpos) << " : " << minpos << '\n'; @@ -316,7 +336,7 @@ try { } } - std::cout << "+--------- special value properties doubledouble vs IEEE-754 --------+\n"; + std::cout << "+--------- special value properties double-double vs IEEE-754 --------+\n"; { float fa; fa = NAN; @@ -331,14 +351,25 @@ try { dd a(fa); if ((a < 0.0f && a > 0.0f && a != 0.0f)) { - std::cout << "doubledouble (dd) is incorrectly implemented\n"; + std::cout << "double-double (dd) is incorrectly implemented\n"; ++nrOfFailedTestCases; } else { - std::cout << "dd NAN has no sign\n"; + std::cout << "double-double (dd) NAN has no sign\n"; } } + std::cout << "---------- Unit in the Last Place --------+\n"; + { + dd a{ 1.0 }; + + ulp_progression("\nULP progression for dd:\n", dd(10.0)); + + using float_type = ::std::enable_if< ::std::is_floating_point::value, float >::type; + using double_type = ::std::enable_if< ::std::is_floating_point::value, double >::type; +// using bla = ::std::enable_if< ::std::is_floating_point::value, int >::type + } + std::cout << "+--------- numeric_limits of double-double vs IEEE-754 --------+\n"; { std::cout << "dd(INFINITY): " << dd(INFINITY) << "\n"; diff --git a/static/native/float/fractionviz.cpp b/static/native/float/fractionviz.cpp index df6deed09..413fda150 100644 --- a/static/native/float/fractionviz.cpp +++ b/static/native/float/fractionviz.cpp @@ -1,6 +1,7 @@ // fractionviz.cpp : fraction bits visualization of native IEEE-754 types // -// Copyright (C) 2017-2023 Stillwater Supercomputing, Inc. +// Copyright (C) 2017 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT // // This file is part of the universal numbers project, which is released under an MIT Open Source license. #include diff --git a/static/native/float/ieee754.cpp b/static/native/float/ieee754.cpp index 111407bd5..81bdafb78 100644 --- a/static/native/float/ieee754.cpp +++ b/static/native/float/ieee754.cpp @@ -1,6 +1,7 @@ // ieee754.cpp : native IEEE-754 operations // -// Copyright (C) 2017-2023 Stillwater Supercomputing, Inc. +// Copyright (C) 2017 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT // // This file is part of the universal numbers project, which is released under an MIT Open Source license. #include diff --git a/static/native/float/manipulators.cpp b/static/native/float/manipulators.cpp index 023c725d0..a76cd70cc 100644 --- a/static/native/float/manipulators.cpp +++ b/static/native/float/manipulators.cpp @@ -1,6 +1,7 @@ // manipulators.cpp : test suite runner for native floating-point manipulators // -// Copyright (C) 2017-2023 Stillwater Supercomputing, Inc. +// Copyright (C) 2017 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT // // This file is part of the universal numbers project, which is released under an MIT Open Source license. #include diff --git a/static/native/float/mathlib.cpp b/static/native/float/mathlib.cpp index 04b5be976..2d46ef8aa 100644 --- a/static/native/float/mathlib.cpp +++ b/static/native/float/mathlib.cpp @@ -1,6 +1,7 @@ // mathlib.cpp : universal math library wrapper // -// Copyright (C) 2017-2023 Stillwater Supercomputing, Inc. +// Copyright (C) 2017 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT // // This file is part of the universal numbers project, which is released under an MIT Open Source license. #include diff --git a/static/native/float/nextafter.cpp b/static/native/float/nextafter.cpp new file mode 100644 index 000000000..10b452163 --- /dev/null +++ b/static/native/float/nextafter.cpp @@ -0,0 +1,146 @@ +// nextafter.cpp : exploration of the nextafter stdlib function to manipulate units in the last place +// +// Copyright (C) 2017 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT +// +// This file is part of the universal numbers project, which is released under an MIT Open Source license. +#include +#include +#include +#include +#include + +// Regression testing guards: typically set by the cmake configuration, but MANUAL_TESTING is an override +#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 +#ifndef REGRESSION_LEVEL_OVERRIDE +#undef REGRESSION_LEVEL_1 +#undef REGRESSION_LEVEL_2 +#undef REGRESSION_LEVEL_3 +#undef REGRESSION_LEVEL_4 +#define REGRESSION_LEVEL_1 1 +#define REGRESSION_LEVEL_2 1 +#define REGRESSION_LEVEL_3 1 +#define REGRESSION_LEVEL_4 1 +#endif + +int main() +try { + using namespace sw::universal; + + std::string test_suite = "nextafter test"; + std::string test_tag = "nextafter"; + bool reportTestCases = false; + int nrOfFailedTestCases = 0; + + ReportTestSuiteHeader(test_suite, reportTestCases); + + { + float from = 0, to = std::nextafter(from, 1.f); + std::cout << "The next representable float after " << std::setprecision(20) << from + << " is " << to + << std::hexfloat << " (" << to << ")\n" << std::defaultfloat; + } + + { + float from = 1, to = std::nextafter(from, 2.f); + std::cout << "The next representable float after " << from << " is " << to + << std::hexfloat << " (" << to << ")\n" << std::defaultfloat; + } + + { + double from = std::nextafter(0.1, 0), to = 0.1; + std::cout << "The number 0.1 lies between two valid doubles:\n" + << std::setprecision(56) << " " << from + << std::hexfloat << " (" << from << ')' << std::defaultfloat + << "\nand " << to << std::hexfloat << " (" << to << ")\n" + << std::defaultfloat << std::setprecision(20); + } + + { + std::cout << "\nDifference between nextafter and nexttoward:\n"; + float from = 0.0f; + long double dir = std::nextafter(from, 1.0L); // first subnormal long double + float x = std::nextafter(from, dir); // first converts dir to float, giving 0 + std::cout << "With nextafter, next float after " << from << " is " << x << '\n'; + x = std::nexttoward(from, dir); + std::cout << "With nexttoward, next float after " << from << " is " << x << '\n'; + } + + std::cout << "\nSpecial values:\n"; + { + // #pragma STDC FENV_ACCESS ON + std::feclearexcept(FE_ALL_EXCEPT); + double from4 = DBL_MAX, to4 = std::nextafter(from4, INFINITY); + std::cout << "The next representable double after " << std::setprecision(6) + << from4 << std::hexfloat << " (" << from4 << ')' + << std::defaultfloat << " is " << to4 + << std::hexfloat << " (" << to4 << ")\n" << std::defaultfloat; + + if (std::fetestexcept(FE_OVERFLOW)) + std::cout << " raised FE_OVERFLOW\n"; + if (std::fetestexcept(FE_INEXACT)) + std::cout << " raised FE_INEXACT\n"; + } // end FENV_ACCESS block + + { + float from = 0.0, to = std::nextafter(from, -0.0); + std::cout << "std::nextafter(+0.0, -0.0) gives " << std::fixed << to << '\n'; + } + + auto precision_loss_demo = [](const auto rem, const Fp start) + { + std::cout << rem; + for (Fp from = start, to, delta; + (delta = (to = std::nextafter(from, +INFINITY)) - from) < Fp(10.0); + from *= Fp(10.0)) { + std::cout << "nextafter(" << std::scientific << std::setprecision(0) << from + << ", INF) gives " << std::fixed << std::setprecision(6) << to + << "; delta = " << delta << '\n'; + } + }; + + precision_loss_demo("\nPrecision loss demo for float:\n", 10.0f); + precision_loss_demo("\nPrecision loss demo for double:\n", 10.0e9); + +#if LONG_DOUBLE_SUPPORT + precision_loss_demo("\nPrecision loss demo for long double:\n", 10.0e17L); + + constexpr long double denorm_min = std::numeric_limits::denorm_min(); + std::cout << "smallest long double: " << to_binary(denorm_min) << " : " << denorm_min << '\n'; +#endif + + + auto ulp_progression = [](const auto rem, const Fp start) + { + std::cout << rem; + for (Fp from = start, to, delta; + (delta = (to = std::nextafter(from, +INFINITY)) - from) < Fp(10.0); + from *= Fp(10.0)) { + std::cout << "ulp(" << std::scientific << std::setprecision(0) << from + << ") gives " << to_binary(ulp(from)) << " : " + << std::fixed << std::setprecision(6) << ulp(from) << '\n'; + } + }; + + ulp_progression("\nULP progression for float:\n", 10.0f); + ulp_progression("\nULP progression for double:\n", 10.0e9); + + std::cout << '\n'; + ReportTestSuiteResults(test_suite, nrOfFailedTestCases); + return (nrOfFailedTestCases > 0 ? EXIT_FAILURE : EXIT_SUCCESS); +} +catch (char const* msg) { + std::cerr << msg << '\n'; + return EXIT_FAILURE; +} +catch (const std::runtime_error& err) { + std::cerr << "Uncaught runtime exception: " << err.what() << std::endl; + return EXIT_FAILURE; +} +catch (...) { + std::cerr << "Caught unknown exception" << '\n'; + return EXIT_FAILURE; +} diff --git a/static/native/float/representable.cpp b/static/native/float/representable.cpp index 679a61734..888594283 100644 --- a/static/native/float/representable.cpp +++ b/static/native/float/representable.cpp @@ -1,6 +1,7 @@ // representable.cpp : check if a ratio is representable without error as a floating-point value // -// Copyright (C) 2017-2023 Stillwater Supercomputing, Inc. +// Copyright (C) 2017 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT // // This file is part of the universal numbers project, which is released under an MIT Open Source license. #include diff --git a/static/native/integer/bit_manipulation.cpp b/static/native/integer/bit_manipulation.cpp index 328e9f501..bed6a8c8e 100644 --- a/static/native/integer/bit_manipulation.cpp +++ b/static/native/integer/bit_manipulation.cpp @@ -1,6 +1,7 @@ // bit_manipulation.cpp : test runner for bit manipulation of native integers // -// Copyright (C) 2017-2022 Stillwater Supercomputing, Inc. +// Copyright (C) 2017 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT // // This file is part of the universal numbers project, which is released under an MIT Open Source license. #include @@ -44,7 +45,7 @@ void TestNLZ() { uint8_t a = 0x1; for (uint32_t i = 0; i < 8; ++i) { int shift = nlz(a); - std::cout << " shift = " << shift << " : " << to_binary(a, 8, true) << '\n'; + std::cout << " shift = " << shift << " : " << to_binary(a, true, 8) << '\n'; a <<= 1; } } @@ -53,7 +54,7 @@ void TestNLZ() { uint16_t a = 0x1; for (uint32_t i = 0; i < 16; ++i) { int shift = nlz(a); - std::cout << " shift = " << shift << " : " << to_binary(a, 16, true) << '\n'; + std::cout << " shift = " << shift << " : " << to_binary(a, true, 16) << '\n'; a <<= 1; } } @@ -61,7 +62,7 @@ void TestNLZ() { uint32_t a = 0x1; for (uint32_t i = 0; i < 32; ++i) { int shift = nlz(a); - std::cout << " shift = " << shift << " : " << to_binary(a, 32, true) << '\n'; + std::cout << " shift = " << shift << " : " << to_binary(a, true, 32) << '\n'; a <<= 1; } } @@ -69,7 +70,7 @@ void TestNLZ() { uint64_t a = 0x1; for (uint32_t i = 0; i < 64; ++i) { int shift = nlz(a); - std::cout << " shift = " << shift << " : " << to_binary(a, 64, true) << '\n'; + std::cout << " shift = " << shift << " : " << to_binary(a, true, 64) << '\n'; a <<= 1; } }