GH-45676: [C++][Python][Compute] Add skew and kurtosis functions (#45677)

### What changes are included in this PR?

* Add "skew" and "kurtosis" scalar aggregate functions
* Add "hash_skew" and "hash_kurtosis" grouped aggregated functions
* Add support for decimal32/decimal64 to "variance", "stddev", "hash_variance" and "hash_stddev"
* Small optimization for "hash_variance" and "hash_stddev": avoid allocating the identity group id mapping when consuming an ExecBatch

### Are these changes tested?

Yes, by included unit tests.

### Are there any user-facing changes?

No, except new APIs.

* GitHub Issue: #45676

Authored-by: Antoine Pitrou <[email protected]>
Signed-off-by: Antoine Pitrou <[email protected]>

Author: pitrou

cpp/src/arrow/acero/hash_aggregate_test.cc

@@ -109,6 +109,9 @@ static auto kVarianceOptionsType = GetFunctionOptionsType<VarianceOptions>(
     DataMember("ddof", &VarianceOptions::ddof),
     DataMember("skip_nulls", &VarianceOptions::skip_nulls),
     DataMember("min_count", &VarianceOptions::min_count));
+static auto kSkewOptionsType = GetFunctionOptionsType<SkewOptions>(
+    DataMember("skip_nulls", &SkewOptions::skip_nulls),
+    DataMember("min_count", &SkewOptions::min_count));
 static auto kQuantileOptionsType = GetFunctionOptionsType<QuantileOptions>(
     DataMember("q", &QuantileOptions::q),
     DataMember("interpolation", &QuantileOptions::interpolation),
@@ -151,6 +154,11 @@ VarianceOptions::VarianceOptions(int ddof, bool skip_nulls, uint32_t min_count)
       min_count(min_count) {}
 constexpr char VarianceOptions::kTypeName[];
 
+SkewOptions::SkewOptions(bool skip_nulls, uint32_t min_count)
+    : FunctionOptions(internal::kSkewOptionsType),
+      skip_nulls(skip_nulls),
+      min_count(min_count) {}
+
 QuantileOptions::QuantileOptions(double q, enum Interpolation interpolation,
                                  bool skip_nulls, uint32_t min_count)
     : FunctionOptions(internal::kQuantileOptionsType),
@@ -203,6 +211,7 @@ void RegisterAggregateOptions(FunctionRegistry* registry) {
   DCHECK_OK(registry->AddFunctionOptionsType(kCountOptionsType));
   DCHECK_OK(registry->AddFunctionOptionsType(kModeOptionsType));
   DCHECK_OK(registry->AddFunctionOptionsType(kVarianceOptionsType));
+  DCHECK_OK(registry->AddFunctionOptionsType(kSkewOptionsType));
   DCHECK_OK(registry->AddFunctionOptionsType(kQuantileOptionsType));
   DCHECK_OK(registry->AddFunctionOptionsType(kTDigestOptionsType));
   DCHECK_OK(registry->AddFunctionOptionsType(kPivotOptionsType));
@@ -271,6 +280,14 @@ Result<Datum> Variance(const Datum& value, const VarianceOptions& options,
   return CallFunction("variance", {value}, &options, ctx);
 }
 
+Result<Datum> Skew(const Datum& value, const SkewOptions& options, ExecContext* ctx) {
+  return CallFunction("skew", {value}, &options, ctx);
+}
+
+Result<Datum> Kurtosis(const Datum& value, const SkewOptions& options, ExecContext* ctx) {
+  return CallFunction("kurtosis", {value}, &options, ctx);
+}
+
 Result<Datum> Quantile(const Datum& value, const QuantileOptions& options,
                        ExecContext* ctx) {
   return CallFunction("quantile", {value}, &options, ctx);

cpp/src/arrow/compute/api_aggregate.cc

@@ -114,6 +114,20 @@ class ARROW_EXPORT VarianceOptions : public FunctionOptions {
   uint32_t min_count;
 };
 
+/// \brief Control Skew and Kurtosis kernel behavior
+class ARROW_EXPORT SkewOptions : public FunctionOptions {
+ public:
+  explicit SkewOptions(bool skip_nulls = true, uint32_t min_count = 0);
+  static constexpr char const kTypeName[] = "SkewOptions";
+  static SkewOptions Defaults() { return SkewOptions{}; }
+
+  /// If true (the default), null values are ignored. Otherwise, if any value is null,
+  /// emit null.
+  bool skip_nulls;
+  /// If less than this many non-null values are observed, emit null.
+  uint32_t min_count;
+};
+
 /// \brief Control Quantile kernel behavior
 ///
 /// By default, returns the median value.
@@ -503,6 +517,34 @@ Result<Datum> Variance(const Datum& value,
                        const VarianceOptions& options = VarianceOptions::Defaults(),
                        ExecContext* ctx = NULLPTR);
 
+/// \brief Calculate the skewness of a numeric array
+///
+/// \param[in] value input datum, expecting Array or ChunkedArray
+/// \param[in] options see SkewOptions for more information
+/// \param[in] ctx the function execution context, optional
+/// \return datum of the computed skewness as a DoubleScalar
+///
+/// \since 20.0.0
+/// \note API not yet finalized
+ARROW_EXPORT
+Result<Datum> Skew(const Datum& value,
+                   const SkewOptions& options = SkewOptions::Defaults(),
+                   ExecContext* ctx = NULLPTR);
+
+/// \brief Calculate the kurtosis of a numeric array
+///
+/// \param[in] value input datum, expecting Array or ChunkedArray
+/// \param[in] options see SkewOptions for more information
+/// \param[in] ctx the function execution context, optional
+/// \return datum of the computed kurtosis as a DoubleScalar
+///
+/// \since 20.0.0
+/// \note API not yet finalized
+ARROW_EXPORT
+Result<Datum> Kurtosis(const Datum& value,
+                       const SkewOptions& options = SkewOptions::Defaults(),
+                       ExecContext* ctx = NULLPTR);
+
 /// \brief Calculate the quantiles of a numeric array
 ///
 /// \param[in] value input datum, expecting Array or ChunkedArray

cpp/src/arrow/compute/api_aggregate.h

@@ -17,16 +17,17 @@
 
 #pragma once
 
+#include <cmath>
+#include <initializer_list>
+
 #include "arrow/compute/kernels/util_internal.h"
 #include "arrow/type.h"
 #include "arrow/type_traits.h"
 #include "arrow/util/bit_run_reader.h"
 #include "arrow/util/int128_internal.h"
 #include "arrow/util/logging.h"
 
-namespace arrow {
-namespace compute {
-namespace internal {
+namespace arrow::compute::internal {
 
 // Find the largest compatible primitive type for a primitive type.
 template <typename I, typename Enable = void>
@@ -254,6 +255,4 @@ SumType SumArray(const ArraySpan& data) {
       data, [](ValueType v) { return static_cast<SumType>(v); });
 }
 
-}  // namespace internal
-}  // namespace compute
-}  // namespace arrow
+}  // namespace arrow::compute::internal

cpp/src/arrow/compute/kernels/aggregate_internal.h

@@ -42,6 +42,7 @@
 
 #include "arrow/testing/gtest_util.h"
 #include "arrow/testing/matchers.h"
+#include "arrow/testing/math.h"
 #include "arrow/testing/random.h"
 #include "arrow/util/logging.h"
 
@@ -3386,6 +3387,9 @@ TEST_F(TestVarStdKernelMergeStability, Basics) {
 #ifndef __MINGW32__  // MinGW has precision issues
   // XXX: The reference value from numpy is actually wrong due to floating
   // point limits. The correct result should equals variance(90, 0) = 4050.
+  // The problem is that the mean is not exactly representable as floating-point,
+  // and that small inaccuracy produces a large deviation when plugged into the M2
+  // calculation.
   std::vector<std::string> chunks = {"[40000008000000490]", "[40000008000000400]"};
   this->AssertVarStdIs(chunks, options, 3904.0);
 #endif
@@ -3430,12 +3434,21 @@ TEST_F(TestVarStdKernelUInt32, Basics) {
   this->AssertVarStdIs("[0, 0, 4294967295]", options, 6.148914688373205e+18);
 }
 
-// https://en.wikipedia.org/wiki/Kahan_summation_algorithm
 void KahanSum(double& sum, double& adjust, double addend) {
-  double y = addend - adjust;
-  double t = sum + y;
-  adjust = (t - sum) - y;
-  sum = t;
+  // Backported enhancement from Neumaier's algorithm: consider case where
+  // sum is small compared to addend.
+  // https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Further_enhancements
+  if (abs(sum) >= abs(addend)) {
+    double y = addend - adjust;
+    double t = sum + y;
+    adjust = (t - sum) - y;
+    sum = t;
+  } else {
+    double y = sum - adjust;
+    double t = addend + y;
+    adjust = (t - addend) - y;
+    sum = t;
+  }
 }
 
 // Calculate reference variance with Welford's online algorithm + Kahan summation
@@ -3534,7 +3547,8 @@ TEST_F(TestVarStdKernelIntegerLength, Basics) {
 
 TEST(TestVarStdKernel, Decimal) {
   // Effectively treated as double, sanity check results here
-  for (const auto& ty : {decimal128(3, 2), decimal256(3, 2)}) {
+  for (const auto& ty :
+       {decimal32(3, 2), decimal64(3, 2), decimal128(3, 2), decimal256(3, 2)}) {
     CheckVarStd(ArrayFromJSON(ty, R"(["1.00"])"), VarianceOptions(), 0);
     CheckVarStd(ArrayFromJSON(ty, R"([null, "1.00", "2.00", "3.00"])"), VarianceOptions(),
                 0.6666666666666666);
@@ -3544,6 +3558,154 @@ TEST(TestVarStdKernel, Decimal) {
   }
 }
 
+//
+// Skew and Kurtosis
+//
+
+constexpr int kSkewUlps = 3;
+constexpr int kKurtosisUlps = 6;
+
+void CheckSkewKurtosis(const Datum& array, const SkewOptions& options,
+                       double expected_skew, double expected_kurtosis, int n_ulps = -1) {
+  ARROW_SCOPED_TRACE("type = ", *array.type());
+  ASSERT_OK_AND_ASSIGN(Datum out_skew, Skew(array, options));
+  ASSERT_OK_AND_ASSIGN(Datum out_kurtosis, Kurtosis(array, options));
+  const auto& skew = checked_cast<const DoubleScalar&>(*out_skew.scalar());
+  const auto& kurtosis = checked_cast<const DoubleScalar&>(*out_kurtosis.scalar());
+  ASSERT_TRUE(skew.is_valid && kurtosis.is_valid);
+  AssertWithinUlp(expected_skew, skew.value, n_ulps >= 0 ? n_ulps : kSkewUlps);
+  AssertWithinUlp(expected_kurtosis, kurtosis.value,
+                  n_ulps >= 0 ? n_ulps : kKurtosisUlps);
+}
+
+class TestSkewKurtosis : public ::testing::Test {
+ public:
+  void AssertSkewKurtosisAre(const Array& array, const SkewOptions& options,
+                             double expected_skew, double expected_kurtosis,
+                             int n_ulps = -1) {
+    CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+  }
+
+  void AssertSkewKurtosisAre(const std::shared_ptr<ChunkedArray>& array,
+                             const SkewOptions& options, double expected_skew,
+                             double expected_kurtosis, int n_ulps = -1) {
+    CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+  }
+
+  void AssertSkewKurtosisAre(const std::shared_ptr<DataType>& type, std::string_view json,
+                             const SkewOptions& options, double expected_skew,
+                             double expected_kurtosis, int n_ulps = -1) {
+    auto array = ArrayFromJSON(type, json);
+    CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+  }
+
+  void AssertSkewKurtosisAre(const std::shared_ptr<DataType>& type,
+                             const std::vector<std::string>& json,
+                             const SkewOptions& options, double expected_skew,
+                             double expected_kurtosis, int n_ulps = -1) {
+    auto array = ChunkedArrayFromJSON(type, json);
+    CheckSkewKurtosis(array, options, expected_skew, expected_kurtosis, n_ulps);
+  }
+
+  void AssertSkewKurtosisInvalid(const Array& array, const SkewOptions& options) {
+    AssertSkewKurtosisInvalidInternal(array, options);
+  }
+
+  void AssertSkewKurtosisInvalid(const std::shared_ptr<ChunkedArray>& array,
+                                 const SkewOptions& options) {
+    AssertSkewKurtosisInvalidInternal(array, options);
+  }
+
+  void AssertSkewKurtosisInvalid(const std::shared_ptr<DataType>& type,
+                                 std::string_view json, const SkewOptions& options) {
+    auto array = ArrayFromJSON(type, json);
+    AssertSkewKurtosisInvalidInternal(array, options);
+  }
+
+  void AssertSkewKurtosisInvalid(const std::shared_ptr<DataType>& type,
+                                 const std::vector<std::string>& json,
+                                 const SkewOptions& options) {
+    auto array = ChunkedArrayFromJSON(type, json);
+    AssertSkewKurtosisInvalidInternal(array, options);
+  }
+
+ private:
+  void AssertSkewKurtosisInvalidInternal(const Datum& array, const SkewOptions& options) {
+    ASSERT_OK_AND_ASSIGN(Datum out_skew, Skew(array, options));
+    ASSERT_OK_AND_ASSIGN(Datum out_kurtosis, Kurtosis(array, options));
+    const auto& skew = checked_cast<const DoubleScalar&>(*out_skew.scalar());
+    const auto& kurtosis = checked_cast<const DoubleScalar&>(*out_kurtosis.scalar());
+    ASSERT_FALSE(skew.is_valid || kurtosis.is_valid);
+  }
+};
+
+TEST_F(TestSkewKurtosis, Basics) {
+  // Test sample from SciPy, with results obtained using numpy.float128
+  auto options = SkewOptions::Defaults();
+  AssertSkewKurtosisAre(float64(), "[1.165, 0.6268, 0.0751, 0.3516, -0.6965]", options,
+                        -0.29322304336607355496, -0.83411431970273759);
+  // Results are slightly different because the input doesn't losslessly convert
+  // to float32.
+  AssertSkewKurtosisAre(float32(), "[1.165, 0.6268, 0.0751, 0.3516, -0.6965]", options,
+                        -0.2932230870440958164, -0.8341143229437093939);
+}
+
+TEST_F(TestSkewKurtosis, Chunked) {
+  auto options = SkewOptions::Defaults();
+  AssertSkewKurtosisAre(float64(), {"[1.165, 0.6268]", "[]", "[0.0751, 0.3516, -0.6965]"},
+                        options, -0.29322304336607355496, -0.83411431970273759);
+  AssertSkewKurtosisAre(float32(), {"[1.165, 0.6268]", "[]", "[0.0751, 0.3516, -0.6965]"},
+                        options, -0.2932230870440958164, -0.8341143229437093939);
+}
+
+TEST_F(TestSkewKurtosis, Decimal) {
+  auto options = SkewOptions::Defaults();
+  for (auto type :
+       {decimal32(5, 4), decimal64(5, 4), decimal128(5, 4), decimal256(5, 4)}) {
+    AssertSkewKurtosisAre(type, R"(["1.1650", "0.6268", "0.0751", "0.3516", "-0.6965"])",
+                          options, -0.29322304336607355496, -0.83411431970273759);
+  }
+}
+
+TEST_F(TestSkewKurtosis, Integral) {
+  auto options = SkewOptions::Defaults();
+  for (auto type : IntTypes()) {
+    AssertSkewKurtosisAre(type, "[1, 2, 3, 5]", options, 0.4346507595746657,
+                          -1.1542857142857144);
+  }
+}
+
+TEST_F(TestSkewKurtosis, SpecialCases) {
+  auto options = SkewOptions::Defaults();
+  for (auto type : {float64(), float32()}) {
+    AssertSkewKurtosisAre(type, "[0, 1, 2]", options, 0.0, -1.5, /*n_ulps=*/0);
+    AssertSkewKurtosisAre(type, "[1]", options, std::nan(""), std::nan(""));
+    AssertSkewKurtosisAre(type, "[1, 1, 1, 1, 1, 1]", options, std::nan(""),
+                          std::nan(""));
+  }
+}
+
+TEST_F(TestSkewKurtosis, Options) {
+  for (auto type : {float64(), float32()}) {
+    auto options = SkewOptions::Defaults();
+    AssertSkewKurtosisInvalid(type, "[]", options);
+    AssertSkewKurtosisInvalid(type, std::vector<std::string>{}, options);
+    AssertSkewKurtosisInvalid(type, {"[]", "[]", "[]"}, options);
+    AssertSkewKurtosisAre(type, "[0, 1, null, 2]", options, 0.0, -1.5);
+    AssertSkewKurtosisAre(type, {"[0, 1]", "[]", "[null, 2]"}, options, 0.0, -1.5);
+    options.min_count = 3;
+    AssertSkewKurtosisAre(type, "[0, 1, null, 2]", options, 0.0, -1.5);
+    AssertSkewKurtosisAre(type, {"[0, 1]", "[]", "[null, 2]"}, options, 0.0, -1.5);
+    options.skip_nulls = false;
+    AssertSkewKurtosisInvalid(type, "[0, 1, null, 2]", options);
+    AssertSkewKurtosisInvalid(type, {"[0, 1]", "[]", "[null, 2]"}, options);
+    options.skip_nulls = true;
+    options.min_count = 4;
+    AssertSkewKurtosisInvalid(type, "[0, 1, null, 2]", options);
+    AssertSkewKurtosisInvalid(type, {"[0, 1]", "[]", "[null, 2]"}, options);
+  }
+}
+
 //
 // Quantile
 //