[C++][Parquet] Scan DELTA_BINARY_PACKED deltas a vector at a time

The value-at-a-time prefix sum becomes a helper that scans whole registers with
a log-step inclusive scan and finishes the remainder one value at a time. Adding
the frame of reference before the scan turns its running multiple into a term the
scan produces rather than a multiply per lane, and the running total is carried
between registers in a vector register rather than through a general-purpose one.

The xsimd include and the vector loop are guarded on ARROW_HAVE_NEON or
ARROW_HAVE_SSE4_2, so a build with neither compiles the scalar loop alone. The
vector loop is compiled only where a register holds four or more values, which at
the 128-bit baseline vectorizes 32-bit values and leaves 64-bit ones on the
value-at-a-time loop.

Also adds a decode benchmark on non-decreasing values, the shape this encoding is
usually chosen for.

Arithmetic is unchanged: every term stays in the unsigned type, so the wrapping
the format specifies is preserved and no decoded value changes.

Author: sfc-gh-pgaur

cpp/src/parquet/decoder.cc

@@ -32,6 +32,10 @@
 #include <utility>
 #include <vector>
 
+#if defined(ARROW_HAVE_NEON) || defined(ARROW_HAVE_SSE4_2)
+#  include <xsimd/xsimd.hpp>
+#endif
+
 #include "arrow/array.h"
 #include "arrow/array/builder_binary.h"
 #include "arrow/array/builder_dict.h"
@@ -1437,6 +1441,71 @@ class DictByteArrayDecoderImpl : public DictDecoderImpl<ByteArrayType> {
 // ----------------------------------------------------------------------
 // DELTA_BINARY_PACKED decoder
 
+namespace {
+
+#if defined(ARROW_HAVE_NEON) || defined(ARROW_HAVE_SSE4_2)
+// One step of an inclusive scan, recursed at compile time over powers of two. Each step
+// adds the vector to a copy of itself slid up by kShift lanes, zero-filling the lanes
+// it vacates, so after the last step lane k holds the sum of lanes 0 through k.
+template <std::size_t kShift, typename Batch>
+Batch InclusiveScanSteps(Batch v) {
+  if constexpr (kShift < Batch::size) {
+    v += xsimd::slide_left<kShift * sizeof(typename Batch::value_type)>(v);
+    return InclusiveScanSteps<kShift * 2>(v);
+  } else {
+    return v;
+  }
+}
+#endif
+
+// Turns a run of deltas into the values they encode, in place: on return element k holds
+// `last + (k + 1) * min_delta + sum of deltas 0..k`, and the return value is the last
+// element. Every term is unsigned, so the wrapping matches the loop this replaces.
+template <typename T>
+std::make_unsigned_t<T> PrefixSumDeltas(T* values, int num_values,
+                                        std::make_unsigned_t<T> min_delta,
+                                        std::make_unsigned_t<T> last) {
+  using UT = std::make_unsigned_t<T>;
+  int i = 0;
+
+#if defined(ARROW_HAVE_NEON) || defined(ARROW_HAVE_SSE4_2)
+  using Batch = xsimd::batch<UT>;
+  constexpr int kLanes = static_cast<int>(Batch::size);
+  // At two lanes the scan loses to the additions it replaces, so it is compiled only
+  // where a register holds four or more values; narrower ones use the loop below.
+  if constexpr (kLanes >= 4) {
+    // Broadcast pattern for the last lane, which carries the running value into the
+    // next vector without a round trip through a general-purpose register.
+    struct LastLane {
+      static constexpr unsigned get(unsigned /*index*/, unsigned size) {
+        return size - 1;
+      }
+    };
+    const auto last_lane =
+        xsimd::make_batch_constant<UT, LastLane, xsimd::default_arch>();
+    const Batch min_delta_v(min_delta);
+    Batch carry(last);
+    for (; i + kLanes <= num_values; i += kLanes) {
+      // Adding the frame before the scan turns its running multiple into a term the
+      // scan produces, rather than a multiply per lane.
+      Batch v = xsimd::bitwise_cast<UT>(xsimd::batch<T>::load_unaligned(values + i));
+      v = InclusiveScanSteps<1>(v + min_delta_v) + carry;
+      xsimd::bitwise_cast<T>(v).store_unaligned(values + i);
+      carry = xsimd::swizzle(v, last_lane);
+    }
+    last = carry.get(0);
+  }
+#endif
+
+  for (; i < num_values; ++i) {
+    last += min_delta + static_cast<UT>(values[i]);
+    values[i] = static_cast<T>(last);
+  }
+  return last;
+}
+
+}  // namespace
+
 template <typename DType>
 class DeltaBitPackDecoder : public TypedDecoderImpl<DType> {
  public:
@@ -1694,17 +1763,9 @@ class DeltaBitPackDecoder : public TypedDecoderImpl<DType> {
             values_decode) {
           ParquetException::EofException();
         }
-        // Held in locals because a store to `buffer` may alias these members, which
-        // would force a reload of each per value.
-        UT last = static_cast<UT>(last_value_);
-        const UT min_delta = static_cast<UT>(min_delta_);
-        for (int j = 0; j < values_decode; ++j) {
-          // Addition between min_delta, packed int and last_value should be treated as
-          // unsigned addition. Overflow is as expected.
-          last += min_delta + static_cast<UT>(buffer[i + j]);
-          buffer[i + j] = last;
-        }
-        last_value_ = static_cast<T>(last);
+        last_value_ = static_cast<T>(PrefixSumDeltas(buffer + i, values_decode,
+                                                     static_cast<UT>(min_delta_),
+                                                     static_cast<UT>(last_value_)));
       }
       // The last miniblock of a coalesced run becomes the current one, fully drained.
       mini_block_idx_ += mini_blocks_coalesced;

cpp/src/parquet/encoding_benchmark.cc

@@ -675,6 +675,22 @@ static auto MakeDeltaBitPackingInputNarrow(size_t length) {
   return numbers;
 }
 
+// Non-decreasing values, the shape this encoding is usually chosen for. The deltas come
+// from the same 1000-wide range as Narrow, so the two inputs differ in the order of the
+// values and not in the width a delta is packed into.
+template <typename DType>
+static auto MakeDeltaBitPackingInputNarrowSorted(size_t length) {
+  using T = typename DType::c_type;
+  auto numbers = std::vector<T>(length);
+  ::arrow::randint<T, T>(length, 0, 1000, &numbers);
+  T value = 0;
+  for (auto& number : numbers) {
+    value = static_cast<T>(value + number);
+    number = value;
+  }
+  return numbers;
+}
+
 template <typename DType>
 static auto MakeDeltaBitPackingInputWide(size_t length) {
   using T = typename DType::c_type;
@@ -713,6 +729,16 @@ static void BM_DeltaBitPackingEncode_Int64_Narrow(benchmark::State& state) {
   BM_DeltaBitPackingEncode<Int64Type>(state, MakeDeltaBitPackingInputNarrow<Int64Type>);
 }
 
+static void BM_DeltaBitPackingEncode_Int32_NarrowSorted(benchmark::State& state) {
+  BM_DeltaBitPackingEncode<Int32Type>(state,
+                                      MakeDeltaBitPackingInputNarrowSorted<Int32Type>);
+}
+
+static void BM_DeltaBitPackingEncode_Int64_NarrowSorted(benchmark::State& state) {
+  BM_DeltaBitPackingEncode<Int64Type>(state,
+                                      MakeDeltaBitPackingInputNarrowSorted<Int64Type>);
+}
+
 static void BM_DeltaBitPackingEncode_Int32_Wide(benchmark::State& state) {
   BM_DeltaBitPackingEncode<Int32Type>(state, MakeDeltaBitPackingInputWide<Int32Type>);
 }
@@ -725,6 +751,8 @@ BENCHMARK(BM_DeltaBitPackingEncode_Int32_Fixed)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingEncode_Int64_Fixed)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingEncode_Int32_Narrow)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingEncode_Int64_Narrow)->Range(MIN_RANGE, MAX_RANGE);
+BENCHMARK(BM_DeltaBitPackingEncode_Int32_NarrowSorted)->Range(MIN_RANGE, MAX_RANGE);
+BENCHMARK(BM_DeltaBitPackingEncode_Int64_NarrowSorted)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingEncode_Int32_Wide)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingEncode_Int64_Wide)->Range(MIN_RANGE, MAX_RANGE);
 
@@ -762,6 +790,16 @@ static void BM_DeltaBitPackingDecode_Int64_Narrow(benchmark::State& state) {
   BM_DeltaBitPackingDecode<Int64Type>(state, MakeDeltaBitPackingInputNarrow<Int64Type>);
 }
 
+static void BM_DeltaBitPackingDecode_Int32_NarrowSorted(benchmark::State& state) {
+  BM_DeltaBitPackingDecode<Int32Type>(state,
+                                      MakeDeltaBitPackingInputNarrowSorted<Int32Type>);
+}
+
+static void BM_DeltaBitPackingDecode_Int64_NarrowSorted(benchmark::State& state) {
+  BM_DeltaBitPackingDecode<Int64Type>(state,
+                                      MakeDeltaBitPackingInputNarrowSorted<Int64Type>);
+}
+
 static void BM_DeltaBitPackingDecode_Int32_Wide(benchmark::State& state) {
   BM_DeltaBitPackingDecode<Int32Type>(state, MakeDeltaBitPackingInputWide<Int32Type>);
 }
@@ -774,6 +812,8 @@ BENCHMARK(BM_DeltaBitPackingDecode_Int32_Fixed)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingDecode_Int64_Fixed)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingDecode_Int32_Narrow)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingDecode_Int64_Narrow)->Range(MIN_RANGE, MAX_RANGE);
+BENCHMARK(BM_DeltaBitPackingDecode_Int32_NarrowSorted)->Range(MIN_RANGE, MAX_RANGE);
+BENCHMARK(BM_DeltaBitPackingDecode_Int64_NarrowSorted)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingDecode_Int32_Wide)->Range(MIN_RANGE, MAX_RANGE);
 BENCHMARK(BM_DeltaBitPackingDecode_Int64_Wide)->Range(MIN_RANGE, MAX_RANGE);
 

cpp/src/parquet/encoding_test.cc

@@ -24,6 +24,7 @@
 #include <limits>
 #include <numeric>
 #include <span>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
@@ -2188,6 +2189,45 @@ TYPED_TEST(TestDeltaBitPackEncoding, MiniblockBitWidthRuns) {
   }
 }
 
+TYPED_TEST(TestDeltaBitPackEncoding, PrefixSumVectorAndTail) {
+  // A decoder may accumulate the running total several deltas at a time and finish the
+  // remainder one at a time. Walk every residual bit width at enough lengths to leave
+  // every remainder such a group can leave, so each width is decoded through the
+  // grouped path, through the remainder, and across the hand-off between them. A
+  // non-zero frame checks the running multiple, and wrapping the total checks that no
+  // term is signed.
+  using T = typename TypeParam::c_type;
+  using UT = std::make_unsigned_t<T>;
+  constexpr int kBits = static_cast<int>(sizeof(T) * 8);
+
+  auto make_values = [](int width, T frame, int num_deltas) {
+    std::vector<T> values;
+    values.reserve(num_deltas + 1);
+    const UT spread = width == kBits ? ~UT{0} : static_cast<UT>((UT{1} << width) - 1);
+    // Two deltas in three sit at the frame, so it is the smallest in every miniblock.
+    UT current = 0;
+    values.push_back(static_cast<T>(current));
+    for (int i = 0; i < num_deltas; ++i) {
+      current = static_cast<UT>(current + static_cast<UT>(frame) +
+                                (i % 3 == 0 ? spread : UT{0}));
+      values.push_back(static_cast<T>(current));
+    }
+    return values;
+  };
+
+  for (int width = 0; width <= kBits; ++width) {
+    for (const T frame : {T{0}, static_cast<T>(-5), T{7}}) {
+      // 16-23 leaves every remainder for group sizes up to eight; 201 crosses a block
+      // boundary with a remainder left.
+      for (const int num_deltas : {16, 17, 18, 19, 20, 21, 22, 23, 201}) {
+        ARROW_SCOPED_TRACE("width = ", width, ", frame = ", static_cast<int64_t>(frame),
+                           ", num_deltas = ", num_deltas);
+        this->CheckRoundtripWithValues(make_values(width, frame, num_deltas));
+      }
+    }
+  }
+}
+
 // ----------------------------------------------------------------------
 // PFOR encode/decode tests.