WIP: precision experiments with rationals

include/universal/internal/blockbinary/blockbinary.hpp

@@ -170,11 +170,11 @@ class blockbinary {
 	explicit operator unsigned long() const      { return (unsigned long)to_ull(); }
 	explicit operator unsigned long long() const { return to_ull(); }
 	// TODO: these need proper implementations that can convert very large integers to the proper scale afforded by the floating-point formats
-	explicit operator float() const              { return float(to_long_long()); }
-	explicit operator double() const             { return double(to_long_long()); }
+	explicit operator float() const              { return to_native<float>(); }
+	explicit operator double() const             { return to_native<double>(); }
 
 #if LONG_DOUBLE_SUPPORT
-	explicit operator long double() const        { return (long double)to_long_long(); }
+	explicit operator long double() const        { return to_native<long double>(); }
 #endif
 
 	// limb access operators
@@ -566,7 +566,11 @@ class blockbinary {
 		_block[MSU] &= MSU_MASK; // assert precondition of properly nulled leading non-bits
 		return *this;
 	}
-	constexpr blockbinary& twosComplement() noexcept { // in-place 2's complement
+	/// <summary>
+	/// in-place 2's complement
+	/// </summary>
+	/// <returns>2's complement of original</returns>
+	constexpr blockbinary& twosComplement() noexcept {
 		blockbinary plusOne(1);
 		if constexpr (NumberType == BinaryNumberType::Signed) {
 			flip();
@@ -780,7 +784,19 @@ class blockbinary {
 		}
 		return ull;
 	}
+	template<typename Real,
+		typename = typename std::enable_if< std::is_floating_point<Real>::value, Real >::type>
+	Real to_native() const {
+		blockbinary tmp(*this);
+		if (isneg()) tmp.twosComplement();
+		Real v{ 1.0 }, base{ 1.0 };
+		for (unsigned i = 0; i < nbits - 1; ++i) {
+			if (tmp.test(i)) v *= base;
+			base *= 2.0;
+		}
 
+		return (isneg() ? -v : v);
+	}
 	// determine the rounding mode: result needs to be rounded up if true
 	bool roundingMode(unsigned targetLsb) const {
 		bool lsb = at(targetLsb);

include/universal/number/cfloat/cfloat_impl.hpp

@@ -467,6 +467,13 @@ class cfloat {
 	BIT_CAST_CONSTEXPR cfloat& operator=(float rhs)     noexcept { return convert_ieee754(rhs); }
 	BIT_CAST_CONSTEXPR cfloat& operator=(double rhs)    noexcept { return convert_ieee754(rhs); }
 
+	// make conversions to native types explicit
+	explicit operator int()                       const noexcept { return to_int(); }
+	explicit operator long()                      const noexcept { return to_long(); }
+	explicit operator long long()                 const noexcept { return to_long_long(); }
+	explicit operator float()                     const noexcept { return to_native<float>(); }
+	explicit operator double()                    const noexcept { return to_native<double>(); }
+
 	// guard long double support to enable ARM and RISC-V embedded environments
 #if LONG_DOUBLE_SUPPORT
 	explicit operator long double()               const noexcept { return to_native<long double>(); }
@@ -1905,13 +1912,6 @@ class cfloat {
 		return v;
 	}
 
-	// make conversions to native types explicit
-	explicit operator int()       const noexcept { return to_int(); }
-	explicit operator long()      const noexcept { return to_long(); }
-	explicit operator long long() const noexcept { return to_long_long(); }
-	explicit operator float()     const noexcept { return to_native<float>(); }
-	explicit operator double()    const noexcept { return to_native<double>(); }
-
 	// convert a cfloat to a blocktriple with the fraction format 1.ffff
 	// we are using the same block type so that we can use block copies to move bits around.
 	// Since we tend to have at least two exponent bits, this will lead to

include/universal/number/rational/numeric_limits.hpp

@@ -35,7 +35,7 @@ class numeric_limits< sw::universal::rational<nbits,bt> >
 		// the pattern is this ratio 0.....1 / 01.....0
 		RationalType r(0, 0);
 		r.setnbit(0);
-		r.setbits(nbits - 2);
+		r.setdbit(nbits - 2);
 		return r;
 	}
 	static constexpr RationalType round_error() { // return largest rounding error

include/universal/number/rational/rational_impl.hpp

@@ -219,6 +219,8 @@ class rational {
 	constexpr void setbits(std::int64_t bits) noexcept { n = bits; d = 1; }
 	constexpr void setnbit(unsigned index) noexcept { n.set(index); }
 	constexpr void setdbit(unsigned index) noexcept { d.set(index); }
+	constexpr void resetnbit(unsigned index) noexcept { n.reset(index); }
+	constexpr void resetdbit(unsigned index) noexcept { d.reset(index); }
 
 	// create specific number system values of interest
 	constexpr rational& maxpos() noexcept {

static/rational/api/precision.cpp

@@ -0,0 +1,88 @@
+// precision.cpp: characterization of rational precision as a function of size
+//
+// 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>
+
+// minimum set of include files to reflect source code dependencies
+// Configure the rational template environment
+// enable/disable arithmetic exceptions
+#define RATIONAL_THROW_ARITHMETIC_EXCEPTION 0
+#include <universal/number/rational/rational.hpp>
+#include <universal/number/cfloat/cfloat.hpp>
+#include <universal/number/posit/posit.hpp>
+#include <universal/verification/test_suite.hpp>
+
+namespace sw {
+	namespace universal {
+
+		template<unsigned nbits, unsigned es>
+		void epsilon() {
+			std::cout << std::setw(5) << nbits << "\t" 
+				<< std::setw(15) << std::numeric_limits<cfloat<nbits, es>>::epsilon() << "\t"
+				<< std::setw(15) << std::numeric_limits<posit<nbits, 2>>::epsilon() << "\t"
+				<< std::setw(15) << std::numeric_limits<rational<nbits>>::epsilon() << '\n';
+		}
+	}
+}
+
+
+int main()
+try {
+	using namespace sw::universal;
+
+	std::string test_suite  = "rational precision characterization";
+	std::string test_tag    = "precision";
+	bool reportTestCases    = true;
+	int nrOfFailedTestCases = 0;
+
+	ReportTestSuiteHeader(test_suite, reportTestCases);
+
+	// what is the progression of precision for increasingly larger binary rational types
+	std::cout << "epsilon for different sizes\n"
+		<< std::setw(5) << "nbits" 
+		<< std::setw(15) << "cfloat"
+		<< std::setw(15) << "posit"
+		<< std::setw(15) << "rational\n";
+	epsilon<4,2>();
+	epsilon<8,2>();
+	epsilon<12,5>();
+	epsilon<16,5>();
+	epsilon<20,8>();
+	epsilon<24,8>();
+	epsilon<28,8>();
+	epsilon<32,8>();
+	epsilon<40,11>();
+	epsilon<48,11>();
+	epsilon<56,11>();
+	epsilon<64,11>();
+	epsilon<80,15>();
+	epsilon<96,15>();
+	epsilon<112,15>();
+	epsilon<128,15>();
+
+	ReportTestSuiteResults(test_suite, nrOfFailedTestCases);
+	return (nrOfFailedTestCases > 0 ? EXIT_FAILURE : EXIT_SUCCESS);
+}
+catch (char const* msg) {
+	std::cerr << "Caught ad-hoc exception: " << msg << std::endl;
+	return EXIT_FAILURE;
+}
+catch (const sw::universal::universal_arithmetic_exception& err) {
+	std::cerr << "Caught unexpected universal arithmetic exception : " << err.what() << std::endl;
+	return EXIT_FAILURE;
+}
+catch (const sw::universal::universal_internal_exception& err) {
+	std::cerr << "Caught unexpected universal internal exception: " << err.what() << std::endl;
+	return EXIT_FAILURE;
+}
+catch (const std::runtime_error& err) {
+	std::cerr << "Caught runtime exception: " << err.what() << std::endl;
+	return EXIT_FAILURE;
+}
+catch (...) {
+	std::cerr << "caught unknown exception" << std::endl;
+	return EXIT_FAILURE;
+}