initial check-in of quad-double skeleton

include/universal/number/qd/attributes.hpp

@@ -0,0 +1,65 @@
+#pragma once
+// attributes.hpp: information functions for quad-double floating-point type and value attributes
+//
+// 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 {
+
+// functions to provide details about the properties of a quad-double (qd) configuration
+	inline bool sign(const qd& a) {
+		return a.sign();
+	}
+
+	inline int scale(const qd& a) {
+		return a.scale();
+	}
+
+	// generate the maxneg through maxpos value range of a quad-double configuration
+	std::string qd_range() {
+		qd v;
+		std::stringstream s;
+		s << std::setw(80) << type_tag(v) << " : [ "
+			<< v.maxneg() << " ... "
+			<< v.minneg() << " "
+			<< "0 "
+			<< v.minpos() << " ... "
+			<< v.maxpos() << " ]";
+		return s.str();
+	}
+
+	/*
+	// report dynamic range of a type, specialized for a quad-double
+	std::string dynamic_range(const qd& a) {
+		std::stringstream s;
+		qd b(SpecificValue::maxneg), c(SpecificValue::minneg), d(SpecificValue::minpos), e(SpecificValue::maxpos);
+		s << type_tag(a) << ": ";
+		s << "minpos scale " << std::setw(10) << d.scale() << "     ";
+		s << "maxpos scale " << std::setw(10) << e.scale() << '\n';
+		s << "[" << b << " ... " << c << ", -0, +0, " << d << " ... " << e << "]\n";
+		s << "[" << to_binary(b) << " ... " << to_binary(c) << ", -0, +0, " << to_binary(d) << " ... " << to_binary(e) << "]\n";
+		qd ninf(SpecificValue::infneg), pinf(SpecificValue::infpos);
+		s << "inclusive range = (" << to_binary(ninf) << ", " << to_binary(pinf) << ")\n";
+		s << "inclusive range = (" << ninf << ", " << pinf << ")\n";
+		return s.str();
+	}
+	*/
+
+	int minpos_scale(const qd& b) {
+		qd c(b);
+		return c.minpos().scale();
+	}
+
+	int maxpos_scale(const qd& b) {
+		qd c(b);
+		return c.maxpos().scale();
+	}
+
+	int max_negative_scale(const qd& b) {
+		qd c(b);
+		return c.maxneg().scale();
+	}
+
+}} // namespace sw::universal

include/universal/number/qd/exceptions.hpp

@@ -0,0 +1,64 @@
+#pragma once
+// exceptions.hpp: definition of arbitrary configuration quad-double exceptions
+//
+// 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/common/exceptions.hpp>
+
+namespace sw { namespace universal {
+
+// base class for quad-double arithmetic exceptions
+struct qd_arithmetic_exception : public universal_arithmetic_exception {
+	qd_arithmetic_exception(const std::string& err) : universal_arithmetic_exception(std::string("quad-double arithmetic exception: ") + err) {};
+};
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+/// specialized exceptions to aid application level exception handling
+
+// invalid_argument is thrown when a mathematical function argument is invalid
+struct qd_invalid_argument : public qd_arithmetic_exception {
+	qd_invalid_argument() : qd_arithmetic_exception("invalid argument") {}
+};
+
+// not_a_number is thrown when a rvar is NaN
+struct qd_not_a_number : public qd_arithmetic_exception {
+	qd_not_a_number() : qd_arithmetic_exception("not a number") {}
+};
+
+// divide by zero arithmetic exception for reals
+struct qd_divide_by_zero : public qd_arithmetic_exception {
+	qd_divide_by_zero() : qd_arithmetic_exception("divide by zero") {}
+};
+
+// divide_by_nan is thrown when the denominator in a division operator is NaN
+struct qd_divide_by_nan : public qd_arithmetic_exception {
+	qd_divide_by_nan() : qd_arithmetic_exception("divide by nan") {}
+};
+
+// operand_is_nan is thrown when an rvar in a binary operator is NaN
+struct qd_operand_is_nan : public qd_arithmetic_exception {
+	qd_operand_is_nan() : qd_arithmetic_exception("operand is nan") {}
+};
+
+// negative argument to sqrt
+struct qd_negative_sqrt_arg : public qd_arithmetic_exception {
+	qd_negative_sqrt_arg() : qd_arithmetic_exception("negative sqrt argument") {}
+};
+
+// negative argument to nroot
+struct qd_negative_nroot_arg : public qd_arithmetic_exception {
+	qd_negative_nroot_arg() : qd_arithmetic_exception("negative nroot argument") {}
+};
+
+
+///////////////////////////////////////////////////////////////////////////////////////////////////
+/// REAL INTERNAL OPERATION EXCEPTIONS
+
+struct qd_internal_exception : public universal_internal_exception {
+	qd_internal_exception(const std::string& err) : universal_internal_exception(std::string("quad-double internal exception: ") + err) {};
+};
+
+
+}} // namespace sw::universal

include/universal/number/qd/manipulators.hpp

@@ -0,0 +1,43 @@
+// manipulators.hpp: definitions of helper functions for quad-double type manipulation
+//
+// 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 <string>
+#include <iomanip>
+#include <universal/internal/blockbinary/blockbinary.hpp>
+#include <universal/number/qd/qd_fwd.hpp>
+// pull in the color printing for shells utility
+#include <universal/utility/color_print.hpp>
+
+namespace sw { namespace universal {
+
+	// Generate a type tag for a quad-double
+	inline std::string type_tag(const qd& = {}) {
+		return std::string("quad-double");
+	}
+
+	// generate a binary, color-coded representation of the quad-double
+	inline std::string color_print(qd const& r, bool nibbleMarker = false) {
+		//constexpr unsigned es = 11;
+		//constexpr unsigned fbits = 106;
+		std::stringstream s;
+
+		/*
+		Color red(ColorCode::FG_RED);
+		Color yellow(ColorCode::FG_YELLOW);
+		Color blue(ColorCode::FG_BLUE);
+		Color magenta(ColorCode::FG_MAGENTA);
+		Color cyan(ColorCode::FG_CYAN);
+		Color white(ColorCode::FG_WHITE);
+		Color def(ColorCode::FG_DEFAULT);
+		*/
+		for (int i = 0; i < 4; ++i) {
+			s << color_print(r[i], nibbleMarker);
+			if (i < 3) s << ", ";
+		}
+		return s.str();
+	}
+
+}} // namespace sw::universal

include/universal/number/qd/math/classify.hpp

@@ -0,0 +1,55 @@
+#pragma once
+// classify.hpp: classification functions for quad-double (qd) 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 {
+
+// STD LIB function for IEEE floats: Categorizes floating point value arg into the following categories: zero, subnormal, normal, infinite, NAN, or implementation-defined category.
+int fpclassify(const qd& a) {
+	return (std::fpclassify(a.high()));
+}
+
+// STD LIB function for IEEE floats: Determines if the given floating point number arg is a cfloative or negative infinity.
+// specialized for quad-double (qd)
+inline bool isinf(const qd& a) {
+	return (std::fpclassify(a.high()) == FP_INFINITE);
+}
+
+// STD LIB function for IEEE floats: Determines if the given floating point number arg is a not-a-number (NaN) value.
+// specialized for quad-double (qd)
+inline bool isnan(const qd& a) {
+	return (std::fpclassify(a.high()) == FP_NAN);
+}
+
+// STD LIB function for IEEE floats: Determines if the given floating point number arg has finite value i.e. it is normal, subnormal or zero, but not infinite or NaN.
+// specialized for quad-double (qd)
+inline bool isfinite(const qd& a) {
+	return (std::fpclassify(a.high()) != FP_INFINITE) && (std::fpclassify(a.high()) != FP_NAN);
+}
+
+// STD LIB function for IEEE floats: Determines if the given floating point number arg is normal, i.e. is neither zero, subnormal, infinite, nor NaN.
+// specialized for quad-double (qd)
+inline bool isnormal(const qd& a) {
+	return (std::fpclassify(a.high()) == FP_NORMAL);
+}
+
+// STD LIB function for IEEE floats: Determines if the given floating point number arg is denormal, i.e. is neither zero, normal, infinite, nor NaN.
+// specialized for quad-double (qd)
+inline bool isdenorm(const qd& a) {
+	return (std::fpclassify(a.high()) == FP_SUBNORMAL);
+}
+
+inline bool iszero(const qd& a) {
+	return (std::fpclassify(a.high()) == FP_ZERO);
+}
+
+bool signbit(qd const& a) {
+	auto signA = std::copysign(1.0, a.high());
+	return signA < 0.0;
+}
+
+}} // namespace sw::universal

include/universal/number/qd/math/complex.hpp

@@ -0,0 +1,12 @@
+#pragma once
+// complex.hpp: complex support 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 {
+
+
+}} // namespace sw::universal

include/universal/number/qd/math/error_and_gamma.hpp

@@ -0,0 +1,21 @@
+#pragma once
+// error_and_gamma.hpp: error/gamma function support for quad-double (qd) 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 {
+
+	// Compute the error function erf(x) = 2 over sqrt(PI) times Integral from 0 to x of e ^ (-t)^2 dt
+	qd erf(qd x) {
+		return qd(std::erf(double(x.high())));
+	}
+
+	// Compute the complementary error function: 1 - erf(x)
+	qd erfc(qd x) {
+		return qd(std::erfc(double(x.high())));
+	}
+
+}} // namespace sw::universal

include/universal/number/qd/math/exponent.hpp

@@ -0,0 +1,105 @@
+#pragma once
+// exponent.hpp: exponent functions for quad-double (qd) floating-point
+//
+// algorithms courtesy Scibuilders, Jack Poulson
+// 
+// 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 {
+
+	// fwd reference
+	qd ldexp(qd const&, int);
+
+static const int n_inv_fact = 15;
+static const double inv_fact[n_inv_fact][2] = {
+  { 1.66666666666666657e-01,  9.25185853854297066e-18},
+  { 4.16666666666666644e-02,  2.31296463463574266e-18},
+  { 8.33333333333333322e-03,  1.15648231731787138e-19},
+  { 1.38888888888888894e-03, -5.30054395437357706e-20},
+  { 1.98412698412698413e-04,  1.72095582934207053e-22},
+  { 2.48015873015873016e-05,  2.15119478667758816e-23},
+  { 2.75573192239858925e-06, -1.85839327404647208e-22},
+  { 2.75573192239858883e-07,  2.37677146222502973e-23},
+  { 2.50521083854417202e-08, -1.44881407093591197e-24},
+  { 2.08767569878681002e-09, -1.20734505911325997e-25},
+  { 1.60590438368216133e-10,  1.25852945887520981e-26},
+  { 1.14707455977297245e-11,  2.06555127528307454e-28},
+  { 7.64716373181981641e-13,  7.03872877733453001e-30},
+  { 4.77947733238738525e-14,  4.39920548583408126e-31},
+  { 2.81145725434552060e-15,  1.65088427308614326e-31}
+};
+
+// Base-e exponential function
+qd exp(const qd& a) {
+	/* Strategy:  We first reduce the size of x by noting that
+     
+		exp(kr + m * log(2)) = 2^m * exp(r)^k
+
+	where m and k are integers.  By choosing m appropriately
+	we can make |kr| <= log(2) / 2 = 0.347.  Then exp(r) is 
+	evaluated using the familiar Taylor series.  Reducing the 
+	argument substantially speeds up the convergence.       */  
+
+	const double k = 512.0;
+	const double inv_k = 1.0 / k;
+
+	if (a.high() <= -709.0) return qd(0.0);
+
+	if (a.high() >=  709.0) return qd(SpecificValue::infpos);
+
+	if (a.iszero()) return qd(1.0);
+
+	if (a.isone()) return qd_e;
+
+	double m = std::floor(a.high() / qd_log2.high() + 0.5);
+	qd r = mul_pwr2(a - qd_log2 * m, inv_k);
+	qd s, t, p;
+
+	p = sqr(r);
+	s = r + mul_pwr2(p, 0.5);
+	p *= r;
+	t = p * qd(inv_fact[0][0], inv_fact[0][1]);
+	int i = 0;
+	do {
+		s += t;
+		p *= r;
+		++i;
+		t = p * qd(inv_fact[i][0], inv_fact[i][1]);
+	} while (std::abs(double(t)) > inv_k * d_eps && i < 5);
+
+	s += t;
+
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s = mul_pwr2(s, 2.0) + sqr(s);
+	s += 1.0;
+
+	return ldexp(s, static_cast<int>(m));
+}
+
+// Base-2 exponential function
+qd exp2(qd x) {
+	return qd(std::exp2(double(x)));
+}
+
+// Base-10 exponential function
+qd exp10(qd x) {
+	return qd(std::exp10(double(x)));
+}
+
+// Base-e exponential function exp(x)-1
+qd expm1(qd x) {
+	return qd(std::expm1(double(x)));
+}
+
+
+}} // namespace sw::universal

include/universal/number/qd/math/fractional.hpp

@@ -0,0 +1,21 @@
+#pragma once
+// fractional.hpp: fractional support for quad-double (qd) 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 {
+
+	// fmod retuns x - n*y where n = x/y with the fractional part truncated
+	qd fmod(qd x, qd y) {
+		return qd(std::fmod(double(x), double(y)));
+	}
+
+	// shim to stdlib
+	qd remainder(qd x, qd y) {
+		return qd(std::remainder(double(x), double(y)));
+	}
+
+}} // namespace sw::universal