From 8736f6160aeafe7a177cb6143fea80157e174e52 Mon Sep 17 00:00:00 2001 From: "Matthias P. Braendli" Date: Sun, 6 Oct 2024 19:47:19 +0200 Subject: Implement fixed-point symbols, FFT and file output --- fpm/fixed.hpp | 490 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 490 insertions(+) create mode 100644 fpm/fixed.hpp (limited to 'fpm/fixed.hpp') diff --git a/fpm/fixed.hpp b/fpm/fixed.hpp new file mode 100644 index 0000000..e2e71bf --- /dev/null +++ b/fpm/fixed.hpp @@ -0,0 +1,490 @@ +#ifndef FPM_FIXED_HPP +#define FPM_FIXED_HPP + +#include +#include +#include +#include +#include +#include + +namespace fpm +{ + +//! Fixed-point number type +//! \tparam BaseType the base integer type used to store the fixed-point number. This can be a signed or unsigned type. +//! \tparam IntermediateType the integer type used to store intermediate results during calculations. +//! \tparam FractionBits the number of bits of the BaseType used to store the fraction +//! \tparam EnableRounding enable rounding of LSB for multiplication, division, and type conversion +template +class fixed +{ + static_assert(std::is_integral::value, "BaseType must be an integral type"); + static_assert(FractionBits > 0, "FractionBits must be greater than zero"); + static_assert(FractionBits <= sizeof(BaseType) * 8 - 1, "BaseType must at least be able to contain entire fraction, with space for at least one integral bit"); + static_assert(sizeof(IntermediateType) > sizeof(BaseType), "IntermediateType must be larger than BaseType"); + static_assert(std::is_signed::value == std::is_signed::value, "IntermediateType must have same signedness as BaseType"); + + // Although this value fits in the BaseType in terms of bits, if there's only one integral bit, this value + // is incorrect (flips from positive to negative), so we must extend the size to IntermediateType. + static constexpr IntermediateType FRACTION_MULT = IntermediateType(1) << FractionBits; + + struct raw_construct_tag {}; + constexpr inline fixed(BaseType val, raw_construct_tag) noexcept : m_value(val) {} + +public: + inline fixed() noexcept = default; + + // Converts an integral number to the fixed-point type. + // Like static_cast, this truncates bits that don't fit. + template ::value>::type* = nullptr> + constexpr inline explicit fixed(T val) noexcept + : m_value(static_cast(val * FRACTION_MULT)) + {} + + // Converts an floating-point number to the fixed-point type. + // Like static_cast, this truncates bits that don't fit. + template ::value>::type* = nullptr> + constexpr inline explicit fixed(T val) noexcept + : m_value(static_cast((EnableRounding) ? + (val >= 0.0) ? (val * FRACTION_MULT + T{0.5}) : (val * FRACTION_MULT - T{0.5}) + : (val * FRACTION_MULT))) + {} + + // Constructs from another fixed-point type with possibly different underlying representation. + // Like static_cast, this truncates bits that don't fit. + template + constexpr inline explicit fixed(fixed val) noexcept + : m_value(from_fixed_point(val.raw_value()).raw_value()) + {} + + // Explicit conversion to a floating-point type + template ::value>::type* = nullptr> + constexpr inline explicit operator T() const noexcept + { + return static_cast(m_value) / FRACTION_MULT; + } + + // Explicit conversion to an integral type + template ::value>::type* = nullptr> + constexpr inline explicit operator T() const noexcept + { + return static_cast(m_value / FRACTION_MULT); + } + + // Returns the raw underlying value of this type. + // Do not use this unless you know what you're doing. + constexpr inline BaseType raw_value() const noexcept + { + return m_value; + } + + //! Constructs a fixed-point number from another fixed-point number. + //! \tparam NumFractionBits the number of bits used by the fraction in \a value. + //! \param value the integer fixed-point number + template FractionBits)>::type* = nullptr> + static constexpr inline fixed from_fixed_point(T value) noexcept + { + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + return (EnableRounding) ? fixed(static_cast( + value / (T(1) << (NumFractionBits - FractionBits)) + + (value / (T(1) << (NumFractionBits - FractionBits - 1)) % 2)), + raw_construct_tag{}) : + fixed(static_cast(value / (T(1) << (NumFractionBits - FractionBits))), + raw_construct_tag{}); + } + + template ::type* = nullptr> + static constexpr inline fixed from_fixed_point(T value) noexcept + { + return fixed(static_cast( + value * (T(1) << (FractionBits - NumFractionBits))), + raw_construct_tag{}); + } + + // Constructs a fixed-point number from its raw underlying value. + // Do not use this unless you know what you're doing. + static constexpr inline fixed from_raw_value(BaseType value) noexcept + { + return fixed(value, raw_construct_tag{}); + } + + // + // Constants + // + static constexpr fixed e() { return from_fixed_point<61>(6267931151224907085ll); } + static constexpr fixed pi() { return from_fixed_point<61>(7244019458077122842ll); } + static constexpr fixed half_pi() { return from_fixed_point<62>(7244019458077122842ll); } + static constexpr fixed two_pi() { return from_fixed_point<60>(7244019458077122842ll); } + + // + // Arithmetic member operators + // + + constexpr inline fixed operator-() const noexcept + { + return fixed::from_raw_value(-m_value); + } + + inline fixed& operator+=(const fixed& y) noexcept + { + m_value += y.m_value; + return *this; + } + + template ::value>::type* = nullptr> + inline fixed& operator+=(I y) noexcept + { + m_value += y * FRACTION_MULT; + return *this; + } + + inline fixed& operator-=(const fixed& y) noexcept + { + m_value -= y.m_value; + return *this; + } + + template ::value>::type* = nullptr> + inline fixed& operator-=(I y) noexcept + { + m_value -= y * FRACTION_MULT; + return *this; + } + + inline fixed& operator*=(const fixed& y) noexcept + { + if (EnableRounding){ + // Normal fixed-point multiplication is: x * y / 2**FractionBits. + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + auto value = (static_cast(m_value) * y.m_value) / (FRACTION_MULT / 2); + m_value = static_cast((value / 2) + (value % 2)); + } else { + auto value = (static_cast(m_value) * y.m_value) / FRACTION_MULT; + m_value = static_cast(value); + } + return *this; + } + + template ::value>::type* = nullptr> + inline fixed& operator*=(I y) noexcept + { + m_value *= y; + return *this; + } + + inline fixed& operator/=(const fixed& y) noexcept + { + assert(y.m_value != 0); + if (EnableRounding){ + // Normal fixed-point division is: x * 2**FractionBits / y. + // To correctly round the last bit in the result, we need one more bit of information. + // We do this by multiplying by two before dividing and adding the LSB to the real result. + auto value = (static_cast(m_value) * FRACTION_MULT * 2) / y.m_value; + m_value = static_cast((value / 2) + (value % 2)); + } else { + auto value = (static_cast(m_value) * FRACTION_MULT) / y.m_value; + m_value = static_cast(value); + } + return *this; + } + + template ::value>::type* = nullptr> + inline fixed& operator/=(I y) noexcept + { + m_value /= y; + return *this; + } + +private: + BaseType m_value; +}; + +// +// Convenience typedefs +// + +using fixed_16_16 = fixed; +using fixed_24_8 = fixed; +using fixed_8_24 = fixed; + +// +// Addition +// + +template +constexpr inline fixed operator+(const fixed& x, const fixed& y) noexcept +{ + return fixed(x) += y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator+(const fixed& x, T y) noexcept +{ + return fixed(x) += y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator+(T x, const fixed& y) noexcept +{ + return fixed(y) += x; +} + +// +// Subtraction +// + +template +constexpr inline fixed operator-(const fixed& x, const fixed& y) noexcept +{ + return fixed(x) -= y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator-(const fixed& x, T y) noexcept +{ + return fixed(x) -= y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator-(T x, const fixed& y) noexcept +{ + return fixed(x) -= y; +} + +// +// Multiplication +// + +template +constexpr inline fixed operator*(const fixed& x, const fixed& y) noexcept +{ + return fixed(x) *= y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator*(const fixed& x, T y) noexcept +{ + return fixed(x) *= y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator*(T x, const fixed& y) noexcept +{ + return fixed(y) *= x; +} + +// +// Division +// + +template +constexpr inline fixed operator/(const fixed& x, const fixed& y) noexcept +{ + return fixed(x) /= y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator/(const fixed& x, T y) noexcept +{ + return fixed(x) /= y; +} + +template ::value>::type* = nullptr> +constexpr inline fixed operator/(T x, const fixed& y) noexcept +{ + return fixed(x) /= y; +} + +// +// Comparison operators +// + +template +constexpr inline bool operator==(const fixed& x, const fixed& y) noexcept +{ + return x.raw_value() == y.raw_value(); +} + +template +constexpr inline bool operator!=(const fixed& x, const fixed& y) noexcept +{ + return x.raw_value() != y.raw_value(); +} + +template +constexpr inline bool operator<(const fixed& x, const fixed& y) noexcept +{ + return x.raw_value() < y.raw_value(); +} + +template +constexpr inline bool operator>(const fixed& x, const fixed& y) noexcept +{ + return x.raw_value() > y.raw_value(); +} + +template +constexpr inline bool operator<=(const fixed& x, const fixed& y) noexcept +{ + return x.raw_value() <= y.raw_value(); +} + +template +constexpr inline bool operator>=(const fixed& x, const fixed& y) noexcept +{ + return x.raw_value() >= y.raw_value(); +} + +namespace detail +{ +// Number of base-10 digits required to fully represent a number of bits +static constexpr int max_digits10(int bits) +{ + // 8.24 fixed-point equivalent of (int)ceil(bits * std::log10(2)); + using T = long long; + return static_cast((T{bits} * 5050445 + (T{1} << 24) - 1) >> 24); +} + +// Number of base-10 digits that can be fully represented by a number of bits +static constexpr int digits10(int bits) +{ + // 8.24 fixed-point equivalent of (int)(bits * std::log10(2)); + using T = long long; + return static_cast((T{bits} * 5050445) >> 24); +} + +} // namespace detail +} // namespace fpm + +// Specializations for customization points +namespace std +{ + +template +struct hash> +{ + using argument_type = fpm::fixed; + using result_type = std::size_t; + + result_type operator()(argument_type arg) const noexcept(noexcept(std::declval>()(arg.raw_value()))) { + return m_hash(arg.raw_value()); + } + +private: + std::hash m_hash; +}; + +template +struct numeric_limits> +{ + static constexpr bool is_specialized = true; + static constexpr bool is_signed = std::numeric_limits::is_signed; + static constexpr bool is_integer = false; + static constexpr bool is_exact = true; + static constexpr bool has_infinity = false; + static constexpr bool has_quiet_NaN = false; + static constexpr bool has_signaling_NaN = false; + static constexpr std::float_denorm_style has_denorm = std::denorm_absent; + static constexpr bool has_denorm_loss = false; + static constexpr std::float_round_style round_style = std::round_to_nearest; + static constexpr bool is_iec_559 = false; + static constexpr bool is_bounded = true; + static constexpr bool is_modulo = std::numeric_limits::is_modulo; + static constexpr int digits = std::numeric_limits::digits; + + // Any number with `digits10` significant base-10 digits (that fits in + // the range of the type) is guaranteed to be convertible from text and + // back without change. Worst case, this is 0.000...001, so we can only + // guarantee this case. Nothing more. + static constexpr int digits10 = 1; + + // This is equal to max_digits10 for the integer and fractional part together. + static constexpr int max_digits10 = + fpm::detail::max_digits10(std::numeric_limits::digits - F) + fpm::detail::max_digits10(F); + + static constexpr int radix = 2; + static constexpr int min_exponent = 1 - F; + static constexpr int min_exponent10 = -fpm::detail::digits10(F); + static constexpr int max_exponent = std::numeric_limits::digits - F; + static constexpr int max_exponent10 = fpm::detail::digits10(std::numeric_limits::digits - F); + static constexpr bool traps = true; + static constexpr bool tinyness_before = false; + + static constexpr fpm::fixed lowest() noexcept { + return fpm::fixed::from_raw_value(std::numeric_limits::lowest()); + }; + + static constexpr fpm::fixed min() noexcept { + return lowest(); + } + + static constexpr fpm::fixed max() noexcept { + return fpm::fixed::from_raw_value(std::numeric_limits::max()); + }; + + static constexpr fpm::fixed epsilon() noexcept { + return fpm::fixed::from_raw_value(1); + }; + + static constexpr fpm::fixed round_error() noexcept { + return fpm::fixed(1) / 2; + }; + + static constexpr fpm::fixed denorm_min() noexcept { + return min(); + } +}; + +template +constexpr bool numeric_limits>::is_specialized; +template +constexpr bool numeric_limits>::is_signed; +template +constexpr bool numeric_limits>::is_integer; +template +constexpr bool numeric_limits>::is_exact; +template +constexpr bool numeric_limits>::has_infinity; +template +constexpr bool numeric_limits>::has_quiet_NaN; +template +constexpr bool numeric_limits>::has_signaling_NaN; +template +constexpr std::float_denorm_style numeric_limits>::has_denorm; +template +constexpr bool numeric_limits>::has_denorm_loss; +template +constexpr std::float_round_style numeric_limits>::round_style; +template +constexpr bool numeric_limits>::is_iec_559; +template +constexpr bool numeric_limits>::is_bounded; +template +constexpr bool numeric_limits>::is_modulo; +template +constexpr int numeric_limits>::digits; +template +constexpr int numeric_limits>::digits10; +template +constexpr int numeric_limits>::max_digits10; +template +constexpr int numeric_limits>::radix; +template +constexpr int numeric_limits>::min_exponent; +template +constexpr int numeric_limits>::min_exponent10; +template +constexpr int numeric_limits>::max_exponent; +template +constexpr int numeric_limits>::max_exponent10; +template +constexpr bool numeric_limits>::traps; +template +constexpr bool numeric_limits>::tinyness_before; + +} + +#endif -- cgit v1.2.3