WIP: implementing nextafter for double-double

Author: Ravenwater

include/universal/native/ieee754_numeric.hpp

@@ -16,7 +16,7 @@ namespace sw { namespace universal {
 		typename = typename ::std::enable_if< ::std::is_floating_point<Real>::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

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();

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

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<double>::max();
@@ -89,6 +89,26 @@ namespace sw {
 			std::cout << c << '\n';
 		}
 
+		template<typename Real,
+			typename = typename ::std::enable_if< ::std::is_floating_point<Real>::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<Real>::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<dd>() << std::endl;
 	}
 
-	std::cout << "+---------    doubledouble subnormal behavior   --------+\n";
+	std::cout << "+---------    double-double subnormal behavior   --------+\n";
 	{
 		constexpr double minpos = std::numeric_limits<double>::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<float>::value, float >::type;
+		using double_type = ::std::enable_if< ::std::is_floating_point<double>::value, double >::type;
+//		using bla = ::std::enable_if< ::std::is_floating_point<int>::value, int >::type
+	}
+
 	std::cout << "+---------    numeric_limits of double-double vs IEEE-754   --------+\n";
 	{
 		std::cout << "dd(INFINITY): " << dd(INFINITY) << "\n";

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 <universal/utility/directives.hpp>

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 <universal/utility/directives.hpp>

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 <universal/utility/directives.hpp>

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 <universal/utility/directives.hpp>

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 <universal/utility/directives.hpp>
+#include <cmath>
+#include <concepts>
+#include <cfenv>
+#include <universal/verification/test_suite.hpp>
+
+// 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 = []<std::floating_point Fp>(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<long double>::denorm_min();
+	std::cout << "smallest long double: " << to_binary(denorm_min) << " : " << denorm_min << '\n';
+#endif
+
+
+    auto ulp_progression = []<std::floating_point Fp>(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;
+}

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 <universal/utility/directives.hpp>