Add floating point helpers, and the fast ilog2

benchmarks/BM_mathfuns.cpp

@@ -22,6 +22,7 @@ public:
         std::uniform_real_distribution<float> dist { 0.1f, 1.0f };
         source = std::vector<float>(state.range(0));
         result = std::vector<float>(state.range(0));
+        intResult = std::vector<int>(state.range(0));
         std::generate(source.begin(), source.end(), [&]() { return dist(gen); });
     }
 
@@ -31,6 +32,7 @@ public:
 
     std::vector<float> source;
     std::vector<float> result;
+    std::vector<int> 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<int>(
+                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();

src/sfizz/MathHelpers.h

@@ -163,3 +163,124 @@ template <class Type>
 constexpr Type sqrtTwo { static_cast<Type>(1.414213562373095048801688724209698078569671875376948073176) };
 template <class Type>
 constexpr Type sqrtTwoInv { static_cast<Type>(0.707106781186547524400844362104849039284835937688474036588) };
+
+/**
+   @brief A fraction which is parameterized by integer type
+ */
+template <class I>
+struct Fraction {
+    typedef I value_type;
+
+    operator double() const noexcept;
+    operator float() const noexcept;
+
+    I num;
+    I den;
+};
+
+template <class I>
+inline Fraction<I>::operator double() const noexcept
+{
+    return static_cast<double>(num) / static_cast<double>(den);
+}
+
+template <class I>
+inline Fraction<I>::operator float() const noexcept
+{
+    return static_cast<float>(num) / static_cast<float>(den);
+}
+
+/**
+   @brief Characteristics of IEEE754 floating point representations
+ */
+template <class F>
+struct FP_traits;
+
+template <> struct FP_traits<double>
+{
+    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<float>
+{
+    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 <class F>
+inline bool fp_sign(F x)
+{
+    typedef FP_traits<F> 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 <class F>
+inline int fp_exponent(F x)
+{
+    typedef FP_traits<F> 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 <class F>
+inline Fraction<uint64_t> fp_mantissa(F x)
+{
+    typedef FP_traits<F> T;
+    union { F real; typename T::same_size_int integer; } u;
+    u.real = x;
+    Fraction<uint64_t> 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 <class F>
+inline F fp_from_parts(bool sgn, int ex, uint64_t mant)
+{
+    typedef FP_traits<F> T;
+    typedef typename T::same_size_int I;
+    union { F real; I integer; } u;
+    u.integer = mant |
+        (static_cast<I>(ex - T::e_offset) << T::m_bits) |
+        (static_cast<I>(sgn) << (T::e_bits + T::m_bits));
+    return u.real;
+}

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})

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 <cmath>
+
+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<int>(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<int>(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<uint64_t> mant = fp_mantissa(f);
+
+        REQUIRE(fp_from_parts<float>(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<uint64_t> mant = fp_mantissa(f);
+
+        REQUIRE(fp_from_parts<double>(sgn, ex, mant.num) == f);
+    }
+}