WIP: trying out cbrt as a test math function for specialization with …

…generics fallback

include/universal/math/functions.hpp

@@ -9,10 +9,11 @@
 // properties of a number
 #include <universal/math/functions/isrepresentable.hpp>
 
+// generic mathematical functions
+#include <universal/math/functions/cbrt.hpp>
+
 // special functions
+#include <universal/math/functions/horners.hpp>
 #include <universal/math/functions/factorial.hpp>
 #include <universal/math/functions/binomial.hpp>
 #include <universal/math/functions/loss.hpp>
-
-// generic mathematical functions
-#include <universal/math/functions/horners.hpp>

include/universal/math/functions/cbrt.hpp

@@ -0,0 +1,18 @@
+#pragma once
+// cbrt.hpp: generic implementation of a cubic root of a Real
+//
+// 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.
+
+namespace sw::universal {
+
+	template<typename Real>
+	Real cbrt(const Real& x) {
+		assert(x >= Real(0));
+		return std::cbrt(x);
+	}
+
+}}  // namespace sw::function
+

include/universal/number/dd/dd_impl.hpp

@@ -1251,38 +1251,6 @@ inline dd reciprocal(const dd& a) {
 /////////////////////////////////////////////////////////////////////////////
 //	power functions
 
-/// <summary>
-/// cbrt is cube root function, that is, x^1/3
-/// </summary>
-/// <param name="a">input</param>
-/// <returns>cube root of a</returns>
-inline dd cbrt(const dd& a) {
-	using std::pow;
-	if (!a.isfinite() || 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;
 

include/universal/number/dd/math/functions/cbrt.hpp

@@ -0,0 +1,43 @@
+#pragma once
+// cbrt.hpp: cbrt function for double-double floating-point
+//
+// 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.
+
+namespace sw { namespace universal {
+
+	/// <summary>
+	/// cbrt is cube root function, that is, x^1/3
+	/// </summary>
+	/// <param name="a">input</param>
+	/// <returns>cube root of a</returns>
+	inline dd cbrt(const dd& a) {
+		using std::pow;
+		if (!a.isfinite() || 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);
+	}
+
+}} // namespace sw::universal

include/universal/number/dd/mathlib.hpp

@@ -21,3 +21,5 @@
 #include <universal/number/dd/math/functions/sqrt.hpp>
 #include <universal/number/dd/math/functions/trigonometry.hpp>
 #include <universal/number/dd/math/functions/truncate.hpp>
+
+#include <universal/number/dd/math/functions/cbrt.hpp>

include/universal/number/qd/qd_impl.hpp

@@ -1131,19 +1131,7 @@ class qd {
 
 ////////////////////////  precomputed constants of note  /////////////////////////////////
 
-// precomputed quad-double constants courtesy of constants example program
-
-
-
-
-
-
-
-
-
-
-
-
+// precomputed quad-double constants 
 
 constexpr qd qd_max(1.79769313486231570815e+308, 9.97920154767359795037e+291);
 

static/dd/math/sqrt.cpp

@@ -8,29 +8,51 @@
 #include <universal/number/dd/dd.hpp>
 #include <universal/verification/test_suite.hpp>
 
-// generate specific test case 
-template<typename Ty>
-void GenerateSqrtTestCase(Ty fa) {
-	unsigned precision = 25;
-	//unsigned width = 30;
-	Ty fref;
-	sw::universal::dd a, ref, v;
-	a = fa;
-	fref = std::sqrt(fa);
-	ref = fref;
-	v = sw::universal::sqrt(a);
-	auto oldPrec = std::cout.precision();
-	std::cout << std::setprecision(precision);
-	std::cout << " -> sqrt(" << fa << ") = " << fref << std::endl;
-	std::cout << " -> sqrt( " << a << ") = " << 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);
-}
 
 namespace sw {
 	namespace universal {
 
+		// generate specific test case 
+		template<typename Ty>
+		void GenerateSqrtTestCase(Ty fa) {
+			unsigned precision = 25;
+			//unsigned width = 30;
+			Ty fref;
+			sw::universal::dd a, ref, v;
+			a = fa;
+			fref = std::sqrt(fa);
+			ref = fref;
+			v = sw::universal::sqrt(a);
+			auto defaultPrecision = std::cout.precision();
+			std::cout << std::setprecision(precision);
+			std::cout << " -> sqrt(" << fa << ") = " << fref << std::endl;
+			std::cout << " -> sqrt( " << a << ") = " << 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(defaultPrecision);
+		}
+
+		// generate specific test case 
+		template<typename Ty>
+		void GenerateCbrtTestCase(Ty fa) {
+			unsigned precision = 25;
+			//unsigned width = 30;
+			Ty fref;
+			sw::universal::dd a, ref, v;
+			a = fa;
+			fref = std::cbrt(fa);
+			ref = fref;
+			v = sw::universal::cbrt(a);
+			auto defaultPrecision = std::cout.precision();
+			std::cout << std::setprecision(precision);
+			std::cout << " -> sqrt(" << fa << ") = " << fref << std::endl;
+			std::cout << " -> sqrt( " << a << ") = " << 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(defaultPrecision);
+		}
+
+
 		template<typename DoubleDouble>
 		int VerifySqrtFunction(bool reportTestCases, DoubleDouble a) {
 			int nrOfFailedTestCases{ 0 };
@@ -46,6 +68,24 @@ namespace sw {
 			}
 			return nrOfFailedTestCases;
 		}
+
+		template<typename DoubleDouble>
+		int VerifyCbrtFunction(bool reportTestCases, DoubleDouble a) {
+			using std::cbrt;
+			int nrOfFailedTestCases{ 0 };
+			DoubleDouble b{ a };
+			for (int i = 0; i < 6; ++i) {
+				a *= a * a;
+				dd c = cbrt(a);
+				if (b != c) {
+					if (reportTestCases) std::cerr << "FAIL : " << b << " != " << c << '\n';
+					++nrOfFailedTestCases;
+				}
+				b *= b * b;
+			}
+			return nrOfFailedTestCases;
+		}
+
 	}
 }
 
@@ -70,8 +110,8 @@ int main()
 try {
 	using namespace sw::universal;
 
-	std::string test_suite  = "double-double mathlib sqrt function validation";
-	std::string test_tag    = "sqrt";
+	std::string test_suite  = "double-double mathlib sqrt/cbrt function validation";
+	std::string test_tag    = "sqrt/cbrt";
 	bool reportTestCases    = true;
 	int nrOfFailedTestCases = 0;
 
@@ -94,9 +134,15 @@ try {
 	return EXIT_SUCCESS;   // ignore errors
 #else
 
+	test_tag = "sqrt";
 	nrOfFailedTestCases += ReportTestResult(VerifySqrtFunction(reportTestCases, dd(2.0)), "sqrt(dd > 1.0)", test_tag);
 	nrOfFailedTestCases += ReportTestResult(VerifySqrtFunction(reportTestCases, dd(0.5)), "sqrt(dd < 1.0)", test_tag);
 
+	test_tag = "cbrt";
+	nrOfFailedTestCases += ReportTestResult(VerifyCbrtFunction(reportTestCases, dd(2.0)), "cbrt(dd > 1.0)", test_tag);
+	nrOfFailedTestCases += ReportTestResult(VerifyCbrtFunction(reportTestCases, dd(0.5)), "cbrt(dd < 1.0)", test_tag);
+
+	nrOfFailedTestCases += ReportTestResult(VerifyCbrtFunction(reportTestCases, 2.0), "cbrt(double > 1.0)", test_tag);
 
 	ReportTestSuiteResults(test_suite, nrOfFailedTestCases);
 	return (nrOfFailedTestCases > 0 ? EXIT_FAILURE : EXIT_SUCCESS);

static/qd/api/constants.cpp

@@ -351,37 +351,38 @@ qd_1_sqrt2      : 7.07106781186547524400844362104849e-01 vs 7.071067811865475244
 		std::cout << std::left << std::setw(15) << e.name << " : " << c << " vs " << e.value << " : " << to_quad(c) << " : " << error << '\n';
 	}
 
+
+	/* 
 	{
 
 		qd a{ 2.0 };
 		qd sqrt2 = sqrt(a);
-		std::cout << "sqrt(2.0) " << to_quad(sqrt2) << '\n';
-
+		std::cout << "sqrt(2.0)   " << to_quad(sqrt2) << '\n';
 		a = 3.0;
-		std::cout << "sqrt(3.0) " << to_quad(sqrt(a)) << '\n';
+		std::cout << "sqrt(3.0)   " << to_quad(sqrt(a)) << '\n';
 		a = 5.0;
-		std::cout << "sqrt(5.0) " << to_quad(sqrt(a)) << '\n';
+		std::cout << "sqrt(5.0)   " << to_quad(sqrt(a)) << '\n';
 
 		std::cout << "1/sqrt(2.0) " << to_quad(reciprocal(sqrt2)) << '\n';
 
 		a = sqrt(qd_pi);
-		std::cout << "2 / sqrtpi " << to_quad(2.0 / a) << '\n';
+		std::cout << "2/sqrt(pi)  " << to_quad(2.0 / a) << '\n';
 	}
-	/* 
+
 		AMD x86	
 		Debug build
-		sqrt(2.0) ( 1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683509e-50)
-		sqrt(3.0) ( 1.7320508075688772, 1.0035084221806903e-16, -1.4959542475733896e-33, 5.3061475632961675e-50)
-		sqrt(5.0) ( 2.2360679774997898, -1.0864230407365012e-16, 5.3086504167631309e-33, -6.6099839950042175e-50)
-		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341755e-50)
-		2 / sqrtpi ( 1.1283791670955126, 1.5335459613165881e-17, -4.7656845966936863e-34, -2.0077946616552625e-50)
+		sqrt(2.0)   ( 1.4142135623730951 , -9.6672933134529135e-17,  4.1386753086994136e-33,  4.9355469914683509e-50)
+		sqrt(3.0)   ( 1.7320508075688772 ,  1.0035084221806903e-16, -1.4959542475733896e-33,  5.3061475632961675e-50)
+		sqrt(5.0)   ( 2.2360679774997898 , -1.0864230407365012e-16,  5.3086504167631309e-33, -6.6099839950042175e-50)
+		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466567264567e-17,  2.0693376543497068e-33,  2.4677734957341755e-50)
+		2/sqrt(pi)  ( 1.1283791670955126 ,  1.5335459613165881e-17, -4.7656845966936863e-34, -2.0077946616552625e-50)
 
 		Release build
-		sqrt(2.0) ( 1.4142135623730951, -9.6672933134529135e-17, 4.1386753203466335e-33, -3.3032885712977947e-49)
-		sqrt(3.0) ( 1.7320508075688772, 1.0035084221806903e-16, -1.4959542883445281e-33, 5.0676801879243325e-50)
-		sqrt(5.0) ( 2.2360679774997898, -1.0864230407365012e-16, 5.3086504310320564e-33, -2.7103246582355688e-49)
-		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466312625432e-17, -3.039266735626984e-34, -1.350504809842679e-50)
-		2 / sqrtpi ( 1.1283791670955126, 1.5335458971746789e-17, 2.6579683555126638e-34, -1.683757146154259e-50)
+		sqrt(2.0)   ( 1.4142135623730951 , -9.6672933134529135e-17,  4.1386753203466335e-33, -3.3032885712977947e-49)
+		sqrt(3.0)   ( 1.7320508075688772 ,  1.0035084221806903e-16, -1.4959542883445281e-33,  5.0676801879243325e-50)
+		sqrt(5.0)   ( 2.2360679774997898 , -1.0864230407365012e-16,  5.3086504310320564e-33, -2.7103246582355688e-49)
+		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466312625432e-17, -3.039266735626984e-34 , -1.350504809842679e-50)
+		2/sqrt(pi)  ( 1.1283791670955126 ,  1.5335458971746789e-17,  2.6579683555126638e-34, -1.683757146154259e-50)
 
 		difference       sqrt2    : -1.16472195516512003859185071422508e-41
 		difference       sqrt3    : +4.07711385546630871610406778813869e-41
@@ -391,19 +392,19 @@ qd_1_sqrt2      : 7.07106781186547524400844362104849e-01 vs 7.071067811865475244
 
 		Apple M2 
 		Debug build
-		sqrt(2.0) ( 1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683509e-50)
-		sqrt(3.0) ( 1.7320508075688772, 1.0035084221806903e-16, -1.4959542475733896e-33, 5.3061475632961675e-50)
-		sqrt(5.0) ( 2.2360679774997898, -1.0864230407365012e-16, 5.3086504167631309e-33, -6.6099839950042175e-50)
-		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341755e-50)
-		2 / sqrtpi ( 1.1283791670955126, 1.5335459613165881e-17, -4.7656845966936863e-34, -2.0077946616552625e-50)
+		sqrt(2.0)   ( 1.4142135623730951 , -9.6672933134529135e-17,  4.1386753086994136e-33,  4.9355469914683509e-50)
+		sqrt(3.0)   ( 1.7320508075688772 ,  1.0035084221806903e-16, -1.4959542475733896e-33,  5.3061475632961675e-50)
+		sqrt(5.0)   ( 2.2360679774997898 , -1.0864230407365012e-16,  5.3086504167631309e-33, -6.6099839950042175e-50)
+		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466567264567e-17,  2.0693376543497068e-33,  2.4677734957341755e-50)
+		2/sqrt(pi)  ( 1.1283791670955126 ,  1.5335459613165881e-17, -4.7656845966936863e-34, -2.0077946616552625e-50)
 
 		Release build
-		sqrt(2.0) ( 1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683509e-50)
-		sqrt(3.0) ( 1.7320508075688772, 1.0035084221806903e-16, -1.4959542475733896e-33, 5.3061475632961675e-50)
-		sqrt(5.0) ( 2.2360679774997898, -1.0864230407365012e-16, 5.3086504167631309e-33, -6.6099839950042175e-50)
-		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466567264567e-17, 2.0693376543497068e-33, 2.4677734957341755e-50)
-		2 / sqrtpi ( 1.1283791670955126, 1.5335459613165881e-17, -4.7656845966936863e-34, -2.0077946616552625e-50)
-	*/
+		sqrt(2.0)   ( 1.4142135623730951 , -9.6672933134529135e-17,  4.1386753086994136e-33,  4.9355469914683509e-50)
+		sqrt(3.0)   ( 1.7320508075688772 ,  1.0035084221806903e-16, -1.4959542475733896e-33,  5.3061475632961675e-50)
+		sqrt(5.0)   ( 2.2360679774997898 , -1.0864230407365012e-16,  5.3086504167631309e-33, -6.6099839950042175e-50)
+		1/sqrt(2.0) ( 0.70710678118654757, -4.8336466567264567e-17,  2.0693376543497068e-33,  2.4677734957341755e-50)
+		2/sqrt(pi)  ( 1.1283791670955126 ,  1.5335459613165881e-17, -4.7656845966936863e-34, -2.0077946616552625e-50)
+
 
 	{
 		qd debug_sqrt2(1.4142135623730951, -9.6672933134529135e-17, 4.1386753086994136e-33, 4.9355469914683509e-50);
@@ -424,6 +425,7 @@ qd_1_sqrt2      : 7.07106781186547524400844362104849e-01 vs 7.071067811865475244
 		std::cout << "difference       1_sqrt2  : " << (debug_1_sqrt2 - release_1_sqrt2) << '\n';
 		std::cout << "difference       2_sqrtpi : " << (debug_2_sqrtpi - release_2_sqrtpi) << '\n';
 	}
+	*/
 
 	std::cout << std::setprecision(oldPrec);