adding a posit oracle as a separate data type to be used to validate …

…specialized posits

CMakeLists.txt

@@ -130,6 +130,7 @@ option(BUILD_NUMBER_AREALS               "Set to ON to build static areal tests"
 option(BUILD_NUMBER_UNUM1S               "Set to ON to build static unum type 1 tests"         OFF)
 option(BUILD_NUMBER_UNUM2S               "Set to ON to build static unum type 2 tests"         OFF)
 option(BUILD_NUMBER_POSITS               "Set to ON to build static unum type 3 posit tests"   OFF)
+option(BUILD_NUMBER_POSITOS              "Set to ON to build static unum type 3 posito tests"  OFF)
 option(BUILD_NUMBER_VALIDS               "Set to ON to build static unum type 3 valid tests"   OFF)
 option(BUILD_NUMBER_LNS                  "Set to ON to build static lns tests"                 OFF)
 option(BUILD_NUMBER_DBNS                 "Set to ON to build static dbns tests"                OFF)
@@ -639,6 +640,7 @@ if(BUILD_NUMBER_STATICS)
 	set(BUILD_NUMBER_UNUM1S ON)
 	set(BUILD_NUMBER_UNUM2S ON)
 	set(BUILD_NUMBER_POSITS ON)
+	set(BUILD_NUMBER_POSITOS ON)
 	set(BUILD_NUMBER_VALIDS ON)
 	set(BUILD_NUMBER_LNS ON)
 	set(BUILD_NUMBER_DBNS ON)
@@ -763,6 +765,10 @@ add_subdirectory("static/posit/specialized")
 add_subdirectory("static/posit2")
 endif(BUILD_NUMBER_POSITS)
 
+if(BUILD_NUMBER_POSITOS)
+add_subdirectory("static/posito")
+endif(BUILD_NUMBER_POSITOS)
+
 if(BUILD_NUMBER_VALIDS)
 add_subdirectory("static/valid")
 endif(BUILD_NUMBER_VALIDS)

include/universal/number/posit/manipulators.hpp

@@ -24,7 +24,7 @@ namespace sw { namespace universal {
 		return str.str();
 	}
 
-	// Generate a type field descriptor for this cfloat
+	// Generate a type field descriptor for this posit
 	template<typename PositType,
 		std::enable_if_t< is_posit<PositType>, bool> = true
 	>
@@ -118,7 +118,7 @@ std::string component_values_to_string(const posit<nbits, es>& p) {
 	return str.str();
 }
 
-// generate a binary string for cfloat
+// generate a binary string for posit
 template<unsigned nbits, unsigned es>
 inline std::string to_hex(const posit<nbits, es>& v, bool nibbleMarker = false, bool hexPrefix = true) {
 	char hexChar[16] = {

include/universal/number/posit/posit.hpp

@@ -36,8 +36,12 @@
 #if !defined(POSIT_THROW_ARITHMETIC_EXCEPTION)
 // default is to use NaR as a signalling error
 #define POSIT_THROW_ARITHMETIC_EXCEPTION 0
+#if !defined(VALUE_THROW_ARITHMETIC_EXCEPTION)
 #define VALUE_THROW_ARITHMETIC_EXCEPTION 0
+#endif
+#if !defined(BITBLOCK_THROW_ARITHMETIC_EXCEPTION)
 #define BITBLOCK_THROW_ARITHMETIC_EXCEPTION 0
+#endif
 #else
 // for the composite value<> class assume the same behavior as requested for posits
 #define VALUE_THROW_ARITHMETIC_EXCEPTION POSIT_THROW_ARITHMETIC_EXCEPTION

include/universal/number/posito/ReadMe.txt

@@ -0,0 +1,86 @@
+========================================================================
+    Posit Project Overview
+========================================================================
+
+Random 32 bits which is intentionally zero heavy
+-----------------
+Sign
+|
+|     
+v    
+1100 1000 1100 0100 0100 0010 0001 0001
+
+How to decode this?
+
+Here is some pseudocode that should work on any bit sequence. Suppose
+
+  n     is the number of bits (which in the example is 32),
+  es    is the number of exponent bits (es = 3 in the example since n = 32),
+  p     is the bit sequence treated as a signed integer.
+
+Also suppose, though it is not the convention, that the bits are numbered from left to right starting with 1, 
+so the sign bit is bit 1 and the rightmost (least significant) bit is bit n. That numbering makes the pseudocode 
+more readable. Assume we have a quick machine instruction FindFirstOne(p) that finds the position of the first 1 bit 
+in p, returning n+1 if all bits are zero, and a similar instruction FindFirstZero(p). The function Bits(j..k, p) 
+extracts bits j through k of p, clipped to bit n. If j is greater than n, Bits(j..k, p) returns zero.
+
+We seek x, the value represented by p. The green text indicates what would happen for the 32 bits in the example.
+
+(* Check for exception values: *)
+If Bits(2..n, p) == 0, (* Check for a string of all zeros after the first bit. *)
+  If Bits(1..1, p) == 0, (* If the first bit is also zero, then x is zero. *)
+    x ← 0;
+  Else, (* The first bit must be one, indicating x is unsigned infinity. *)
+    x ← ±∞;
+  End If;
+  End; (* Can cut out early. *)
+End If;
+
+The example sequence fails the exception value tests, so continue.
+
+(* Process sign bit: *)
+If Bits(1..1, p) == 1, sign ← –1; p ← –p; (* 2's complement negate p. *)
+  Else sign ← 1;
+End If;
+
+The first bit is a 1, so sign = -1 and p gets negated. Flip all the bits and add 1:
+0011 0111 0011 1011 1011 1101 1110 1111
+
+(* Process regime bits. The second bit, r, determines the sign of the regime. *)
+r ← Bits(2..2, p)
+Bits(1..1, p) ← r (* Duplicate first regime bit into sign bit so the "FindFirst..." instructions work. *)
+If r == 0, (* regime value is negative *)
+  k = FindFirstOne(p);
+  shift ← –(2^es) × (k – 2); (* Actually a shift of k-2 by a shift by es, so this is fast. *)
+Else (* regime value is zero or positive *)
+  k = FindFirstZero(p);
+  shift ← (2^es) × (k – 3);
+End If;
+
+The example 001... has r = 0 so the regime value is negative. The first 1 occurs in bit 3, 
+so the regime represents –1 and contributes a shift of –(2^es) = –8. This leaves k pointing to bit 3. 
+The exponent will be in the next es bits after bit 3, which will fine-tune the shift roughed out by the regime bits.
+
+(* Add the exponent bits, as an unsigned integer, to the shift. *)
+shift ← shift + Bits(k+1..k+es, p);
+
+The three exponent bits in the example are the underlined ones: 0011 0111 0011 1011 1011 1101 1110 1111. 
+As an unsigned binary integer, 101 means decimal 5, so the shift is –8 + 5 = –3.
+
+(* The remaining bits are the fraction, with a hidden bit of 1, so we're done. It looks like a regular float: *)
+x ← sign × 2^shift × 1.Bits(k+es+1..n) 
+End
+
+The fraction bits in the example are underlined: 0011 0111 0011 1011 1011 1101 1110 1111. That fraction, 
+plus the hidden bit 1, is 121355759/(2^26), or exactly 1.80834172666072845458984375. 
+Scale this by 2^(shift) = 2^(–3) = 1/8, and the sign, to get –0.22604271583259105682373046875. 
+Notice that there are three more bits in the fraction than there are in an IEEE 32-bit float, 
+giving almost a full decimal more accuracy.
+
+I didn't actually run the above pseudocode so it may contain a typo or two, but it is closely based on the 
+Mathematica code for decoding a posit and that has been debugged and in use for months. 
+It is possible to decode the bit string _without_ doing 2's complement negation if the sign bit is zero 
+(the hidden bit for negative posits represents –2, not 1!), which according to Isaac Yonemoto leads to 
+simpler hardware... but it creates a more cryptic pseudocode explanation, so I used the 2's complement negation here.
+
+/////////////////////////////////////////////////////////////////////////////

include/universal/number/posito/attributes.hpp

@@ -0,0 +1,28 @@
+#pragma once
+// attributes.hpp: functions to query number system attributes
+//
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+//
+// 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 properties of a posit configuration
+
+using namespace sw::universal::internal;
+
+// calculate the scale of a posit
+template<unsigned nbits, unsigned es>
+inline int scale(const posito<nbits, es>& p) {
+	regime<nbits, es>    _regime;
+	exponent<nbits, es>  _exponent;
+	internal::bitblock<nbits> tmp(p.get());
+	tmp = sign(p) ? internal::twos_complement(tmp) : tmp;
+	int k = decode_regime(tmp);
+	unsigned nrRegimeBits = _regime.assign_regime_pattern(k);	// get the regime bits
+	_exponent.extract_exponent_bits(tmp, nrRegimeBits);							// get the exponent bits
+	// return the scale
+	return _regime.scale() + _exponent.scale();
+}
+
+}} // namespace sw::universal

include/universal/number/posito/exceptions.hpp

@@ -0,0 +1,74 @@
+#pragma once
+// exceptions.hpp: exception hierarchy for exceptions during posit calculations
+//
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+//
+// 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 {
+
+///////////////////////////////////////////////////////////////////////////////////////////////////
+/// POSIT ARITHMETIC EXCEPTIONS
+
+// base class for posit arithmetic exceptions
+struct posito_arithmetic_exception : public universal_arithmetic_exception {			
+	posito_arithmetic_exception(const std::string& error) 
+		: universal_arithmetic_exception(std::string("posit arithmetic exception: ") + error) {};
+};
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+/// specialized exceptions to aid application level exception handling
+
+// not_a_real is thrown when a rvar is NaR
+struct posito_nar : public posito_arithmetic_exception {
+	posito_nar() : posito_arithmetic_exception("nar (not a real)") {}
+};
+
+// divide_by_zero is thrown when the denominator in a division operator is 0
+struct posito_divide_by_zero : public posito_arithmetic_exception {
+	posito_divide_by_zero() : posito_arithmetic_exception("divide by zero") {}
+};
+
+// divide_by_nar is thrown when the denominator in a division operator is NaR
+struct posito_divide_by_nar : public posito_arithmetic_exception {
+	posito_divide_by_nar() : posito_arithmetic_exception("divide by nar") {}
+};
+
+// numerator_is_nar is thrown when the numerator in a division operator is NaR
+struct posito_numerator_is_nar : public posito_arithmetic_exception {
+	posito_numerator_is_nar() : posito_arithmetic_exception("numerator is nar") {}
+};
+
+// operand_is_nar is thrown when an rvar in a binary operator is NaR
+struct posito_operand_is_nar : public posito_arithmetic_exception {
+	posito_operand_is_nar() : posito_arithmetic_exception("operand is nar") {}
+};
+
+// thrown when division yields no signficant fraction bits
+struct posito_division_result_is_zero : public posito_arithmetic_exception {
+	posito_division_result_is_zero() : posito_arithmetic_exception("division yielded no significant bits") {}
+};
+
+// thrown when division yields NaR
+struct posito_division_result_is_infinite : public posito_arithmetic_exception {
+	posito_division_result_is_infinite() : posito_arithmetic_exception("division yielded infinite") {}
+};
+
+///////////////////////////////////////////////////////////////////////////////////////////////////
+/// POSIT INTERNAL OPERATION EXCEPTIONS
+
+struct posito_internal_exception	: public universal_internal_exception {
+	posito_internal_exception(const std::string& error) 
+		: universal_internal_exception(std::string("posit internal exception: ") + error) {};
+};
+
+struct posito_hpos_too_large : public posito_internal_exception {
+	posito_hpos_too_large() : posito_internal_exception("position of hidden bit too large for given posit") {}
+};
+
+struct posito_rbits_too_large : public posito_internal_exception {
+	posito_rbits_too_large() :posito_internal_exception("number of remaining bits too large for this fraction") {}
+};
+
+}} // namespace sw::universal

include/universal/number/posito/fdp.hpp

@@ -0,0 +1,87 @@
+#pragma once
+// fdp.hpp :  include file containing templated C++ interfaces to fused dot product
+//
+// Copyright (C) 2017 Stillwater Supercomputing, Inc.
+//
+// This file is part of the universal numbers project, which is released under an MIT Open Source license.
+#include <iostream>
+#include <vector>
+#include <universal/traits/posito_traits.hpp>
+
+namespace sw { namespace universal {
+
+/// //////////////////////////////////////////////////////////////////
+/// fused dot product operators
+/// fdp_qc         fused dot product with quire continuation
+/// fdp_stride     fused dot product with non-negative stride
+/// fdp            fused dot product of two vectors
+
+// Fused dot product with quire continuation
+template<typename Qy, typename Vector>
+void fdp_qc(Qy& sum_of_products, size_t n, const Vector& x, size_t incx, const Vector& y, size_t incy) {
+	size_t ix, iy;
+	for (ix = 0, iy = 0; ix < n && iy < n; ix = ix + incx, iy = iy + incy) {
+		sum_of_products += sw::universal::quire_mul(x[ix], y[iy]);
+	}
+}
+
+// Resolved fused dot product, with the option to control capacity bits in the quire
+template<typename Vector>
+enable_if_posito<value_type<Vector>, value_type<Vector> > // as return type
+fdp_stride(size_t n, const Vector& x, size_t incx, const Vector& y, size_t incy) {
+	constexpr unsigned nbits = Vector::value_type::nbits;
+	constexpr unsigned es = Vector::value_type::es;
+	constexpr unsigned capacity = 20; // support vectors up to 1M elements
+	quire<nbits, es, capacity> q = 0;
+	size_t ix, iy;
+	for (ix = 0, iy = 0; ix < n && iy < n; ix = ix + incx, iy = iy + incy) {
+		q += sw::universal::quire_mul(x[ix], y[iy]);
+		if (sw::universal::_trace_quire_add) std::cout << q << '\n';
+	}
+	typename Vector::value_type sum;
+	convert(q.to_value(), sum);     // one and only rounding step of the fused-dot product
+	return sum;
+}
+
+#if defined(_MSC_VER)
+/* Microsoft Visual Studio. --------------------------------- */
+// Specialized resolved fused dot product that assumes unit stride and a standard vector
+template<typename Vector>
+enable_if_posito<value_type<Vector>, value_type<Vector> > // as return type
+fdp(const Vector& x, const Vector& y) {
+	constexpr unsigned nbits = Vector::value_type::nbits;
+	constexpr unsigned es = Vector::value_type::es;
+	constexpr unsigned capacity = 20; // support vectors up to 1M elements
+	quire<nbits, es, capacity> q(0);
+	size_t ix, iy, n = size(x);
+	for (ix = 0, iy = 0; ix < n && iy < n; ++ix, ++iy) {
+		q += sw::universal::quire_mul(x[ix], y[iy]);
+	}
+	typename Vector::value_type sum;
+	convert(q.to_value(), sum);     // one and only rounding step of the fused-dot product
+	return sum;
+}
+#else
+
+// template<typename Scalar> constexpr auto size(const std::vector<Scalar>& v) -> decltype(v.size()) { return (v.size()); }
+
+// Specialized resolved fused dot product that assumes unit stride and a standard vector
+template<typename Vector>
+enable_if_posito<value_type<Vector>, value_type<Vector> > // as return type
+fdp(const Vector& x, const Vector& y) {
+	constexpr unsigned nbits = Vector::value_type::nbits;
+	constexpr unsigned es = Vector::value_type::es;
+	constexpr unsigned capacity = 20; // support vectors up to 1M elements
+	quire<nbits, es, capacity> q(0);
+	size_t ix, iy, n = size(x);
+	for (ix = 0, iy = 0; ix < n && iy < n; ++ix, ++iy) {
+		q += sw::universal::quire_mul(x[ix], y[iy]);
+	}
+	typename Vector::value_type sum;
+	convert(q.to_value(), sum);     // one and only rounding step of the fused-dot product
+	return sum;
+}
+#endif
+
+}} // namespace sw::universal
+