From db321e6a00d48325657cf1e6e2b34eacf355ee0d Mon Sep 17 00:00:00 2001 From: Jean Pierre Cimalando Date: Thu, 13 Feb 2020 07:37:06 +0100 Subject: [PATCH] Add floating point helpers, and the fast ilog2 --- benchmarks/BM_mathfuns.cpp | 27 +++++++++ src/sfizz/MathHelpers.h | 121 +++++++++++++++++++++++++++++++++++++ tests/CMakeLists.txt | 1 + tests/FloatHelpersT.cpp | 53 ++++++++++++++++ 4 files changed, 202 insertions(+) create mode 100644 tests/FloatHelpersT.cpp diff --git a/benchmarks/BM_mathfuns.cpp b/benchmarks/BM_mathfuns.cpp index 1c64b553..245a283f 100644 --- a/benchmarks/BM_mathfuns.cpp +++ b/benchmarks/BM_mathfuns.cpp @@ -22,6 +22,7 @@ public: std::uniform_real_distribution dist { 0.1f, 1.0f }; source = std::vector(state.range(0)); result = std::vector(state.range(0)); + intResult = std::vector(state.range(0)); std::generate(source.begin(), source.end(), [&]() { return dist(gen); }); } @@ -31,6 +32,7 @@ public: std::vector source; std::vector result; + std::vector intResult; }; BENCHMARK_DEFINE_F(MyFixture, Dummy) @@ -133,6 +135,29 @@ BENCHMARK_DEFINE_F(MyFixture, SIMDCos) } } +BENCHMARK_DEFINE_F(MyFixture, ScalarLibmFloorLog2) +(benchmark::State& state) +{ + for (auto _ : state) { + for (size_t i = 0, n = source.size(); i < n; ++i) { + intResult[i] = static_cast( + std::floor(std::log2(std::fabs(source[i])))); + } + benchmark::DoNotOptimize(intResult); + } +} + +BENCHMARK_DEFINE_F(MyFixture, ScalarFastFloorLog2) +(benchmark::State& state) +{ + for (auto _ : state) { + for (size_t i = 0, n = source.size(); i < n; ++i) { + intResult[i] = fp_exponent(source[i]); + } + benchmark::DoNotOptimize(intResult); + } +} + BENCHMARK_REGISTER_F(MyFixture, Dummy)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); BENCHMARK_REGISTER_F(MyFixture, ScalarExp)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); BENCHMARK_REGISTER_F(MyFixture, SIMDExp)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); @@ -144,5 +169,7 @@ BENCHMARK_REGISTER_F(MyFixture, ScalarSin)->RangeMultiplier(4)->Range(1 << 6, 1 BENCHMARK_REGISTER_F(MyFixture, SIMDSin)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); BENCHMARK_REGISTER_F(MyFixture, ScalarCos)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); BENCHMARK_REGISTER_F(MyFixture, SIMDCos)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); +BENCHMARK_REGISTER_F(MyFixture, ScalarLibmFloorLog2)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); +BENCHMARK_REGISTER_F(MyFixture, ScalarFastFloorLog2)->RangeMultiplier(4)->Range(1 << 6, 1 << 10); BENCHMARK_MAIN(); diff --git a/src/sfizz/MathHelpers.h b/src/sfizz/MathHelpers.h index 07229ca8..17e7bb8e 100644 --- a/src/sfizz/MathHelpers.h +++ b/src/sfizz/MathHelpers.h @@ -163,3 +163,124 @@ template constexpr Type sqrtTwo { static_cast(1.414213562373095048801688724209698078569671875376948073176) }; template constexpr Type sqrtTwoInv { static_cast(0.707106781186547524400844362104849039284835937688474036588) }; + +/** + @brief A fraction which is parameterized by integer type + */ +template +struct Fraction { + typedef I value_type; + + operator double() const noexcept; + operator float() const noexcept; + + I num; + I den; +}; + +template +inline Fraction::operator double() const noexcept +{ + return static_cast(num) / static_cast(den); +} + +template +inline Fraction::operator float() const noexcept +{ + return static_cast(num) / static_cast(den); +} + +/** + @brief Characteristics of IEEE754 floating point representations + */ +template +struct FP_traits; + +template <> struct FP_traits +{ + typedef double type; + typedef uint64_t same_size_int; + static_assert(sizeof(type) == sizeof(same_size_int), + "Unexpected size of floating point type"); + static constexpr int e_bits = 11; + static constexpr int m_bits = 52; + static constexpr int e_offset = -1023; +}; + +template <> struct FP_traits +{ + typedef float type; + typedef uint32_t same_size_int; + static_assert(sizeof(type) == sizeof(same_size_int), + "Unexpected size of floating point type"); + static constexpr int e_bits = 8; + static constexpr int m_bits = 23; + static constexpr int e_offset = -127; +}; + +/** + @brief Get the sign part of a IEEE754 floating point number. + + The number is reconstructed as `(-1^sign)*(1+mantissa)*(2^exponent)`. + See also `fp_exponent` and `fp_mantissa`. + */ +template +inline bool fp_sign(F x) +{ + typedef FP_traits T; + union { F real; typename T::same_size_int integer; } u; + u.real = x; + return ((u.integer >> (T::e_bits + T::m_bits)) & 1) != 0; +} + +/** + @brief Get the exponent part of a IEEE754 floating point number. + + The number is reconstructed as `(-1^sign)*(1+mantissa)*(2^exponent)`. + See also `fp_sign` and `fp_mantissa`. + + It is a faster way of computing `floor(log2(abs(x)))`. + */ +template +inline int fp_exponent(F x) +{ + typedef FP_traits T; + union { F real; typename T::same_size_int integer; } u; + u.real = x; + int ex = (u.integer >> T::m_bits) & ((1u << T::e_bits) - 1); + return ex + T::e_offset; +} + +/** + @brief Get the mantissa part of a IEEE754 floating point number. + The number is reconstructed as `(-1^sign)*(1+mantissa)*(2^exponent)`. + See also `fp_sign` and `fp_exponent`. + */ +template +inline Fraction fp_mantissa(F x) +{ + typedef FP_traits T; + union { F real; typename T::same_size_int integer; } u; + u.real = x; + Fraction f; + f.den = uint64_t{1} << T::m_bits; + f.num = u.integer & (f.den - 1); + return f; +} + +/** + @brief Reconstruct a IEEE754 floating point number from its parts. + + The parts must be in their range of validity. + */ +template +inline F fp_from_parts(bool sgn, int ex, uint64_t mant) +{ + typedef FP_traits T; + typedef typename T::same_size_int I; + union { F real; I integer; } u; + u.integer = mant | + (static_cast(ex - T::e_offset) << T::m_bits) | + (static_cast(sgn) << (T::e_bits + T::m_bits)); + return u.real; +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2bb008a5..72149427 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -24,6 +24,7 @@ set(SFIZZ_TEST_SOURCES MainT.cpp SynthT.cpp RegionTriggersT.cpp + FloatHelpersT.cpp ) add_executable(sfizz_tests ${SFIZZ_TEST_SOURCES}) diff --git a/tests/FloatHelpersT.cpp b/tests/FloatHelpersT.cpp new file mode 100644 index 00000000..58b92acc --- /dev/null +++ b/tests/FloatHelpersT.cpp @@ -0,0 +1,53 @@ +// SPDX-License-Identifier: BSD-2-Clause + +// This code is part of the sfizz library and is licensed under a BSD 2-clause +// license. You should have receive a LICENSE.md file along with the code. +// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz + +#include "catch2/catch.hpp" +#include "sfizz/MathHelpers.h" +#include + +TEST_CASE("[FloatMath] Fast ilog2 (float)") +{ + for (float x = -100.0f; x < +100.0f; x += 0.01f) { + int ex1 = fp_exponent(x); + int ex2 = static_cast(std::floor(std::log2(std::fabs(x)))); + REQUIRE(ex1 == ex2); + } +} + +TEST_CASE("[FloatMath] Fast ilog2 (double)") +{ + for (double x = -100.0; x < +100.0; x += 0.01) { + int ex1 = fp_exponent(x); + int ex2 = static_cast(std::floor(std::log2(std::fabs(x)))); + REQUIRE(ex1 == ex2); + } +} + +TEST_CASE("[FloatMath] Break apart and reconstruct (float)") +{ + for (int p = 0; p < 128; ++p) { + float f = 440.0 * std::pow(2.0, (p - 69.0) / 12.0); + + bool sgn = fp_sign(f); + int ex = fp_exponent(f); + Fraction mant = fp_mantissa(f); + + REQUIRE(fp_from_parts(sgn, ex, mant.num) == f); + } +} + +TEST_CASE("[FloatMath] Break apart and reconstruct (double)") +{ + for (int p = 0; p < 128; ++p) { + double f = 440.0 * std::pow(2.0, (p - 69.0) / 12.0); + + bool sgn = fp_sign(f); + int ex = fp_exponent(f); + Fraction mant = fp_mantissa(f); + + REQUIRE(fp_from_parts(sgn, ex, mant.num) == f); + } +}