PERF: Add SIMD instructions with xsimd to reduce moments

LICENSES/XSIMD_LICENSE

@@ -0,0 +1,29 @@
+Copyright (c) 2016, Johan Mabille, Sylvain Corlay, Wolf Vollprecht and Martin Renou
+Copyright (c) 2016, QuantStack
+Copyright (c) 2018, Serge Guelton
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+
+* 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.
+
+* Neither the name of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+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 HOLDER 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.

meson.build

@@ -10,7 +10,7 @@ project(
     default_options: [
         'buildtype=release',
         'c_std=c17',
-        'cpp_std=c++17',
+        'cpp_std=c++20',
         'warning_level=2',
     ],
 )
@@ -37,15 +37,17 @@ add_project_arguments(
 )
 
 cc = meson.get_compiler('c')
+cxx = meson.get_compiler('cpp')
 if cc.get_id() == 'msvc'
     # Tracking issue: https://github.com/pandas-dev/pandas/issues/63701
     # Ignore some MSVC specific warnings:
     # C4244: possible loss of data in conversion. Reproductible with `-Wconversion`.
     # C4267: conversion from `size_t` to smaller type.
     # C4551: occurs due to Cython generating code with (void)func.
     #        https://github.com/cython/cython/issues/3579
+    # C4146: unary minus operator applied to unsigned type. Occurs in xsimd.
     add_project_arguments(
-        ['/wd4244', '/wd4267', '/wd4551'],
+        ['/wd4244', '/wd4267', '/wd4551', '/wd4146'],
         language: ['c', 'cpp'],
     )
 endif

pandas/_libs/algos.pxd

@@ -1,11 +1,13 @@
 cimport cython
+from cython cimport size_t
 from libc.math cimport (
     NAN,
     sqrt,
 )
 from numpy cimport (
     float64_t,
     int64_t,
+    uint8_t,
 )
 
 from pandas._libs.dtypes cimport (
@@ -50,6 +52,22 @@ cdef inline void moments_add_value(
     mean[0] += delta_n
 
 
+cdef extern from "pandas/moments.h":
+    ctypedef struct Moments:
+        float64_t mean
+        float64_t m2
+        float64_t m3
+        float64_t m4
+        size_t n
+
+    Moments moments_reduce(
+            const double *values,
+            size_t n,
+            bint skipna,
+            const uint8_t *mask,
+            int max_moment) noexcept nogil
+
+
 @cython.cdivision(True)
 cdef inline float64_t calc_skew(
     int64_t nobs, float64_t m2, float64_t m3

pandas/_libs/algos.pyx

@@ -1603,9 +1603,9 @@ def diff_2d(
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef void accumulate_moments_scalar(
-    const float64_t[:] values,
+    const float64_t[::1] values,
     bint skipna,
-    const uint8_t[:] mask,
+    const uint8_t[::1] mask,
     int64_t* nobs,
     float64_t* mean,
     float64_t* m2,
@@ -1614,17 +1614,20 @@ cdef void accumulate_moments_scalar(
     int max_moment,
 ) noexcept nogil:
     cdef:
-        Py_ssize_t i, n = len(values)
-        bint uses_mask = mask is not None
-        float64_t val
+        Moments moments
+        const float64_t* values_ptr = &values[0]
+        const uint8_t* mask_ptr = &mask[0] if mask is not None else NULL
+        size_t n = <size_t>values.shape[0]
 
-    for i in range(n):
-        val = values[i]
-        if uses_mask and mask[i]:
-            val = NaN
-        if skipna and isnan(val):
-            continue
-        moments_add_value(val, nobs, mean, m2, m3, m4, max_moment)
+    moments = moments_reduce(values_ptr, n, skipna, mask_ptr, max_moment)
+    if max_moment >= 4:
+        m4[0] = moments.m4
+    if max_moment >= 3:
+        m3[0] = moments.m3
+
+    m2[0] = moments.m2
+    mean[0] = moments.mean
+    nobs[0] = <int64_t>moments.n
 
 
 @cython.boundscheck(False)
@@ -1676,9 +1679,9 @@ cdef void accumulate_moments_axis(
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def scalar_skew(
-    const float64_t[:] values,
+    const float64_t[::1] values,
     bint skipna,
-    const uint8_t[:] mask,
+    const uint8_t[::1] mask,
 ) -> float:
     cdef:
         int64_t nobs = 0
@@ -1693,9 +1696,9 @@ def scalar_skew(
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def scalar_kurt(
-    const float64_t[:] values,
+    const float64_t[::1] values,
     bint skipna,
-    const uint8_t[:] mask,
+    const uint8_t[::1] mask,
 ) -> float:
     cdef:
         int64_t nobs = 0

pandas/_libs/include/pandas/moments.h

@@ -0,0 +1,33 @@
+/*
+Copyright (c) 2026, PyData Development Team
+All rights reserved.
+
+Distributed under the terms of the BSD Simplified License.
+
+The full license is in the LICENSE file, distributed with this software.
+*/
+
+#pragma once
+
+#include <stdbool.h>
+#include <stddef.h>
+#include <stdint.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct {
+  double mean;
+  double m2;
+  double m3;
+  double m4;
+  size_t n;
+} Moments;
+
+/// Compute central moments until `max_moment` using `n` elements from `values`.
+Moments moments_reduce(const double *values, size_t n, bool skipna,
+                       const uint8_t *mask, int max_moment);
+#ifdef __cplusplus
+}
+#endif