1// Numerically stable reductions (FP64 reference implementations).2// sum_naive is provided only to *measure* error; library code uses sum_neumaier or sum_pairwise.3#pragma once4#include <span>5#include <cmath>6#include <vector>7#include <algorithm>8#include "raid/core/types.hpp"9#include "raid/core/error.hpp"1011namespace raid::stats {1213template <class T>14T sum_naive(std::span<const T> x) { T s = 0; for (T v : x) s += v; return s; }1516/// Neumaier's improved Kahan–Babuška summation.17template <class T>18T sum_neumaier(std::span<const T> x) {19 T s = 0, c = 0;20 for (T v : x) {21 const T t = s + v;22 if (std::abs(s) >= std::abs(v)) c += (s - t) + v; else c += (v - t) + s;23 s = t;24 }25 return s + c;26}2728/// Pairwise (cascade) summation with a base-case block of 128; error O(log n) ulp.29template <class T>30T sum_pairwise(std::span<const T> x) {31 if (x.size() <= 128) { T s = 0; for (T v : x) s += v; return s; }32 const usize m = x.size() / 2;33 return sum_pairwise(x.subspan(0, m)) + sum_pairwise(x.subspan(m));34}3536template <class T>37T dot_neumaier(std::span<const T> a, std::span<const T> b) {38 RAID_REQUIRE(a.size() == b.size(), invalid_argument, "dot: size mismatch");39 T s = 0, c = 0;40 for (usize i = 0; i < a.size(); ++i) {41 const T v = a[i] * b[i], t = s + v;42 if (std::abs(s) >= std::abs(v)) c += (s - t) + v; else c += (v - t) + s;43 s = t;44 }45 return s + c;46}4748template <class T>49T sum_squares(std::span<const T> x) { return dot_neumaier(x, x); }5051template <class T>52T mean(std::span<const T> x) { RAID_REQUIRE(!x.empty(), insufficient_samples, "mean of empty span"); return sum_neumaier(x) / static_cast<T>(x.size()); }5354/// Welford one-pass mean/variance (ddof = 1 by default).55struct Welford {56 usize n{0}; f64 mean{0}; f64 m2{0};57 void add(f64 x) { ++n; const f64 d = x - mean; mean += d / static_cast<f64>(n); m2 += d * (x - mean); }58 f64 variance(int ddof = 1) const { return n > static_cast<usize>(ddof) ? m2 / static_cast<f64>(n - static_cast<usize>(ddof)) : std::numeric_limits<f64>::quiet_NaN(); }59 /// Chan–Golub–LeVeque parallel combine.60 void merge(const Welford& o) {61 if (o.n == 0) return; if (n == 0) { *this = o; return; }62 const f64 d = o.mean - mean; const usize nn = n + o.n;63 mean += d * static_cast<f64>(o.n) / static_cast<f64>(nn);64 m2 += o.m2 + d * d * static_cast<f64>(n) * static_cast<f64>(o.n) / static_cast<f64>(nn);65 n = nn;66 }67};6869template <class T>70T variance(std::span<const T> x, int ddof = 1) { Welford w; for (T v : x) w.add(static_cast<f64>(v)); return static_cast<T>(w.variance(ddof)); }7172struct Moments { usize n; f64 mean, variance, skewness, kurtosis /* excess */, min, max; };73/// Two-pass moments; skewness/kurtosis use the biased (population) normalisation g1, g2 − 3 (as scipy with bias=True).74Moments moments(std::span<const f64> x);7576/// Autocovariance at lag k with the 1/n normalisation (as statsmodels acovf(adjusted=False)).77f64 autocovariance(std::span<const f64> x, usize k, bool demean = true);78std::vector<f64> autocorrelation(std::span<const f64> x, usize max_lag);7980/// Maximum drawdown of a cumulative-return path (fractional, ≥ 0), and the indices of the peak and trough.81struct Drawdown { f64 max_drawdown; usize peak_index; usize trough_index; };82Drawdown max_drawdown(std::span<const f64> cumulative);8384} // namespace raid::stats85