Added a cumsum helper

benchmarks/BM_cumsum.cpp

@@ -0,0 +1,84 @@
+// Copyright (c) 2019, Paul Ferrand
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice, this
+//    list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+#include <benchmark/benchmark.h>
+#include <random>
+#include <numeric>
+#include <vector>
+#include <cmath>
+#include <iostream>
+#include "../sfizz/SIMDHelpers.h"
+
+class CumsumArray : public benchmark::Fixture {
+public:
+  void SetUp(const ::benchmark::State& state) {
+    std::random_device rd { };
+    std::mt19937 gen { rd() };
+    std::uniform_real_distribution<float> dist { 0, 1 };
+    input = std::vector<float>(state.range(0));
+    output = std::vector<float>(state.range(0));
+    std::generate(input.begin(), input.end(), [&]() { return dist(gen); });
+  }
+
+  void TearDown(const ::benchmark::State& state [[maybe_unused]]) {
+
+  }
+
+  std::vector<float> input;
+  std::vector<float> output;
+};
+
+
+BENCHMARK_DEFINE_F(CumsumArray, Scalar)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        cumsum<float, false>(input, absl::MakeSpan(output));
+    }
+}
+
+BENCHMARK_DEFINE_F(CumsumArray, SIMD)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        cumsum<float, true>(input, absl::MakeSpan(output));
+    }
+}
+
+BENCHMARK_DEFINE_F(CumsumArray, Scalar_Unaligned)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        cumsum<float, false>(absl::MakeSpan(input).subspan(1), absl::MakeSpan(output).subspan(1));
+    }
+}
+
+BENCHMARK_DEFINE_F(CumsumArray, SIMD_Unaligned)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        cumsum<float, true>(absl::MakeSpan(input).subspan(1), absl::MakeSpan(output).subspan(1));
+    }
+}
+
+BENCHMARK_REGISTER_F(CumsumArray, Scalar)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
+BENCHMARK_REGISTER_F(CumsumArray, SIMD)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
+BENCHMARK_REGISTER_F(CumsumArray, Scalar_Unaligned)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
+BENCHMARK_REGISTER_F(CumsumArray, SIMD_Unaligned)->RangeMultiplier(4)->Range(1 << 2, 1 << 12);
+BENCHMARK_MAIN();

benchmarks/CMakeLists.txt

@@ -75,6 +75,9 @@ target_link_libraries(bm_mean benchmark absl::span absl::algorithm)
 add_executable(bm_meanSquared BM_meanSquared.cpp ${SFIZZ_SIMD_SOURCES})
 target_link_libraries(bm_meanSquared benchmark absl::span absl::algorithm)
 
+add_executable(bm_cumsum BM_cumsum.cpp ${SFIZZ_SIMD_SOURCES})
+target_link_libraries(bm_cumsum benchmark absl::span absl::algorithm)
+
 add_custom_target(sfizz_benchmarks)
 add_dependencies(sfizz_benchmarks 
 	bm_opf_high_vs_low 
@@ -83,6 +86,7 @@ add_dependencies(sfizz_benchmarks
 	bm_mean
 	bm_meanSquared
 	bm_fill
+	bm_cumsum
 	bm_mathfuns
 	bm_gain
 	bm_looping

sfizz/Config.h

@@ -61,6 +61,7 @@ namespace SIMDConfig {
     constexpr bool multiplyAdd { false };
     constexpr bool copy { false };
     constexpr bool pan { true };
+    constexpr bool cumsum { true };
     constexpr bool mean { false };
     constexpr bool meanSquared { false };
 }

sfizz/SIMDDummy.cpp

@@ -143,3 +143,9 @@ float meanSquared<float, true>(absl::Span<const float> vector) noexcept
 {
     return meanSquared<float, false>(vector);
 }
+
+template <>
+void cumsum<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept
+{
+    cumsum<float, false>(input, output);
+}

sfizz/SIMDHelpers.h

@@ -449,3 +449,29 @@ T meanSquared(absl::Span<const T> vector) noexcept
 
 template <>
 float meanSquared<float, true>(absl::Span<const float> vector) noexcept;
+
+template <class T>
+inline void snippetCumsum(const T*& input, T*& output)
+{
+    *output = *(output - 1) + *input++;
+    output++;
+}
+
+template <class T, bool SIMD = SIMDConfig::cumsum>
+void cumsum(absl::Span<const float> input, absl::Span<float> output) noexcept
+{
+    ASSERT(output.size() >= input.size());
+    if (input.size() == 0)
+        return;
+
+    auto out = output.data();
+    auto in = input.data();
+    const auto sentinel = in + std::min(input.size(), output.size());
+
+    *out++ = *in++;
+    while (in < sentinel)
+        snippetCumsum(in, out);
+}
+
+template <>
+void cumsum<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept;

sfizz/SIMDSSE.cpp

@@ -613,7 +613,7 @@ float mean<float, true>(absl::Span<const float> vector) noexcept
     auto* sentinel = vector.end();
     const auto* lastAligned = prevAligned(sentinel);
 
-    while (unaligned(value))
+    while (unaligned(value) && value < lastAligned)
         result += *value++;
     
     auto mmSums = _mm_setzero_ps();
@@ -645,7 +645,7 @@ float meanSquared<float, true>(absl::Span<const float> vector) noexcept
     auto* sentinel = vector.end();
     const auto* lastAligned = prevAligned(sentinel);
 
-    while (unaligned(value)) {
+    while (unaligned(value) && value < lastAligned) {
         result += (*value) * (*value);
         value++;
     }
@@ -670,3 +670,35 @@ float meanSquared<float, true>(absl::Span<const float> vector) noexcept
 
     return result / static_cast<float>(vector.size());
 }
+
+template <>
+void cumsum<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept
+{
+    ASSERT(output.size() >= input.size());
+    if (input.size() == 0)
+        return;
+
+    auto out = output.data();
+    auto in = input.data();
+    const auto sentinel = in + std::min(input.size(), output.size());
+    const auto lastAligned = prevAligned(sentinel);
+
+    *out++ = *in++;
+    while (unaligned(in, out) && in < lastAligned)
+        snippetCumsum(in, out);
+
+    auto mmOutput = _mm_set_ps1(*(out - 1));
+    while (in < lastAligned) {
+        auto mmOffset = _mm_load_ps(in);
+        mmOffset = _mm_add_ps(mmOffset, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(mmOffset), 4)));
+        mmOffset = _mm_add_ps(mmOffset, _mm_shuffle_ps(_mm_setzero_ps(), mmOffset, 0x40));
+        mmOutput = _mm_add_ps(mmOutput, mmOffset);
+        _mm_store_ps(out, mmOutput);
+        mmOutput = _mm_shuffle_ps(mmOutput, mmOutput, _MM_SHUFFLE(3, 3, 3, 3));
+        in += TypeAlignment;
+        out += TypeAlignment;
+    }
+
+    while (in < sentinel)
+        snippetCumsum(in, out);
+}

tests/SIMDHelpersT.cpp

@@ -712,3 +712,23 @@ TEST_CASE("[Helpers] Mean Squared (SIMD vs scalar)")
     absl::c_iota(input, 0.0);
     REQUIRE(meanSquared<float, false>(input) == meanSquared<float, true>(input));
 }
+
+TEST_CASE("[Helpers] Cumulative sum ")
+{
+    std::array<float, 6> input { 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f }; // 1.1 2.3 3.6 5.0 6.5 8.1
+    std::array<float, 6> output; 
+    std::array<float, 6> expected { 1.1f, 2.3f, 3.6f, 5.0f, 6.5f, 8.1f };
+    cumsum<float, false>(input, absl::MakeSpan(output));
+    REQUIRE(approxEqual<float>(output, expected));
+}
+
+TEST_CASE("[Helpers] Cumulative sum (SIMD vs Scalar)")
+{
+    std::vector<float> input(bigBufferSize);
+    std::vector<float> outputScalar(bigBufferSize);
+    std::vector<float> outputSIMD(bigBufferSize);
+    linearRamp<float>(absl::MakeSpan(input), 0.0, 0.1);
+    cumsum<float, false>(input, absl::MakeSpan(outputScalar));
+    cumsum<float, true>(input, absl::MakeSpan(outputSIMD));
+    REQUIRE(approxEqual<float>(outputScalar, outputSIMD));
+}