Adding Kahan's Smooth Surprise to Applications. Also moved Rump to th…

…e accuracy/mathemetics area

applications/accuracy/mathematics/kahan_smooth_surprise.cpp

@@ -0,0 +1,77 @@
+// harmonic_series.cpp: experiments with mixed-precision representations of the Harmonic Series
+//
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+// SPDX-License-Identifier: MIT
+//
+// This file is part of the UNIVERSAL project, which is released under an MIT Open Source license.
+#include <universal/utility/directives.hpp>
+#include <iostream>
+#include <iomanip>
+
+#include <universal/number/posit/posit.hpp>
+#include <universal/number/cfloat/cfloat.hpp>
+#include <universal/number/dd/dd.hpp>
+#include <universal/number/qd/qd.hpp>
+
+#include <universal/number/rational/rational.hpp>
+
+namespace sw {
+	namespace universal {
+
+
+
+		template<typename Scalar>
+		Scalar f(Scalar x) {
+			using std::log, std::abs;
+			return log(abs(3 * (1 - x) + 1)) / 80 + x * x + 1;
+		}
+
+		template<typename Scalar>
+		void report_on_f(Scalar x) {
+			Scalar f_of_x = f(x);
+			std::cout << "f(" << x << ") = " << to_binary(f_of_x) << " : " << f_of_x << '\n';
+		}
+
+		// Kahan's Smooth Surprise: minimize log(abs(3(1-x) + 1))/80 + x^2 + 1 in the interval [0.8, 2.0]
+		template<typename Scalar>
+		Scalar smooth_surprise(int samples = 10) {
+
+			Scalar x{ 0.8 };
+			Scalar dx = Scalar(1.2) / Scalar(samples);
+			Scalar min = std::numeric_limits<Scalar>::max();
+			for (int i = 0; i < samples; ++i) {
+				Scalar y = f(x);
+				if (y < min) min = y;
+				//std::cout << std::setw(10) << i << " : " << x << " : " << y << '\n';
+				x += dx;
+			}
+			return min;
+
+		}
+
+
+	}
+}
+
+int main()
+try {
+	using namespace sw::universal;
+
+	int samples = 1024 * 512;
+	std::cout << "minimum = " << smooth_surprise<float>(samples) << '\n';
+
+	float f_4{ 4.0 }, f_3{ 3.0 }, f_4_3 = f_4 / f_3;
+	report_on_f(f_4_3);
+	double d_4{ 4.0 }, d_3{ 3.0 }, d_4_3 = d_4 / d_3;
+	report_on_f(d_4_3);
+	dd dd_4{ 4.0 }, dd_3{ 3.0 }, dd_4_3 = dd_4 / dd_3;
+	report_on_f(dd_4_3);
+	qd qd_4{ 4.0 }, qd_3{ 3.0 }, qd_4_3 = qd_4 / qd_3;
+	report_on_f(qd_4_3);
+
+	rational<64> r_4{ 4 }, r_3{ 3 }, r_4_3 = r_4 / r_3;
+	report_on_f(r_4_3);
+}
+catch(...) {
+}
+

applications/precision/numeric/rump_equation.cpp → applications/accuracy/mathematics/rump_equation.cpp

@@ -16,6 +16,8 @@
 #include <universal/number/cfloat/cfloat.hpp>
 #include <universal/number/posit/posit.hpp>
 #include <universal/number/lns/lns.hpp>
+#include <universal/number/dd/dd.hpp>
+#include <universal/number/qd/qd.hpp>
 
 #include <universal/blas/blas.hpp>
 
@@ -195,12 +197,12 @@ There are a couple fpbench versions of it:  https://fpbench.org/benchmarks.html#
 		"posit80",
 		"posit128",
 		"posit156",
-		"bfloat16",
-		"bfloat32",
-		"bfloat64",
-		"bfloat80",
-		"bfloat100",
-		"bfloat120"
+		"cfloat16",
+		"cfloat32",
+		"cfloat64",
+		"cfloat80",
+		"dd",
+		"qd"
 	};
 	Labels column = {
 		"Rump1",
@@ -221,6 +223,12 @@ There are a couple fpbench versions of it:  https://fpbench.org/benchmarks.html#
 	GenerateRow<posit<80, 2>>(a, b, table, 8);
 	GenerateRow<posit<128, 2>>(a, b, table, 9);
 	GenerateRow<posit<156, 2>>(a, b, table, 10);
+	GenerateRow<cfloat<16, 11, std::uint16_t, true>>(a, b, table, 11);
+	GenerateRow<cfloat<32, 11, std::uint32_t, true>>(a, b, table, 12);
+	GenerateRow<cfloat<64, 11, std::uint64_t, true>>(a, b, table, 13);
+	GenerateRow<cfloat<80, 11, std::uint64_t, true>>(a, b, table, 14);
+	GenerateRow<dd>(a, b, table, 15);
+	GenerateRow<qd>(a, b, table, 16);
 
 	// print the table
 	constexpr size_t COLUMN_WIDTH = 20;

include/universal/number/rational/math/logarithm.hpp

@@ -0,0 +1,35 @@
+#pragma once
+// logarithm.hpp: logarithm functions for rational
+//
+// 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 {
+
+// Natural logarithm of x
+template<unsigned nbits, typename bt>
+rational<nbits, bt> log(rational<nbits, bt> x) {
+	return rational<nbits, bt>(std::log(double(x)));
+}
+
+// Binary logarithm of x
+template<unsigned nbits, typename bt>
+rational<nbits, bt> log2(rational<nbits, bt> x) {
+	return rational<nbits, bt>(std::log2(double(x)));
+}
+
+// Decimal logarithm of x
+template<unsigned nbits, typename bt>
+rational<nbits, bt> log10(rational<nbits, bt> x) {
+	return rational<nbits, bt>(std::log10(double(x)));
+}
+
+// Natural logarithm of 1+x
+template<unsigned nbits, typename bt>
+rational<nbits, bt> log1p(rational<nbits, bt> x) {
+	return rational<nbits, bt>(std::log1p(double(x)));
+}
+
+}} // namespace sw::universal

include/universal/number/rational/math/minmax.hpp

@@ -0,0 +1,26 @@
+#pragma once
+// minmax.hpp: min/max functions for rational
+//
+// 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 {
+
+// minimum of two values
+template<unsigned nbits, typename bt>
+rational<nbits, bt>
+min(rational<nbits, bt> x, rational<nbits, bt> y) {
+	return rational<nbits, bt>(std::min(double(x), double(y)));
+}
+
+// maximum of two values
+template<unsigned nbits, typename bt>
+rational<nbits, bt>
+max(rational<nbits, bt> x, rational<nbits, bt> y) {
+	return rational<nbits, bt>(std::max(double(x), double(y)));
+}
+
+
+}} // namespace sw::universal

include/universal/number/rational/math/sqrt.hpp

@@ -0,0 +1,71 @@
+#pragma once
+// sqrt.hpp: sqrt functions for rational
+//
+// 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/native/ieee754.hpp>
+//#include <universal/number/lns/math/sqrt_tables.hpp>
+
+#ifndef LNS_NATIVE_SQRT
+#define LNS_NATIVE_SQRT 0
+#endif
+
+namespace sw { namespace universal {
+
+	/*
+	- Consider the function argument, x, in floating-point form, with a base
+	(or radix) B, exponent e, and a fraction, f , such that 1/B <= f < 1.
+	Then we have x = f Be. The number of bits in the exponent and
+	fraction, and the value of the base, depends on the particular floating
+	point arithmetic system chosen.
+
+	- Use properties of the elementary function to range reduce the argument
+	x to a small fixed interval.
+
+	- Use a small polynomial approximation to produce an initial estimate,
+	y0, of the function on the small interval. Such an estimate may
+	be good to perhaps 5 to 10 bits.
+
+	- Apply Newton iteration to refine the result. This takes the form yk =
+	yk?1/2 + (f /2)/yk?1. In base 2, the divisions by two can be done by
+	exponent adjustments in floating-point computation, or by bit shifting
+	in fixed-point computation.
+
+	Convergence of the Newton method is quadratic, so the number of
+	correct bits doubles with each iteration. Thus, a starting point correct
+	to 7 bits will produce iterates accurate to 14, 28, 56, ... bits. Since the
+	number of iterations is very small, and known in advance, the loop is
+	written as straight-line code.
+
+	- Having computed the function value for the range-reduced argument,
+	make whatever adjustments are necessary to produce the function value
+	for the original argument; this step may involve a sign adjustment,
+	and possibly a single multiplication and/or addition.
+	*/
+
+
+	// sqrt for arbitrary rational
+	template<unsigned nbits, typename bt>
+	inline rational<nbits, bt> sqrt(const rational<nbits, bt>& a) {
+#if RATIONAL_THROW_ARITHMETIC_EXCEPTION
+		if (a.isneg()) throw rational_negative_sqrt_arg();
+#else
+		if (a.isneg()) std::cerr << "rationalns argument to sqrt is negative: " << a << std::endl;
+#endif
+		if (a.iszero()) return a;
+		return rational<nbits, bt>(std::sqrt((double)a));  // TBD
+	}
+
+	// reciprocal sqrt
+	template<unsigned nbits, typename bt>
+	inline rational<nbits, bt> rsqrt(const rational<nbits, bt>& a) {
+		rational<nbits, bt> v = sqrt(a);
+		return v.reciprocate();
+	}
+
+	///////////////////////////////////////////////////////////////////
+	// specialized sqrt configurations
+
+}} // namespace sw::universal

include/universal/number/rational/mathlib.hpp

@@ -6,19 +6,19 @@
 //
 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
 
-/*
-#include <universal/number/rational/math/classify.hpp>
-#include <universal/number/rational/math/complex.hpp>
-#include <universal/number/rational/math/error_and_gamma.hpp>
-#include <universal/number/rational/math/exponent.hpp>
-#include <universal/number/rational/math/fractional.hpp>
-#include <universal/number/rational/math/hyperbolic.hpp>
-#include <universal/number/rational/math/hypot.hpp>
+
+//#include <universal/number/rational/math/classify.hpp>
+//#include <universal/number/rational/math/complex.hpp>
+//#include <universal/number/rational/math/error_and_gamma.hpp>
+//#include <universal/number/rational/math/exponent.hpp>
+//#include <universal/number/rational/math/fractional.hpp>
+//#include <universal/number/rational/math/hyperbolic.hpp>
+//#include <universal/number/rational/math/hypot.hpp>
 #include <universal/number/rational/math/logarithm.hpp>
 #include <universal/number/rational/math/minmax.hpp>
-#include <universal/number/rational/math/next.hpp>
-#include <universal/number/rational/math/pow.hpp>
+//#include <universal/number/rational/math/next.hpp>
+//#include <universal/number/rational/math/pow.hpp>
 #include <universal/number/rational/math/sqrt.hpp>
-#include <universal/number/rational/math/trigonometry.hpp>
-#include <universal/number/rational/math/truncate.hpp>
-*/
+//#include <universal/number/rational/math/trigonometry.hpp>
+//#include <universal/number/rational/math/truncate.hpp>
+