diff --git a/include/sw/universal/number/unum2/CMakeLists.txt b/include/sw/universal/number/unum2/CMakeLists.txt deleted file mode 100644 index 38068a2e0..000000000 --- a/include/sw/universal/number/unum2/CMakeLists.txt +++ /dev/null @@ -1,3 +0,0 @@ -file (GLOB SOURCES "./*.cpp") - -compile_all("true" "unum2" "${SOURCES}") diff --git a/include/sw/universal/number/unum2/common.hpp b/include/sw/universal/number/unum2/common.hpp new file mode 100644 index 000000000..e97a18bbd --- /dev/null +++ b/include/sw/universal/number/unum2/common.hpp @@ -0,0 +1,28 @@ +// Common implementations +// +// Copyright (C) 2017-2026 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT + +#pragma once + +#include + + +// No of operations supported by the op matrix +#define OP_MATRIX_TOTAL_SUPPORTED_OPS 2 + + +static inline uint64_t _horizontal_invert(uint64_t i, uint64_t MASK) { + // Invert index according to 2's compliment. + return ((~i & MASK) + 1) & MASK; +} + + +namespace sw { namespace universal { + +enum op_matrix_type { + OP_MATRIX_TYPE_ADD = 0, + OP_MATRIX_TYPE_MUL = 1 +}; + +}} diff --git a/include/sw/universal/number/unum2/lattice.hpp b/include/sw/universal/number/unum2/lattice.hpp new file mode 100644 index 000000000..54ad5d5c6 --- /dev/null +++ b/include/sw/universal/number/unum2/lattice.hpp @@ -0,0 +1,169 @@ +// Lattice implementation for Unum 2.0 +// +// Copyright (C) 2017-2026 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + + +static inline bool _is_power_of_two(uint64_t value) { + int count = value & 0x01; + while(value = value >> 1) + count += value & 0x01; + + return count == 1; +} + + +namespace sw { namespace universal { + +template +class lattice final { +private: + std::vector _exacts = { exacts... }; + uint64_t _N = sizeof... (exacts) << 3; // Number of points including non-exacts + uint64_t _N_half = _N >> 1; + uint64_t _N_quarter = _N_half >> 1; + uint64_t _MASK = _N - 1; + +public: + lattice() { + if(!_is_power_of_two(_N)) + throw std::invalid_argument("Number of elements in the lattice has to be power of 2!"); + + // First element in the lattice must be 1. + if(_exacts[0] != 1) + throw std::invalid_argument("First element in lattice must be 1"); + + int last = 1; + for(std::vector::iterator it = _exacts.begin() + 1; it < _exacts.end(); ++it) { + if(*it < 1 || *it <= last) + throw std::invalid_argument("Lattice always be in ascending order"); + last = *it; + } + } + + lattice(const lattice&) = delete; + lattice& operator=(const lattice&) = delete; + + std::string get_exact(uint64_t i) const { + if(i >= _N) + throw std::out_of_range("Lattice index out of range"); + + // Return nothing if not exact. + if(i & 0x01) + return ""; + + // Check for infinity, zero, -1 or 1. + if(i == _N_half) + return "inf"; + else if(i == 0) + return "0"; + else if(i == _N_quarter) + return "1"; + else if(i == 3 * _N_quarter) + return "-1"; + + std::ostringstream oss; + + // Negative + if(i > _N_half) { + oss << '-'; + i = _horizontal_invert(i, _MASK); + } + + // Vertical invert + if(i >= _N_quarter) + oss << _exacts[(i - _N_quarter) >> 1]; + else oss << '/' << _exacts[_exacts.size() - (i >> 1)]; + + return oss.str(); + } + + void print() const { + std::cout << "inf <-->"; + + size_t size = _exacts.size(); + + for(int i = size - 1; ~i; i--) + std::cout << " -" << _exacts[i] << " <-->"; + for(int i = 1; i < size; i++) + std::cout << " -/" << _exacts[i] << " <-->"; + + std::cout << " 0 <-->"; + + for(int i = size - 1; i; i--) + std::cout << " /" << _exacts[i] << " <-->"; + for(int i = 0; i < size; i++) + std::cout << " " << _exacts[i] << " <-->"; + + std::cout << " inf" << std::endl; + } + + double exactvalue(uint64_t i) const { + if(i >= _N) + throw std::out_of_range("Lattice index out of range"); + + // Return nothing if not exact. + if(i & 0x01) + return 0; + + // Check for infinity, zero, -1 or 1. + if(i == _N >> 1) + return INFINITY; + else if(i == 0) + return 0; + else if(i == _N_half) + return 1; + else if(i == 3 * _N_quarter) + return -1; + + double res = 1; + + // Negative + if(i > _N_half) { + res = -1; + i = _horizontal_invert(i, _MASK); + } + + // Vertical invert + if(i >= _N_quarter) + res *= static_cast(_exacts[(i - _N_quarter) >> 1]); + else res *= 1 / static_cast(_exacts[_exacts.size() - (i >> 1)]); + + return res; + } + + static lattice& instance() { + static lattice _lat; + return _lat; + } + + template + static op_matrix& op_matrix_instance() { + static op_matrix _op_mat = op_matrix(sizeof... (exacts) << 3); + return _op_mat; + } + + template friend class sw::universal::unum2; +}; + + +// SOME DEFAULT LATTICES + +// 8bit linear lattice, 256 points +using linear_8bit = lattice<1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, + 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32>; + +}} diff --git a/include/sw/universal/number/unum2/manipulators.hpp b/include/sw/universal/number/unum2/manipulators.hpp deleted file mode 100644 index 579db009a..000000000 --- a/include/sw/universal/number/unum2/manipulators.hpp +++ /dev/null @@ -1,69 +0,0 @@ -#pragma once -// manipulators.hpp: definitions of helper functions for unum2 type manipulation -// -// Copyright (C) 2017-2021 Stillwater Supercomputing, Inc. -// -// This file is part of the universal numbers project, which is released under an MIT Open Source license. -#include -#include -#include // for frexp/frexpf -#include // for typeid() - -// pull in the color printing for shells utility -#include - -// This file contains functions that use the unum2 type. -// If you have helper functions that the unum2 type could use, but dofsize not depend on -// the unum2 type, you can add them to the file unum_helpers.hpp. - -namespace sw { namespace universal { - -// DEBUG/REPORTING HELPERS - -template -std::string unum2_range() { - std::stringstream ss; - ss << " unum2<" << std::setw(3) << essize << "," << fsize << "> "; - //ss << "minpos scale " << std::setw(10) << minpos_scale() << " "; - //ss << "maxpos scale " << std::setw(10) << maxpos_scale() << " "; - ss << "minimum " << std::setw(12) << std::numeric_limits>::min() << " "; - ss << "maximum " << std::setw(12) << std::numeric_limits>::max() ; - return ss.str(); -} - -// Generate a type tag for this unum, for example, unum<8,1> -template -std::string type_tag(const unum2& = {}) { - std::stringstream ss; - ss << "unum2<" << essize << "," << fsize << ">"; - return ss.str(); -} - -// generate a unum2 format ASCII format essize.fsizexNN...NNp -template -inline std::string hex_print(const unum2& p) { - std::stringstream ss; - ss << essize << '.' << fsize << 'x' << to_hex(p.get()) << 'p'; - return ss.str(); -} - -template -std::string pretty_print(const unum2& p, int printPrecision = std::numeric_limits::max_digits10) { - std::stringstream ss; - return ss.str(); -} - -template -std::string info_print(const unum2& p, int printPrecision = 17) { - std::stringstream ss; - return ss.str(); -} - -template -std::string color_print(const unum2& p) { - std::stringstream ss; - return ss.str(); -} - -}} // namespace sw::universal - diff --git a/include/sw/universal/number/unum2/math_functions.hpp b/include/sw/universal/number/unum2/math_functions.hpp deleted file mode 100644 index 5b18ba457..000000000 --- a/include/sw/universal/number/unum2/math_functions.hpp +++ /dev/null @@ -1,41 +0,0 @@ -#pragma once -// math_functions.hpp: definition of arbitrary real mathematical functions -// -// Copyright (C) 2017-2021 Stillwater Supercomputing, Inc. -// -// This file is part of the universal numbers project, which is released under an MIT Open Source license. - -#if defined(__clang__) -/* Clang/LLVM. ---------------------------------------------- */ - - -#elif defined(__ICC) || defined(__INTEL_COMPILER) -/* Intel ICC/ICPC. ------------------------------------------ */ - - -#elif defined(__GNUC__) || defined(__GNUG__) -/* GNU GCC/G++. --------------------------------------------- */ - - -#elif defined(__HP_cc) || defined(__HP_aCC) -/* Hewlett-Packard C/aC++. ---------------------------------- */ - -#elif defined(__IBMC__) || defined(__IBMCPP__) -/* IBM XL C/C++. -------------------------------------------- */ - -#elif defined(_MSC_VER) -/* Microsoft Visual Studio. --------------------------------- */ - - -#elif defined(__PGI) -/* Portland Group PGCC/PGCPP. ------------------------------- */ - -#elif defined(__SUNPRO_C) || defined(__SUNPRO_CC) -/* Oracle Solaris Studio. ----------------------------------- */ - -#endif - -namespace sw { namespace universal { - - -}} // namespace sw::universal diff --git a/include/sw/universal/number/unum2/op_matrix.hpp b/include/sw/universal/number/unum2/op_matrix.hpp new file mode 100644 index 000000000..e737d867d --- /dev/null +++ b/include/sw/universal/number/unum2/op_matrix.hpp @@ -0,0 +1,89 @@ +// Operation matrix/table for unum2 operations. +// +// Copyright (C) 2017-2026 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT + +#pragma once + +#include +#include + +#include +#include + +namespace sw { namespace universal { + +template +class op_matrix final { +private: + // Buffer of op matrices. + unum2 *_matrices[OP_MATRIX_TOTAL_SUPPORTED_OPS]; + bool *_present_flags[OP_MATRIX_TOTAL_SUPPORTED_OPS]; + size_t _N; + +public: + op_matrix(size_t N) : _N { N } { + for(int i = 0; i < OP_MATRIX_TOTAL_SUPPORTED_OPS; i++) { + _matrices[i] = nullptr; + _present_flags[i] = nullptr; + } + + try { + for(int i = 0; i < OP_MATRIX_TOTAL_SUPPORTED_OPS; i++) { + _matrices[i] = new unum2[N * N]; + _present_flags[i] = new bool[N * N]; + } + } catch(const std::bad_alloc& e) { + cleanup(); + throw; // rethrow + } + + for(int i = 0; i < OP_MATRIX_TOTAL_SUPPORTED_OPS; i++) { + for(int j = 0; j < N * N; j++) + _present_flags[i][j] = false; + } + } + + op_matrix(const op_matrix&) = delete; + op_matrix& operator=(const op_matrix&) = delete; + op_matrix(op_matrix&&) = delete; + op_matrix& operator=(op_matrix&&) = delete; + + ~op_matrix() { + cleanup(); + } + + // Check if op_matrix has an operation + bool has(uint64_t i, uint64_t j, op_matrix_type type) { + return _present_flags[type][i * _N + j]; + } + + // Set a op_matrix operation + void set(uint64_t i, uint64_t j, op_matrix_type type, unum2 num) { + if(has(i, j, type)) + return; + + _matrices[type][i * _N + j] = num; + _present_flags[type][i * _N + j] = true; + } + + unum2 get(uint64_t i, uint64_t j, op_matrix_type type) { + if(!has(i, j, type)) + throw std::runtime_error("The given operation has to been set/interned yet!"); + + return _matrices[type][i * _N + j]; + } + +private: + void cleanup() { + for(int i = 0; i < OP_MATRIX_TOTAL_SUPPORTED_OPS; i++) { + delete[] _matrices[i]; + delete[] _present_flags[i]; + + _matrices[i] = nullptr; + _present_flags[i] = nullptr; + } + } +}; + +}} diff --git a/include/sw/universal/number/unum2/unum2.hpp b/include/sw/universal/number/unum2/unum2.hpp new file mode 100644 index 000000000..84b285216 --- /dev/null +++ b/include/sw/universal/number/unum2/unum2.hpp @@ -0,0 +1,8 @@ +// The Unum2 number system +// +// Copyright (C) 2017-2026 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT + +#pragma once + +#include diff --git a/include/sw/universal/number/unum2/unum2_impl.hpp b/include/sw/universal/number/unum2/unum2_impl.hpp index d75bb8f92..3c1653ed4 100644 --- a/include/sw/universal/number/unum2/unum2_impl.hpp +++ b/include/sw/universal/number/unum2/unum2_impl.hpp @@ -1,320 +1,643 @@ -#pragma once -// unum2.hpp: definition of the flexible configuration universal number system -// -// Copyright (C) 2017-2021 Stillwater Supercomputing, Inc. +// Universal Number 2.0 implementation including SORN (Sets of Real Numbers). // -// This file is part of the universal numbers project, which is released under an MIT Open Source license. -#include -#include +// Copyright (C) 2017-2026 Stillwater Supercomputing, Inc. +// SPDX-License-Identifier: MIT + +#pragma once + + +#ifndef UNUM2_USE_OP_MATRIX +#define UNUM2_USE_OP_MATRIX 1 // Use op matrix +#endif // UNUM2_USE_OP_MATRIX + -#include -#include -#include -#include +#include + +#include +#include +#include +#include +#include +#include namespace sw { namespace universal { - -// Forward definitions -template class unum2; -template unum2 abs(const unum2& v); - -// template class reprfsizesizeenting a value in scientific notation, using a template size for the number of fraction bits -template -class unum2 { -public: - static constexpr unsigned UTAGSIZE = 1 + esizesize + fsizesize; - static constexpr unsigned UTAGMASK = unsigned(~(int64_t(-1) << UTAGSIZE)); - static constexpr unsigned EBITSMASK = 1; - static constexpr unsigned FBITSMASK = 2; - - unum2() {} - - // specific value constructor - constexpr unum2(const SpecificValue code) { - switch (code) { - case SpecificValue::maxpos: - maxpos(); - break; - case SpecificValue::minpos: - minpos(); - break; - default: - zero(); - break; - case SpecificValue::minneg: - minneg(); - break; - case SpecificValue::maxneg: - maxneg(); - break; - } - } - - unum2(signed char initial_value) { *this = initial_value; } - unum2(short initial_value) { *this = initial_value; } - unum2(int initial_value) { *this = initial_value; } - unum2(long long initial_value) { *this = initial_value; } - unum2(unsigned long long initial_value) { *this = initial_value; } - unum2(float initial_value) { *this = initial_value; } - unum2(double initial_value) { *this = initial_value; } - unum2(long double initial_value) { *this = initial_value; } - unum2(const unum2& rhs) { *this = rhs; } - - // assignment operators - unum2& operator=(signed char rhs) { - return *this = (long long)(rhs); - } - unum2& operator=(short rhs) { - return *this = (long long)(rhs); - } - unum2& operator=(int rhs) { - return *this = (long long)(rhs); - } - unum2& operator=(long long rhs) { - return *this; - } - unum2& operator=(unsigned long long rhs) { - return *this; - } - unum2& operator=(float rhs) { - - return *this; - } - unum2& operator=(double rhs) { - - return *this; - } - unum2& operator=(long double rhs) { - - return *this; - } - - // arithmetic operators - // prefix operator - unum2 operator-() const { - return *this; - } - - unum2& operator+=(const unum2& rhs) { - return *this; - } - unum2& operator+=(double rhs) { - return *this += unum2(rhs); - } - unum2& operator-=(const unum2& rhs) { - - return *this; - } - unum2& operator-=(double rhs) { - return *this -= unum2(rhs); - } - unum2& operator*=(const unum2& rhs) { - - return *this; - } - unum2& operator*=(double rhs) { - return *this *= unum2(rhs); - } - unum2& operator/=(const unum2& rhs) { - - return *this; - } - unum2& operator/=(double rhs) { - return *this /= unum2(rhs); - } - unum2& operator++() { - return *this; - } - unum2 operator++(int) { - unum2 tmp(*this); - operator++(); - return tmp; - } - unum2& operator--() { - return *this; - } - unum2 operator--(int) { - unum2 tmp(*this); - operator--(); - return tmp; - } - - // modifiers - - /// - /// clear the content of this bfloat to zero - /// - /// void - inline constexpr void clear() noexcept { - - } - /// - /// set the number to +0 - /// - /// void - inline constexpr void setzero() noexcept { clear(); } - /// - /// set the number to +inf - /// - /// boolean to make it + or - infinity, default is -inf - /// void - inline constexpr void setinf(bool sign = true) noexcept { - - } - /// - /// set the number to a quiet NaN (+nan) or a signalling NaN (-nan, default) - /// - /// boolean to make it + or - infinity, default is -inf - /// void - inline constexpr void setnan(int NaNType = NAN_TYPE_SIGNALLING) noexcept { - } - // specific number system values of interest - inline constexpr unum2& maxpos() noexcept { - - return *this; - } - inline constexpr unum2& minpos() noexcept { - - return *this; - } - inline constexpr unum2& zero() noexcept { - - return *this; - } - inline constexpr unum2& minneg() noexcept { - - return *this; - } - inline constexpr unum2& maxneg() noexcept { - - return *this; - } - - // selectors - inline bool isneg() const { return false; } - inline bool iszero() const { return false; } - inline bool isinf() const { return false; } - inline bool isnan() const { return false; } - inline bool issnan() const { return false; } - inline bool isqnan() const { return false; } - inline bool sign() const { return false; } - inline int32_t scale() const { return false; } // 2^+-2^31 should be enough to capture empirical use cases - inline std::string get() const { return std::string("tbd"); } - - - long double to_long_double() const { - return 0.0l; - } - double to_double() const { - return 0.0; - } - float to_float() const { - return 0.0f; - } - // Maybe remove explicit - explicit operator long double() const { return to_long_double(); } - explicit operator double() const { return to_double(); } - explicit operator float() const { return to_float(); } +// Op matrix forward declaration +template +class op_matrix; + +template +class unum2 final { private: + const T& _lattice; + + // SORN + uint64_t _sorn_length; + std::bitset<1 << (sizeof(S) * 8)> _sorn; // empty - // template parameters need namfsizesize different from class template parameters (for gcc and clang) - template - friend std::ostream& operator<< (std::ostream& ostr, const unum2& r); - template - friend std::istream& operator>> (std::istream& istr, unum2& r); - - template - friend bool operator==(const unum2& lhs, const unum2& rhs); - template - friend bool operator!=(const unum2& lhs, const unum2& rhs); - template - friend bool operator< (const unum2& lhs, const unum2& rhs); - template - friend bool operator> (const unum2& lhs, const unum2& rhs); - template - friend bool operator<=(const unum2& lhs, const unum2& rhs); - template - friend bool operator>=(const unum2& lhs, const unum2& rhs); +public: + unum2(S index) : _lattice(T::instance()), _sorn_length(_lattice._N) { + if(_lattice._N > (1 << (sizeof(S) * 8))) + throw std::invalid_argument("Lattice point count overflows given storage type!"); + + // Do bitwise or with with half the lattice size. SORN index starts from 2's compliment. That is + // the first bit in SORN represents 'inf' instead of '0'. Next SORN bit represents (inf, -en) and + // so on so forth. + _sorn = std::bitset<1 << (sizeof(S) * 8)>().set((index ^ _lattice._N_half) & _lattice._MASK); + } + + unum2() : _lattice(T::instance()), _sorn_length(_lattice._N) { + _sorn.reset(); + } + + unum2(const unum2& num) : + _lattice(num._lattice), + _sorn(num._sorn), + _sorn_length(num._sorn_length) + {} + + unum2(unum2&& num) noexcept : + _lattice(num._lattice), + _sorn(std::move(num._sorn)), + _sorn_length(num._sorn_length) + {} + + unum2& operator = (const unum2& num) { + if (this != &num) { + _sorn = num._sorn; + _sorn_length = num._sorn_length; + } + return *this; + } + + unum2& operator = (unum2&& num) noexcept { + if (this != &num) { + _sorn = std::move(num._sorn); + _sorn_length = num._sorn_length; + } + return *this; + } + + friend std::ostream& operator << (std::ostream& os, const unum2& u) { + std::ostringstream oss; + + int64_t left_bound = -1; + bool bound = false; // Series of continuous 1s in the SORN bitset + bool written = false; // Has something been written to sstream? + for(S i = 0; i < u._sorn_length; i++) { + if(u._sorn[i] == 1) { + // If already bound, continue + if(bound) continue; + + // Set left bound + bound = true; + left_bound = i; + + // If something has been written, that means there were other bounds. Add Union sign. + if(written) + oss << " U "; + } else { + // Single bit bound. Can be exact or inexact. + if(bound) { + // End bound + bound = false; + written = true; + + if(left_bound == i - 1) { + // If inexact + if(left_bound & 0x01) { // Check ubit + oss << "(" << u._lattice.get_exact(u._conv_idx(left_bound - 1)) + << ", " << u._lattice.get_exact(u._conv_idx(left_bound + 1)) + << ")"; + } else { + if(left_bound == 0) { + if(u._sorn[u._lattice._N - 1] != 1) + oss << "inf"; + else written = false; + } else oss << u._lattice.get_exact(u._conv_idx(left_bound)); + } + } else { // Multiple bit bound + // Check if left_bound is inexact. If so, get previous exact. + if(left_bound & 0x01) { + left_bound--; + oss << "("; + } else oss << "["; + + int64_t right_bound; + char brace = ']'; + if((i - 1) & 0x01) { + right_bound = i; + brace = ')'; + } else right_bound = i - 1; + + oss << u._lattice.get_exact(u._conv_idx(left_bound)) << ", " + << u._lattice.get_exact(u._conv_idx(right_bound)) << brace; + } + } + } + } + + // No bounds. + if(left_bound == -1) + oss << "[EMPTY]"; + else if(bound == true && left_bound == 0) + oss << "[EVERYTHING]"; + + // Bit equal 0 code over again. + else if(bound) { + // Final bit should be 1 if there is a bound. + int64_t i = u._sorn_length - 1; + + if(left_bound == i) { + // Final bit index in SORN is always inexact + oss << "(" << u._lattice.get_exact(u._conv_idx(i - 1)); + + // If only the first SORN bit is set, that infers infinity is included. + if(u._sorn[0] == 1) + oss << ", inf]"; + else oss << ", " << u._lattice.get_exact(u._conv_idx(i + 1)) << ")"; + } else { // Multiple bit bound + if(left_bound & 0x01) { + left_bound--; + oss << "("; + } else oss << "["; + + int64_t right_bound; + char brace = ')'; + if(u._sorn[0] == 1) { + right_bound = 0; + brace = ']'; + } else right_bound = i + 1; + + oss << u._lattice.get_exact(u._conv_idx(left_bound)) << ", " + << u._lattice.get_exact(u._conv_idx(right_bound)) << brace; + } + } + + return (os << oss.str()); + } + + static unum2 empty() { + unum2 res(0); + res._sorn.reset(); + + return res; + } + + static unum2 everything() { + unum2 res(0); + res._sorn.set(); + + return res; + } + + template + static unum2 from(TT value) { + return unum2(_from_index(value)); + } + + template + static unum2 interval(TT a, TT b) { + if(a == b) + return unum2::from(a); + + S ai = unum2::_from_index(a); + S bi = unum2::_from_index(b); + + if(a < b) return unum2::_bound(ai, bi); + return unum2::_bound_inverse(ai, bi); + } + + // Addition + unum2 operator + (unum2& other) { + auto res = unum2::empty(); + for(S i = 0; i < _sorn_length; i++) { + for(S j = 0; j < _sorn_length; j++) { + if(_sorn[i] == 1 && other._sorn[j] == 1) + res._sorn |= unum2::_sumpoint(_conv_idx(i), _conv_idx(j), _lattice)._sorn; + } + } + + return res; + } + + // Multiplication + unum2 operator * (unum2& other) { + auto res = unum2::empty(); + for(S i = 0; i < _sorn_length; i++) { + for(S j = 0; j < _sorn_length; j++) { + if(_sorn[i] == 1 && other._sorn[j] == 1) + res._sorn |= unum2::_mulpoint(_conv_idx(i), _conv_idx(j), _lattice)._sorn; + } + } + + return res; + } + + // Negation + unum2 operator - () { + auto res = unum2::empty(); + for(uint64_t i = 0; i < _sorn_length; i++) { + if(_sorn[i] == 1) { + // Horizontal invert + uint64_t neg_idx = _horizontal_invert(_conv_idx(i), _lattice._MASK); + res._sorn.set(_conv_idx(neg_idx)); + } + } + + return res; + } + + // Subtraction + unum2 operator - (unum2& other) { + auto neg = -other; + return this->unum2::operator + (neg); + } + + // Invert + unum2 operator ~ () { + uint64_t msb_mask = _lattice._N >> 1; + auto res = unum2::empty(); + for(uint64_t i = 0; i < _sorn_length; i++) { + if(_sorn[i] == 1) { + // Vertical invert + uint64_t ix = _conv_idx(i); + uint64_t msb_set = ix & msb_mask; + uint64_t inverted_idx = ~ix & _lattice._MASK; + + if(!msb_set) inverted_idx &= (msb_mask - 1); + else inverted_idx |= msb_set; + inverted_idx = (inverted_idx + 1) & _lattice._MASK; + + res._sorn.set(_conv_idx(inverted_idx)); + } + } + + return res; + } + + unum2 operator / (unum2& other) { + auto inv = ~other; + return this->unum2::operator * (inv); + } + + bool operator == (unum2& other) { + return this->_sorn == other._sorn; + } + + // Raise to a power + unum2 operator ^ (double n) { + auto res = unum2::empty(); + for(S i = 0; i < _sorn_length; i++) { + if(_sorn[i] == 1) + res._sorn |= unum2::_powpoint(_conv_idx(i), n, _lattice)._sorn; + } + + return res; + } + + // Absolute + unum2 abs() { + auto res = unum2::empty(); + for(int i = 0; i < _sorn_length; i++) { + if(_sorn[i] == 1) { + uint64_t idx = _conv_idx(i); + + // Negative point + if(idx > _lattice._N_half) + // Horizontal invert + res._sorn.set(_conv_idx(_horizontal_invert(idx, _lattice._MASK))); + else res._sorn.set(i); + } + } + + return res; + } + + S _conv_idx(S idx) const { + // In Unum2 SORN, 0 starts from _N_half, instead of 0 for compatibility + // reasons. This function converts adjusted index (e.g. 15 -> 0) to + // absolute index (e.g. 0 -> 0) and vice versa. + return (idx ^ _lattice._N_half) & _lattice._MASK; + } + +private: + template + static S _from_index(TT value) { + T lattice = T(); + + if(!std::isfinite(value)) + return lattice._N >> 1; + else if(value == 0) + return 0; + else if(value == 1) + return lattice._N_quarter; + else if(value == -1) + return lattice._N_quarter * 3; + + TT absolute_value = std::abs(value); + size_t exact_size = lattice._exacts.size(); + + // Try exact values. + int64_t index = -1; + for(int64_t i = 1; i < lattice._exacts.size(); i++) { + int e = lattice._exacts[i]; + + if(absolute_value == static_cast(e)) { + // Vertical flip + index = (i << 1) + lattice._N_quarter; + break; + } + else if(absolute_value == (1 / static_cast(e))) { + index = lattice._N_quarter - (i << 1); + break; + } + + // No exacts :( + // Compare bounds. + if(absolute_value > static_cast(lattice._exacts[i - 1]) && + absolute_value < static_cast(e)) + { + // Vertical flip. + index = (i << 1) + lattice._N_quarter - 1; + break; + } else { + TT reciprocal_right = 1 / static_cast(lattice._exacts[exact_size - i - 1]); + TT reciprocal_left = 1 / static_cast(lattice._exacts[exact_size - i]); + + if(absolute_value > reciprocal_left && absolute_value < reciprocal_right) { + index = (i << 1) + 1; + break; + } + } + } + + if(index == -1) { + // Check (0, e1) or (en, inf) bounds. + if(absolute_value > 0 && absolute_value < (1 / static_cast(lattice._exacts[exact_size - 1]))) + index = 1; + else if(absolute_value > static_cast(lattice._exacts[exact_size - 1])) + index = lattice._N_half - 1; + } + + // Horizontal flip + if(value < 0) + index = ((~index & lattice._MASK) + 1) & lattice._MASK; + + return index & lattice._MASK; + } + + static unum2 _bound(S a, S b) { + // Given unum has to be points. + // and a < b + // Note: Unum 'a' is operated on, thus changes. + + if(a == b) + return a; + + auto au = unum2(a); + auto bu = unum2(b); + auto res = au._sorn; + auto criterion = res; + + // If b is infinity + if(b == au._lattice._N_half) { + int shift_count = au._lattice._N - au._sorn._Find_first(); + while(shift_count--) { + criterion <<= 1; + res |= criterion; + } + + res |= bu._sorn; + } else { + while(criterion != bu._sorn) { + criterion <<= 1; + res |= criterion; + } + } + + au._sorn = res; + return au; + } + + static unum2 _bound_inverse(S a, S b) { + // When bound a > bound b + auto res = unum2::_bound(b, a); + res._sorn = res._sorn ^ std::bitset<1 << (sizeof(S) * 8)>().set(); + // Lattice should move coutner-clockwise in this case, but including the bounded + // unums. + res._sorn.set(a ^ res._lattice._N_half).set(b ^ res._lattice._N_half); + return res; + } + + static unum2 _sumpoint(S i, S j, const T& lattice) { + op_matrix* op_mat; + + if(UNUM2_USE_OP_MATRIX) { + op_mat = &T::template op_matrix_instance(); + + if(op_mat.has(i, j, OP_MATRIX_TYPE_ADD)) + return op_mat->get(i, j, OP_MATRIX_TYPE_ADD); + } + + unum2 res; + + // i and j both represent infinity + if(i == lattice._N_half && j == lattice._N_half) { + res = unum2::everything(); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_ADD, res); + return res; + } + // i or j represent infinity + else if(i == lattice._N_half || j == lattice._N_half) { + res = unum2(lattice._N_half); // inf + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_ADD, res); + return res; + } + // i represents 0 + else if(i == 0 || j == 0) { + res = unum2(j); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_ADD, res); + return res; + } + + // is exact + bool ie = !(i & 0x01); + bool je = !(j & 0x01); + + double i_left; + double i_right; + double j_left; + double j_right; + + if(ie && je) { + res = unum2::from(lattice.exactvalue(i & lattice._MASK) + lattice.exactvalue(j & lattice._MASK)); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_ADD, res); + return res; + } else { + j_left = lattice.exactvalue((j - 1) & lattice._MASK); + j_right = lattice.exactvalue((j + 1) & lattice._MASK); + } + + if(ie) { // only i is exact + i_left = lattice.exactvalue(i & lattice._MASK); + i_right = i_left; + } + // only j is exact + else if(je) { + res = _sumpoint(j, i, lattice); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_ADD, res); + return res; + } else { + // None is exact + i_left = lattice.exactvalue((i - 1) & lattice._MASK); + i_right = lattice.exactvalue((i + 1) & lattice._MASK); + } + + S res_left_idx = unum2::_from_index(i_left + j_left); + S res_right_idx = unum2::_from_index(i_right + j_right); + + res = unum2::_bound(res_left_idx, res_right_idx); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_ADD, res); + return res; + } + + static unum2 _mulpoint(S i, S j, const T& lattice) { + op_matrix* op_mat; + + if(UNUM2_USE_OP_MATRIX) { + op_mat = &T::template op_matrix_instance(); + + if(op_mat->has(i, j, OP_MATRIX_TYPE_MUL)) { + return op_mat->get(i, j, OP_MATRIX_TYPE_MUL); + } + } + + unum2 res; + + // inf * 0 = everything + if((i == lattice._N_half && j == 0) || (i == 0 && j == lattice._N_half)) { + res = unum2::everything(); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + // inf * 1 = inf + else if((i == lattice._N_half && j == lattice._N_quarter) || (i == lattice._N_quarter && j == lattice._N_half)) { + res = unum2(lattice._N_half); // inf + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + // i represents 1 + else if(i == lattice._N_quarter) { + res = unum2(j); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + // j represents 1 + else if(j == lattice._N_quarter) { + res = unum2(i); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + // i represents 0 + else if(i == 0 || j == 0) { + res = unum2(0); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + + // is exact + bool ie = !(i & 0x01); + bool je = !(j & 0x01); + + double i_left; + double i_right; + double j_left; + double j_right; + + if(ie && je) { + res = unum2::from(lattice.exactvalue(i & lattice._MASK) * lattice.exactvalue(j & lattice._MASK)); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + else { + j_left = lattice.exactvalue((j - 1) & lattice._MASK); + j_right = lattice.exactvalue((j + 1) & lattice._MASK); + } + + if(ie) { // only i is exact + i_left = lattice.exactvalue(i & lattice._MASK); + i_right = i_left; + } + // only j is exact + else if(je) { + res = _mulpoint(j, i, lattice); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } else { + // None is exact + i_left = lattice.exactvalue((i - 1) & lattice._MASK); + i_right = lattice.exactvalue((i + 1) & lattice._MASK); + } + + // candidates + double candidates[] = { i_left * j_left, i_left * j_right, i_right * j_left, i_right * j_right }; + + // check for NaNs. That should occur only when we encounter things like inf * 0. Return everything + // in this case. + for(int i = 0; i < 4; i++) { + if(std::isnan(candidates[i])) { + res = unum2::everything(); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + } + + double res_left = *std::min_element(candidates, candidates + 4); + double res_right = *std::max_element(candidates, candidates + 4); + S res_left_idx = unum2::_from_index(res_left); + S res_right_idx = unum2::_from_index(res_right); + + if(!(res_left_idx & 0x01)) // exact + res_left_idx = (res_left_idx + 1) & lattice._MASK; + if(!(res_right_idx & 0x01)) // exact + res_right_idx = (res_right_idx - 1) & lattice._MASK; + + res = unum2::_bound(res_left_idx, res_right_idx); + if(UNUM2_USE_OP_MATRIX) op_mat->set(i, j, OP_MATRIX_TYPE_MUL, res); + return res; + } + + static unum2 _powpoint(S i, double n, const T& lattice) { + if(!(i & 0x01)) { + double value = lattice.exactvalue(i); + value = std::pow(value, n); + + // e.g. sqrt(-2) which will result complex number. + if(std::isnan(value)) + return unum2::empty(); + + return unum2::from(value); + } + + double left_bound = lattice.exactvalue((i - 1) & lattice._MASK); + double right_bound = lattice.exactvalue((i + 1) & lattice._MASK); + + left_bound = std::pow(left_bound, n); + right_bound = std::pow(right_bound, n); + if(std::isnan(left_bound) || std::isnan(right_bound)) + return unum2::empty(); + + S left_idx = unum2::_from_index(std::min(left_bound, right_bound)); + S right_idx = unum2::_from_index(std::max(left_bound, right_bound)); + + if(!(left_idx & 0x01)) // exact + left_idx = (left_idx + 1) & lattice._MASK; + if(!(right_idx & 0x01)) // exact + right_idx = (right_idx - 1) & lattice._MASK; + + return unum2::_bound(left_idx, right_idx); + } }; -////////////////////// operators -template -inline std::ostream& operator<<(std::ostream& ostr, const unum2& v) { - - return ostr; -} - -template -inline std::istream& operator>>(std::istream& istr, const unum2& v) { - istr >> v._fraction; - return istr; -} - -template -inline bool operator==(const unum2& lhs, const unum2& rhs) { return false; } -template -inline bool operator!=(const unum2& lhs, const unum2& rhs) { return !operator==(lhs, rhs); } -template -inline bool operator< (const unum2& lhs, const unum2& rhs) { return false; } -template -inline bool operator> (const unum2& lhs, const unum2& rhs) { return operator< (rhs, lhs); } -template -inline bool operator<=(const unum2& lhs, const unum2& rhs) { return !operator> (lhs, rhs); } -template -inline bool operator>=(const unum2& lhs, const unum2& rhs) { return !operator< (lhs, rhs); } - -// posit - posit binary arithmetic operators -// BINARY ADDITION -template -inline unum2 operator+(const unum2& lhs, const unum2& rhs) { - unum2 sum(lhs); - sum += rhs; - return sum; -} -// BINARY SUBTRACTION -template -inline unum2 operator-(const unum2& lhs, const unum2& rhs) { - unum2 diff(lhs); - diff -= rhs; - return diff; -} -// BINARY MULTIPLICATION -template -inline unum2 operator*(const unum2& lhs, const unum2& rhs) { - unum2 mul(lhs); - mul *= rhs; - return mul; -} -// BINARY DIVISION -template -inline unum2 operator/(const unum2& lhs, const unum2& rhs) { - unum2 ratio(lhs); - ratio /= rhs; - return ratio; -} - - -template -inline std::string components(const unum2& v) { - std::stringstream s; - if (v.iszero()) { - s << " zero b" << std::setw(esizesize) << v.fraction(); - return s.str(); - } - else if (v.isinf()) { - s << " infinite b" << std::setw(esizesize) << v.fraction(); - return s.str(); - } - s << "(" << (v.sign() ? "-" : "+") << "," << v.scale() << "," << v.fraction() << ")"; - return s.str(); -} - -/// Magnitude of a scientific notation value (equivalent to turning the sign bit off). -template -unum2 abs(const unum2& v) { - return unum2(false, v.scale(), v.fraction(), v.isZero()); -} - - -}} // namespace sw::universal +}}