ILKernelGenerator.Reduction.cs

using System;
using System.Collections.Concurrent;
using System.Globalization;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;
using NumSharp.Backends.Iteration;
using NumSharp.Utilities;

// =============================================================================
// ILKernelGenerator.Reduction.cs — per-chunk axis-reduction kernels (NDIter-driven)
// =============================================================================
//
// MODEL
// -----
// These kernels implement the inner loop of a 2-operand REDUCE iterator built by
// NDIterRef.NewReduce ([input, output], REDUCE_OK | EXTERNAL_LOOP). The iterator
// owns the outer walk and advances the operand pointers between calls; each kernel
// invocation folds ONE contiguous inner stripe into the output:
//
//     do { inner(dataptrs, strides, count, aux); } while (iternext(iter));
//
// Two output modes per call, distinguished by the output byte stride:
//   • outStride == 0  → PINNED: the reduced axis is innermost. Fold the whole
//     `count`-element stripe into the single output slot (read it, accumulate,
//     write back once).
//   • outStride != 0  → SLAB: a kept axis is innermost. Each inner element targets
//     a distinct output slot; fold `in[c]` into `out[c]`. The same slab is
//     revisited across outer iterations, so the output MUST be pre-seeded with the
//     reduction identity (DefaultEngine seeds it via SeedReduceIdentity).
//
// COMPLEX (Phase 1)
// -----------------
// System.Numerics.Complex is two contiguous doubles [real, imaginary]. A complex
// Sum slab-fold is therefore a plain f64 elementwise add over 2·count lanes, and a
// complex Sum pinned-reduce is a 2-lane (re,im) f64 accumulator — both ride the
// Vector256<double> machinery. Prod uses scalar complex multiply (no SIMD form);
// Min/Max use lexicographic (Real, Imaginary) pick with NaN-first-wins, matching
// DirectILKernelGenerator's existing ComplexLexPick semantics exactly.
//
// This is the migration TARGET model (see ILKernelGenerator.cs). The legacy
// whole-array complex path lives in DirectILKernelGenerator.Reduction.Axis.cs and
// is retired per-dtype as families move here.
// =============================================================================

namespace NumSharp.Backends.Kernels
{
    public static partial class ILKernelGenerator
    {
        /// <summary>
        /// Cache key for a per-chunk reduction kernel: operation + input dtype +
        /// accumulator (output) dtype. Layout is handled at runtime inside the
        /// kernel body (pinned vs slab, contiguous vs strided inner loop), so it is
        /// NOT part of the key.
        /// </summary>
        public readonly record struct ReduceKernelKey(ReductionOp Op, NPTypeCode InType, NPTypeCode AccType);

        internal static readonly ConcurrentDictionary<ReduceKernelKey, NDInnerLoopFunc> _reduceCache = new();

        /// <summary>
        /// Returns the cached per-chunk reduction kernel for the given
        /// (op, input, accumulator) triple, or null when no NDIter-driven kernel
        /// exists yet (caller falls back to the DirectILKernelGenerator path).
        /// The returned delegate matches <see cref="NDInnerLoopFunc"/>; hand it to
        /// an iterator built by <see cref="NDIterRef.NewReduce(NDArray, NDArray, int, NDIterGlobalFlags)"/>.
        /// </summary>
        public static NDInnerLoopFunc GetReduceInnerLoop(ReduceKernelKey key) =>
            _reduceCache.GetOrAdd(key, CreateReduceInnerLoop);

        private static unsafe NDInnerLoopFunc CreateReduceInnerLoop(ReduceKernelKey key)
        {
            // Complex: dedicated double-pair kernels (SIMD sum, lex min/max).
            if (key.InType == NPTypeCode.Complex && key.AccType == NPTypeCode.Complex)
            {
                // Sum: prefer the IL-emitted Vector128 pairwise kernel (bit-exact NumPy
                // complex128; the hand-written ComplexSumKernel is a flat accumulator that
                // diverges in the low bits). Fall back to it only when emission is unavailable.
                if (key.Op == ReductionOp.Sum)
                {
                    var emitted = TryEmitPairwiseSumKernel(NPTypeCode.Complex);
                    if (emitted != null) return emitted;
                }
                return key.Op switch
                {
                    ReductionOp.Sum  => ComplexSumKernel,
                    ReductionOp.Prod => ComplexProdKernel,
                    ReductionOp.Min  => ComplexMinKernel,
                    ReductionOp.Max  => ComplexMaxKernel,
                    _ => null
                };
            }

            // Half mean / Decimal: generic INumber scalar kernels fed CONTIGUOUS stripes by the
            // per-chunk iterator (no cache-hostile column gather like the legacy coordinate walk).
            //   - Half MEAN accumulates in Double (input Half → Double); ReduceMean casts back to
            //     Half. (Half sum/prod/min/max stay on the Direct path — see UseNDIterReduce:
            //     their Half-accumulator serial chain can't beat it on the pinned axis.)
            //   - Decimal accumulates in full-precision Decimal (no NumPy reference type).
            //     CreateTypedReduceKernel can build a same-type Half kernel too; it just isn't
            //     routed here today.
            // Half SUM: accumulate in a FLOAT32 shadow (ReduceAdd allocates Single, narrows to Half
            // once at the end). The kernel reproduces NumPy's per-inner-loop narrow via an in-place
            // RoundToF16 (== widen(narrow), no half↔memory round-trip): PINNED folds the stripe in
            // float32 pairwise then RoundToF16; SLAB accumulates across outer steps with a per-step
            // RoundToF16 (saturates, matching NumPy). See ILKernelGenerator.Reduction.Half.cs. The
            // (Half,Double) route below stays for MEAN (its internal sum accumulates in Double).
            if (key.InType == NPTypeCode.Half && key.AccType == NPTypeCode.Single)
                return key.Op == ReductionOp.Sum ? (NDInnerLoopFunc)HalfSumKernel : null;
            if (key.InType == NPTypeCode.Half && key.AccType == NPTypeCode.Double)
                return key.Op == ReductionOp.Sum ? CreateTypedReduceKernel<Half, double>(ReductionOp.Sum) : null;
            if (key.InType == NPTypeCode.Decimal && key.AccType == NPTypeCode.Decimal)
                return CreateTypedReduceKernel<decimal, decimal>(key.Op);

            // Phase 6 — numeric same-type Sum/Mean (PairwiseFold, bit-for-bit NumPy
            // pairwise_sum) AND Min/Max (SimdMinMaxSameType, NaN-propagating). Double AND
            // float32 route here so f64/f32 reductions ride the stride-ordered NDIter reduce
            // path — SIMD on EVERY layout (transpose/F/strided/negcol/broadcast), fixing the
            // C-contiguity-gated Direct axis kernel's collapse on those layouts. Integer
            // widening sums and Prod stay on the Direct path (CreateSimdReduceKernel returns
            // null for those → caller falls back).
            if (key.InType == key.AccType && (key.InType == NPTypeCode.Double || key.InType == NPTypeCode.Single))
            {
                // Sum: prefer the IL-EMITTED SIMD pairwise kernel — bit-for-bit identical to
                // NumPy's pairwise_sum yet width-native SIMD (the scalar PairwiseFold below
                // forfeits vectorization; see ILKernelGenerator.Reduction.Pairwise.cs). Falls
                // back to the generic scalar fold only if emission is unavailable (IL disabled).
                if (key.Op == ReductionOp.Sum)
                {
                    var emitted = TryEmitPairwiseSumKernel(key.InType);
                    if (emitted != null) return emitted;
                }
                // Min/Max: width-native SIMD with explicit per-lane NaN propagation
                // (Vector256.Min/Max alone do not propagate NaN). Same-type, no widening.
                else if (key.Op == ReductionOp.Min)
                    return key.InType == NPTypeCode.Double ? (NDInnerLoopFunc)SimdMinKernel<double> : SimdMinKernel<float>;
                else if (key.Op == ReductionOp.Max)
                    return key.InType == NPTypeCode.Double ? (NDInnerLoopFunc)SimdMaxKernel<double> : SimdMaxKernel<float>;
                if (key.InType == NPTypeCode.Double) return CreateSimdReduceKernel<double>(key.Op);
                return CreateSimdReduceKernel<float>(key.Op);
            }

            return null;
        }

        /// <summary>
        /// Builds a scalar per-chunk reduce kernel over arbitrary numeric
        /// <typeparamref name="TIn"/> → <typeparamref name="TAccum"/> via .NET generic math.
        /// The JIT monomorphizes one tight body per (TIn,TAccum,op); CreateTruncating folds to
        /// a reinterpret/convert and the arithmetic to native ops. Min/Max propagate NaN
        /// (NumPy parity) via <c>TAccum.IsNaN</c>; Decimal has no NaN so those checks fold away.
        /// </summary>
        private static unsafe NDInnerLoopFunc CreateTypedReduceKernel<TIn, TAccum>(ReductionOp op)
            where TIn : unmanaged, INumberBase<TIn>
            where TAccum : unmanaged, INumber<TAccum>
        {
            return op switch
            {
                ReductionOp.Sum  => TypedSumKernel<TIn, TAccum>,
                ReductionOp.Prod => TypedProdKernel<TIn, TAccum>,
                ReductionOp.Min  => TypedMinKernel<TIn, TAccum>,
                ReductionOp.Max  => TypedMaxKernel<TIn, TAccum>,
                _ => null
            };
        }

        // =====================================================================
        // Generic same-type SUM kernels (Phase 6 — numeric: double/float/…)
        // =====================================================================
        //
        // For TIn == TAccum (no NEP50 widening: float/double Sum/Mean) one generic
        // body monomorphizes per dtype. The two reduce orientations map exactly onto
        // NumPy's two summation behaviors (loops_arithm_fp.dispatch.c.src reduce
        // branch + loops_utils.h.src pairwise_sum):
        //   • PINNED (reduced axis is the contiguous inner loop): NumPy folds the
        //     stripe with pairwise_sum. We do the same — PairwiseFold below is ported
        //     1:1 and is BIT-FOR-BIT identical to np.add.reduce for float/double
        //     (validated).
        //   • SLAB (a kept axis is the contiguous inner loop; reduced axis is outer):
        //     NumPy accumulates rows sequentially (out[c] += in[c]); so do we (the
        //     Vector256 streaming add), which already bit-matches NumPy on that
        //     orientation. Accuracy is orientation-dependent — exactly like NumPy.
        // NaN propagates naturally (a NaN element makes the whole reduction NaN).
        //
        // The pairwise PINNED path is what makes float32 sum/mean SAFE to route
        // (its earlier exclusion was due to a flat 8-accumulator divergence of ~24;
        // pairwise removes that — exact parity). NDIter supplies the substrate:
        // NewReduce orders the contiguous axis innermost (→ PINNED for a contiguous
        // reduced axis), EXTERNAL_LOOP delivers the full stripe as (ptr, count,
        // stride), and the +0 identity seed + slot-accumulate reproduce NumPy's
        // `*acc += pairwise_sum(...)` (so buffered/chunked stripes stay correct).

        private const int PairwiseBlock = 128; // NumPy PW_BLOCKSIZE

        /// <summary>
        /// NumPy's pairwise_sum (loops_utils.h.src) ported 1:1; <paramref name="stride"/>
        /// is in ELEMENTS. n&lt;8 naive (seed -0 to preserve signed zero); n≤128 eight
        /// independent accumulators unrolled by 8 (+ software prefetch on the contiguous
        /// path) then the exact tree-combine; n&gt;128 split (kept a multiple of 8) and
        /// recurse. Bit-for-bit identical to <c>np.add.reduce</c> for float/double.
        /// </summary>
        /// <remarks>
        /// Scalar: the IL-emitted SIMD folds (<see cref="TryGetPairwiseFoldDouble"/> /
        /// <see cref="TryGetPairwiseFoldSingle"/>) compute the same bits faster and are preferred. This one is the
        /// fallback where those are unavailable (no SIMD, no dynamic code), which is why it is internal: the flat
        /// add.reduce engine (<c>NumPyAddReduce</c>) falls back to it on such hosts.
        /// </remarks>
        /// <typeparam name="T">float or double — the accumulator is the element type, as in NumPy's FLOAT / DOUBLE loops.</typeparam>
        /// <param name="a">Address of logical element 0.</param>
        /// <param name="n">Element count (0 returns -0.0).</param>
        /// <param name="stride">Distance between consecutive elements, in ELEMENTS (may be negative).</param>
        /// <returns>The pairwise sum, WITHOUT the reduce loop's addition to the identity.</returns>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        internal static unsafe T PairwiseFold<T>(T* a, long n, long stride)
            where T : unmanaged, INumber<T>
        {
            if (n < 8)
            {
                T res = -T.Zero; // preserve -0.0 (NumPy: summing only -0 must stay -0)
                for (long i = 0; i < n; i++) res += a[i * stride];
                return res;
            }
            if (n <= PairwiseBlock)
            {
                // 8 independent accumulators (NumPy's r[0..7]); they break the FP-add
                // dependency chain for ILP. The JIT schedules the 8 independent adds; for the
                // routed axis case the stripe recurses into many small leaves so the cost is
                // recursion + combine, not this inner loop (an explicit Vector256 lane form
                // measured identical, so we keep the simpler scalar body). Stride is in elements.
                T r0 = a[0],          r1 = a[stride],     r2 = a[2 * stride], r3 = a[3 * stride],
                  r4 = a[4 * stride], r5 = a[5 * stride], r6 = a[6 * stride], r7 = a[7 * stride];
                long i;
                bool pf = stride == 1 && Sse.IsSupported;
                long ahead = 512 / sizeof(T);
                for (i = 8; i < n - (n % 8); i += 8)
                {
                    if (pf) Sse.Prefetch0((void*)(a + i + ahead)); // NPY_PREFETCH equivalent
                    r0 += a[(i + 0) * stride]; r1 += a[(i + 1) * stride];
                    r2 += a[(i + 2) * stride]; r3 += a[(i + 3) * stride];
                    r4 += a[(i + 4) * stride]; r5 += a[(i + 5) * stride];
                    r6 += a[(i + 6) * stride]; r7 += a[(i + 7) * stride];
                }
                T res = ((r0 + r1) + (r2 + r3)) + ((r4 + r5) + (r6 + r7));
                for (; i < n; i++) res += a[i * stride];
                return res;
            }
            long n2 = n / 2; n2 -= n2 % 8;
            return PairwiseFold(a, n2, stride) + PairwiseFold(a + n2 * stride, n - n2, stride);
        }

        /// <summary>
        /// Builds a generic same-type (<typeparamref name="T"/> == accumulator) per-chunk
        /// reduce kernel. Currently only Sum is emitted (PairwiseSumSameType); Min/Max/Prod
        /// return null so the caller falls back to the Direct path until they are migrated.
        /// </summary>
        private static unsafe NDInnerLoopFunc CreateSimdReduceKernel<T>(ReductionOp op)
            where T : unmanaged, INumber<T>
        {
            return op switch
            {
                ReductionOp.Sum => PairwiseSumSameType<T>,
                _ => null
            };
        }

        /// <summary>
        /// Same-type sum, dual mode. PINNED (reduced axis inner): pairwise-fold the stripe
        /// and accumulate into the slot — bit-for-bit NumPy parity. SLAB (kept axis inner):
        /// a 4×-unrolled Vector256 elementwise add of the contiguous row into the output
        /// accumulator row (sequential over the outer reduced axis, matching NumPy). Strided
        /// inner loops fall to the scalar/pairwise-with-stride path. NaN propagates naturally.
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static unsafe void PairwiseSumSameType<T>(void** dataptrs, long* strides, long count, void* auxdata)
            where T : unmanaged, INumber<T>
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];
            int sz = sizeof(T);

            if (outS == 0)
            {
                // PINNED: pairwise-fold the stripe into the slot. The slot is seeded to the
                // +0 identity, so a single stripe gives *out = pairwise(stripe) == NumPy;
                // chunked/buffered stripes accumulate (NumPy's `*acc += pairwise_sum(...)`).
                *(T*)outp += PairwiseFold((T*)inp, count, inS / sz);
            }
            else
            {
                // SLAB: out[c] += in[c].
                if (inS == sz && outS == sz)
                {
                    T* id = (T*)inp; T* od = (T*)outp; long i = 0; int W = Vector256<T>.Count;
                    if (Vector256.IsHardwareAccelerated)
                    {
                        long lim = count - count % (W * 4);
                        for (; i < lim; i += W * 4)
                        {
                            Vector256.Store(Vector256.Add(Vector256.Load(od + i),         Vector256.Load(id + i)),         od + i);
                            Vector256.Store(Vector256.Add(Vector256.Load(od + i + W),     Vector256.Load(id + i + W)),     od + i + W);
                            Vector256.Store(Vector256.Add(Vector256.Load(od + i + W * 2), Vector256.Load(id + i + W * 2)), od + i + W * 2);
                            Vector256.Store(Vector256.Add(Vector256.Load(od + i + W * 3), Vector256.Load(id + i + W * 3)), od + i + W * 3);
                        }
                        for (; i + W <= count; i += W)
                            Vector256.Store(Vector256.Add(Vector256.Load(od + i), Vector256.Load(id + i)), od + i);
                    }
                    for (; i < count; i++) od[i] += id[i];
                }
                else
                {
                    for (long k = 0; k < count; k++) *(T*)(outp + k * outS) += *(T*)(inp + k * inS);
                }
            }
        }

        /// <summary>
        /// Pre-fill <paramref name="output"/> with the reduction identity for
        /// <paramref name="op"/> before driving a REDUCE iterator. Required because
        /// the per-chunk kernels fold into the existing output slot(s). Any output
        /// layout (contiguous fresh alloc or user-supplied view) is honored.
        /// </summary>
        /// <remarks>
        /// A writeable output whose elements occupy one dense memory block (C- or F-contiguous —
        /// every fresh reduction result, and most user <c>out=</c> arrays) takes the fast path: the
        /// identity is written to logical element 0 through <see cref="NDArray.SetAtIndex(object, long)"/>
        /// — so its exact unboxing conversion and read-only check are unchanged — and that element's
        /// bytes are then replicated across the block with a vectorized fill. The per-element boxed
        /// <c>SetAtIndex</c> loop it replaces cost ~2 ns per output (a 3×1M axis-0 sum spent ~2 ms
        /// seeding 1M slots). A strided or read-only output keeps that loop (which raises the same
        /// read-only error at its first write).
        /// </remarks>
        /// <param name="output">The reduction output to seed; its dtype selects the identity.</param>
        /// <param name="op">The reduction whose identity is written (Sum/Prod/Min/Max, Any/All for bool).</param>
        /// <exception cref="NotSupportedException"><paramref name="op"/> has no identity for the output dtype (Complex supports Sum/Prod/Min/Max only; see <see cref="ReductionOpExtensions.GetIdentity"/>).</exception>
        /// <exception cref="InvalidCastException">The boxed identity does not unbox to the output's element type (a caller pairing the op with a foreign dtype).</exception>
        public static unsafe void SeedReduceIdentity(NDArray output, ReductionOp op)
        {
            long n = output.size;
            if (n == 0) return;

            NPTypeCode tc = output.GetTypeCode;
            object id;
            if (tc == NPTypeCode.Complex)
            {
                id = op switch
                {
                    ReductionOp.Sum  => System.Numerics.Complex.Zero,
                    ReductionOp.Prod => System.Numerics.Complex.One,
                    // (±inf, ±inf) so the first finite element displaces the identity
                    // under lexicographic comparison — mirrors GetIdentityValueTyped<Complex>.
                    // (GetIdentity's GetMaxValue gives (+inf,0), which is NOT the lex-largest.)
                    ReductionOp.Min  => new System.Numerics.Complex(double.PositiveInfinity, double.PositiveInfinity),
                    ReductionOp.Max  => new System.Numerics.Complex(double.NegativeInfinity, double.NegativeInfinity),
                    _ => throw new NotSupportedException($"SeedReduceIdentity: op {op} unsupported for Complex")
                };
            }
            else
            {
                // Scalar numerics (Half/Double/Decimal/...): the shared identity table is correct
                // for total-ordering min/max (min seed = +inf/MaxValue, max seed = -inf/MinValue).
                id = op.GetIdentity(tc);
            }

            if (TryGetDenseOutputBlock(output, out byte* block, out int itemSize))
            {
                // Element 0 goes through the boxed setter (its conversion + writeable check are the
                // contract); the remaining n-1 slots receive a copy of its bytes.
                output.SetAtIndex(id, 0);
                ReplicateFirstElement(block, n, itemSize);
                return;
            }

            for (long i = 0; i < n; i++) output.SetAtIndex(id, i);
        }

        /// <summary>
        /// Divide every element of <paramref name="output"/> by <paramref name="count"/> in
        /// place — the post-pass that turns an accumulated axis Sum into a Mean. For Complex
        /// this is NumPy's complex division by <c>count+0j</c> (np.mean's <c>true_divide(sum, n)</c>
        /// runs the CDOUBLE_divide loop, Smith's algorithm: <c>((re + im*0)*(1/n), (im - re*0)*(1/n))</c>,
        /// see <see cref="DefaultEngine.ComplexDivideByCountLikeNumPy"/>), so an
        /// infinite or NaN part poisons the other one and the rounding is NumPy's; it used to divide each
        /// component by the count. Any output layout is honored.
        /// </summary>
        /// <remarks>
        /// A writeable dense-block output (C- or F-contiguous) is divided through typed pointers — a
        /// <see cref="Vector{T}"/> loop for float64/float32, a scalar loop for complex/decimal — with the
        /// IDENTICAL arithmetic of the boxed loop (<c>x / (double)count</c>, <c>x / (float)count</c>,
        /// NumPy's complex division, decimal): IEEE division is correctly rounded per lane, so the vector
        /// and scalar forms are bit-identical. The boxed <c>GetAtIndex</c>/<c>SetAtIndex</c> loop it
        /// replaces cost ~4 ns per output (np.mean(axis=0) of 3×1M spent ~4 of its 6.7 ms here). A
        /// strided or read-only output keeps the boxed loop.
        /// </remarks>
        /// <param name="output">The accumulated sums, divided in place.</param>
        /// <param name="count">The reduced-axis length (0 leaves the output untouched, as before).</param>
        /// <exception cref="NotSupportedException">The output dtype is not Complex/Double/Single/Decimal (the only mean accumulators).</exception>
        public static void MeanDivideByCount(NDArray output, long count)
        {
            long n = output.size;
            if (n == 0 || count == 0) return;

            NPTypeCode tc = output.GetTypeCode;
            if (TryDivideDenseBlockByCount(output, tc, n, count))
                return;

            switch (tc)
            {
                case NPTypeCode.Complex:
                {
                    // NumPy's complex division by n+0j (see the dense-block twin below), not component-wise.
                    for (long i = 0; i < n; i++)
                    {
                        var c = (System.Numerics.Complex)output.GetAtIndex(i);
                        output.SetAtIndex(DefaultEngine.ComplexDivideByCountLikeNumPy(c, count), i);
                    }
                    return;
                }
                case NPTypeCode.Double: // Half mean accumulates here (Half→Double); ReduceMean casts back.
                {
                    double d = count;
                    for (long i = 0; i < n; i++)
                        output.SetAtIndex((double)output.GetAtIndex(i) / d, i);
                    return;
                }
                case NPTypeCode.Single: // mean(float32) stays float32 (NumPy: f32 sum / count → f32).
                {
                    float d = count;
                    for (long i = 0; i < n; i++)
                        output.SetAtIndex((float)output.GetAtIndex(i) / d, i);
                    return;
                }
                case NPTypeCode.Decimal:
                {
                    decimal d = count;
                    for (long i = 0; i < n; i++)
                        output.SetAtIndex((decimal)output.GetAtIndex(i) / d, i);
                    return;
                }
                default:
                    throw new NotSupportedException($"MeanDivideByCount not implemented for {tc}");
            }
        }

        /// <summary>
        /// Whether a reduction output's elements occupy ONE dense memory block that may be written
        /// through a raw pointer: writeable, not broadcast, and C- or F-contiguous (every fresh
        /// reduction result qualifies; a strided or read-only user <c>out=</c> does not). For such an
        /// output the logical element 0 is the block's first element and the n elements are the
        /// block's n consecutive slots — in an order that does not matter to the whole-block
        /// elementwise post-passes that use this (a fill, a divide).
        /// </summary>
        /// <param name="output">The reduction output.</param>
        /// <param name="block">Receives the address of the block's first element (base + view offset), or null.</param>
        /// <param name="itemSize">Receives the element size in bytes.</param>
        /// <returns>True when the output is a writeable dense block; false keeps the caller on its per-element path.</returns>
        private static unsafe bool TryGetDenseOutputBlock(NDArray output, out byte* block, out int itemSize)
        {
            var s = output.Shape;
            itemSize = output.GetTypeCode.SizeOf();
            // IsWriteable is the same flag the boxed setter's read-only check reads, so a read-only
            // output stays on the path that raises; a broadcast view is never a dense block.
            if (s.IsWriteable && !s.IsBroadcasted && (s.IsContiguous || s.IsFContiguous))
            {
                block = (byte*)output.Address + s.offset * itemSize;
                return true;
            }

            block = null;
            return false;
        }

        /// <summary>
        /// Copy the bytes of the element at <paramref name="block"/> into the following
        /// <paramref name="n"/> − 1 slots of a dense block — a width-typed, vectorized fill (the
        /// element is replicated as raw bits, so any dtype of the width round-trips exactly).
        /// </summary>
        /// <param name="block">The dense block; its first element holds the value to replicate.</param>
        /// <param name="n">The number of elements in the block (≥ 1).</param>
        /// <param name="itemSize">The element size in bytes (1, 2, 4, 8 or 16 — every NumSharp dtype).</param>
        /// <exception cref="NotSupportedException"><paramref name="itemSize"/> is not a NumSharp element width.</exception>
        private static unsafe void ReplicateFirstElement(byte* block, long n, int itemSize)
        {
            switch (itemSize)
            {
                case 1: FillWith(block + 1, n - 1, *block); return;
                case 2: FillWith((ushort*)block + 1, n - 1, *(ushort*)block); return;
                case 4: FillWith((uint*)block + 1, n - 1, *(uint*)block); return;
                case 8: FillWith((ulong*)block + 1, n - 1, *(ulong*)block); return;
                case 16: FillWith((Block16*)block + 1, n - 1, *(Block16*)block); return;
                default: throw new NotSupportedException(string.Create(CultureInfo.InvariantCulture, $"ReplicateFirstElement: element size {itemSize}"));
            }
        }

        /// <summary>
        /// Fill <paramref name="count"/> elements at <paramref name="p"/> with <paramref name="value"/>
        /// through <see cref="Span{T}.Fill"/> (vectorized for primitive widths), in <see cref="int"/>-sized
        /// chunks so a block past 2^31 elements is filled completely.
        /// </summary>
        /// <typeparam name="T">The element type (width carrier).</typeparam>
        /// <param name="p">The first slot to fill.</param>
        /// <param name="count">The number of slots (≤ 0 fills nothing).</param>
        /// <param name="value">The value every slot receives.</param>
        private static unsafe void FillWith<T>(T* p, long count, T value) where T : unmanaged
        {
            while (count > 0)
            {
                int chunk = (int)Math.Min(count, int.MaxValue);
                new Span<T>(p, chunk).Fill(value);
                p += chunk;
                count -= chunk;
            }
        }

        /// <summary>
        /// A 16-byte raw carrier (Complex, Decimal) for <see cref="ReplicateFirstElement"/>: copying it
        /// moves bits and nothing else, where a <see cref="decimal"/> or <see cref="System.Numerics.Complex"/>
        /// carrier would suggest value semantics the fill does not have.
        /// </summary>
        private struct Block16
        {
            /// <summary>The low eight bytes.</summary>
            public ulong Lo;

            /// <summary>The high eight bytes.</summary>
            public ulong Hi;
        }

        /// <summary>
        /// <see cref="MeanDivideByCount"/>'s fast path: divide a writeable dense-block output by the count
        /// through typed pointers with the SAME arithmetic as the boxed loop — float64 <c>x / (double)count</c>
        /// and float32 <c>x / (float)count</c> over <see cref="Vector{T}"/> lanes (IEEE division is correctly
        /// rounded per lane, so bit-identical to the scalar form, NaN payloads included), complex
        /// through NumPy's complex division by <c>count+0j</c>, decimal <c>x / (decimal)count</c>.
        /// </summary>
        /// <param name="output">The accumulated sums.</param>
        /// <param name="tc">The output dtype.</param>
        /// <param name="n">The output element count (≥ 1).</param>
        /// <param name="count">The reduced-axis length (≥ 1).</param>
        /// <returns>True when divided here; false when the output is not a writeable dense block or the dtype is not a mean accumulator (the caller's boxed loop then runs, raising for an unsupported dtype).</returns>
        private static unsafe bool TryDivideDenseBlockByCount(NDArray output, NPTypeCode tc, long n, long count)
        {
            if (!TryGetDenseOutputBlock(output, out byte* block, out _))
                return false;

            switch (tc)
            {
                case NPTypeCode.Double:
                {
                    double* p = (double*)block;
                    double d = count;
                    long i = 0;
                    if (Vector.IsHardwareAccelerated)
                    {
                        var vd = new Vector<double>(d);
                        int w = Vector<double>.Count;
                        for (; i <= n - w; i += w)
                            Unsafe.WriteUnaligned(p + i, Unsafe.ReadUnaligned<Vector<double>>(p + i) / vd);
                    }

                    for (; i < n; i++) p[i] = p[i] / d;
                    return true;
                }
                case NPTypeCode.Single:
                {
                    float* p = (float*)block;
                    float d = count; // the boxed loop's `float d = count` — the same (rounding) conversion
                    long i = 0;
                    if (Vector.IsHardwareAccelerated)
                    {
                        var vd = new Vector<float>(d);
                        int w = Vector<float>.Count;
                        for (; i <= n - w; i += w)
                            Unsafe.WriteUnaligned(p + i, Unsafe.ReadUnaligned<Vector<float>>(p + i) / vd);
                    }

                    for (; i < n; i++) p[i] = p[i] / d;
                    return true;
                }
                case NPTypeCode.Complex:
                {
                    // np.mean's true_divide(sum, count) runs NumPy's COMPLEX division loop with the divisor n+0j
                    // (Smith's algorithm: rat = 0, scl = 1/n, ((re + im*0)*scl, (im - re*0)*scl)): a component-wise
                    // re / n rounds differently (re * (1/n) is not re / n for most n) and keeps an infinite or NaN part
                    // from poisoning the other one (np.mean([5+infj]) is nan+infj).
                    var p = (System.Numerics.Complex*)block;
                    for (long i = 0; i < n; i++)
                        p[i] = DefaultEngine.ComplexDivideByCountLikeNumPy(p[i], count);

                    return true;
                }
                case NPTypeCode.Decimal:
                {
                    var p = (decimal*)block;
                    decimal d = count;
                    for (long i = 0; i < n; i++) p[i] = p[i] / d;
                    return true;
                }
                default:
                    return false;
            }
        }

        // =====================================================================
        // Complex kernels
        // =====================================================================

        private const long ComplexBytes = 16; // sizeof(System.Numerics.Complex)

        // out += sum(in)  — double-pair SIMD for the contiguous fast path.
        [MethodImpl(OptimizeAndInline)]
        private static unsafe void ComplexSumKernel(void** dataptrs, long* strides, long count, void* auxdata)
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];

            if (outS == 0)
            {
                // PINNED: fold the whole stripe into one (re,im) slot.
                double* o = (double*)outp;
                double re = o[0], im = o[1];

                if (inS == ComplexBytes)
                {
                    double* d = (double*)inp;
                    long n2 = count * 2;
                    long i = 0;
                    if (Vector256.IsHardwareAccelerated && n2 >= 4)
                    {
                        var acc = Vector256<double>.Zero;
                        for (; i + 4 <= n2; i += 4)
                            acc = Vector256.Add(acc, Vector256.Load(d + i));
                        double* tmp = stackalloc double[4];
                        Vector256.Store(acc, tmp);
                        re += tmp[0] + tmp[2];
                        im += tmp[1] + tmp[3];
                    }
                    for (; i < n2; i += 2) { re += d[i]; im += d[i + 1]; }
                }
                else
                {
                    for (long k = 0; k < count; k++)
                    {
                        double* c = (double*)(inp + k * inS);
                        re += c[0]; im += c[1];
                    }
                }

                o[0] = re; o[1] = im;
                return;
            }

            // SLAB: out[c] += in[c].
            if (inS == ComplexBytes && outS == ComplexBytes)
            {
                double* id = (double*)inp;
                double* od = (double*)outp;
                long n2 = count * 2;
                long i = 0;
                if (Vector256.IsHardwareAccelerated)
                    for (; i + 4 <= n2; i += 4)
                        Vector256.Store(Vector256.Add(Vector256.Load(od + i), Vector256.Load(id + i)), od + i);
                for (; i < n2; i++) od[i] += id[i];
            }
            else
            {
                for (long k = 0; k < count; k++)
                {
                    double* c = (double*)(inp + k * inS);
                    double* o = (double*)(outp + k * outS);
                    o[0] += c[0]; o[1] += c[1];
                }
            }
        }

        // out *= prod(in)  — scalar complex multiply (no SIMD form for complex mul).
        // Inlined naive (ac-bd, ad+bc) on raw doubles: System.Numerics.Complex's
        // operator* carries .NET's infinity-rescue branch (Smith-style) and doesn't
        // inline through `*=` here, which made prod ~5 ns/elem. The naive formula is
        // also what NumPy uses, so this is closer parity on inf/0 mixes too.
        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static unsafe void ComplexProdKernel(void** dataptrs, long* strides, long count, void* auxdata)
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];

            if (outS == 0)
            {
                double* o = (double*)outp;
                double pr = o[0], pi = o[1];
                for (long k = 0; k < count; k++)
                {
                    double* c = (double*)(inp + k * inS);
                    double br = c[0], bi = c[1];
                    double nr = pr * br - pi * bi;
                    pi = pr * bi + pi * br;
                    pr = nr;
                }
                o[0] = pr; o[1] = pi;
            }
            else
            {
                for (long k = 0; k < count; k++)
                {
                    double* c = (double*)(inp + k * inS);
                    double* o = (double*)(outp + k * outS);
                    double ar = o[0], ai = o[1], br = c[0], bi = c[1];
                    double nr = ar * br - ai * bi;
                    o[1] = ar * bi + ai * br;
                    o[0] = nr;
                }
            }
        }

        private static unsafe void ComplexMinKernel(void** dataptrs, long* strides, long count, void* auxdata)
            => ComplexMinMax(dataptrs, strides, count, pickGreater: false);

        private static unsafe void ComplexMaxKernel(void** dataptrs, long* strides, long count, void* auxdata)
            => ComplexMinMax(dataptrs, strides, count, pickGreater: true);

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static unsafe void ComplexMinMax(void** dataptrs, long* strides, long count, bool pickGreater)
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];

            if (outS == 0)
            {
                double* o = (double*)outp;
                double ar = o[0], ai = o[1];
                for (long k = 0; k < count; k++)
                {
                    double* c = (double*)(inp + k * inS);
                    LexFold(ref ar, ref ai, c[0], c[1], pickGreater);
                }
                o[0] = ar; o[1] = ai;
            }
            else
            {
                for (long k = 0; k < count; k++)
                {
                    double* c = (double*)(inp + k * inS);
                    double* o = (double*)(outp + k * outS);
                    double ar = o[0], ai = o[1];
                    LexFold(ref ar, ref ai, c[0], c[1], pickGreater);
                    o[0] = ar; o[1] = ai;
                }
            }
        }

        /// <summary>
        /// NumPy-parity Complex Min/Max fold on raw (re,im) doubles: a NaN-containing
        /// accumulator wins (stays), otherwise the incoming value replaces it when it
        /// is lexicographically more extreme on (Real, Imaginary). Same semantics as
        /// DirectILKernelGenerator.ComplexLexPick (the legacy path) — NaN tested via
        /// the branch-light <c>x != x</c> idiom.
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static void LexFold(ref double ar, ref double ai, double br, double bi, bool pickGreater)
        {
            if (ar != ar || ai != ai) return;           // accumulator is NaN → keep it
            if (br != br || bi != bi) { ar = br; ai = bi; return; } // incoming NaN → take it
            bool aGreater = ar > br || (ar == br && ai > bi);
            bool takeA = pickGreater ? aGreater : !aGreater;
            if (!takeA) { ar = br; ai = bi; }
        }

        // =====================================================================
        // Generic scalar kernels (Half / Decimal — and any INumber numeric)
        // =====================================================================

        /// <summary>
        ///     Reads one input element and converts it to the accumulator type the way NumPy's loop casts its inputs: a
        ///     plain read when the types match, NumPy's C cast for a float / complex / decimal element into an integer
        ///     accumulator, a widening <c>CreateTruncating</c> otherwise.
        /// </summary>
        /// <typeparam name="TIn">The input element type.</typeparam>
        /// <typeparam name="TAccum">The accumulator (loop) type.</typeparam>
        /// <param name="p">The element's address.</param>
        /// <returns>The element as <typeparamref name="TAccum"/>.</returns>
        /// <remarks>
        ///     <para>
        ///     Every <c>typeof</c> test is a JIT constant per instantiation, so each (TIn, TAccum) pair compiles to one of
        ///     the three branches: the same-type read stays a reinterpret (the difference between Half pinned-reduce at
        ///     ~parity vs +30%), and integer widening stays <c>CreateTruncating</c>.
        ///     </para>
        ///     <para>
        ///     A floating source into an integer accumulator is the cumsum / nancumsum <c>dtype=int</c> case. NumPy's
        ///     iterator casts each element with C's conversion (truncate toward zero, then wrap; NaN / inf / out of range
        ///     give the integer minimum truncated to the width, 0 for 8- and 16-bit), the cast <c>astype</c> runs, while
        ///     .NET's <c>CreateTruncating</c> SATURATES: <c>np.cumsum([-3.0, 1.0], axis=0, dtype=uint8)</c> gave
        ///     <c>[0, 1]</c> for NumPy's <c>[253, 254]</c> (parity QA-06). <see cref="Converts.ChangeType{TIn, TOut}"/> is
        ///     that NumPy cast (typed for float32 / float64 / decimal, boxed for float16 / complex128).
        ///     </para>
        /// </remarks>
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static unsafe TAccum ConvIn<TIn, TAccum>(byte* p)
            where TIn : unmanaged, INumberBase<TIn>
            where TAccum : unmanaged, INumber<TAccum>
        {
            if (typeof(TIn) == typeof(TAccum))
                return *(TAccum*)p;
            if (IsFloatingSource<TIn>() && IsIntegerAccumulator<TAccum>())
                return Converts.ChangeType<TIn, TAccum>(*(TIn*)p);
            return TAccum.CreateTruncating(*(TIn*)p);
        }

        /// <summary>
        ///     Whether <typeparamref name="T"/> is an element type NumPy casts to an integer loop with C's float conversion
        ///     (float16 / float32 / float64, NumSharp's decimal, and complex128, whose real part is cast).
        /// </summary>
        /// <typeparam name="T">The input element type.</typeparam>
        /// <returns>True for the floating element types; a JIT constant per instantiation.</returns>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static bool IsFloatingSource<T>()
            => typeof(T) == typeof(double) || typeof(T) == typeof(float) || typeof(T) == typeof(Half)
               || typeof(T) == typeof(decimal) || typeof(T) == typeof(System.Numerics.Complex);

        /// <summary>
        ///     Whether <typeparamref name="T"/> is an integer accumulator (the signed and unsigned integers and char, the
        ///     uint16 lane), whose float conversion must wrap like C rather than saturate like .NET.
        /// </summary>
        /// <typeparam name="T">The accumulator type.</typeparam>
        /// <returns>True for the integer accumulators; a JIT constant per instantiation.</returns>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        private static bool IsIntegerAccumulator<T>()
            => typeof(T) == typeof(sbyte) || typeof(T) == typeof(byte) || typeof(T) == typeof(short)
               || typeof(T) == typeof(ushort) || typeof(T) == typeof(int) || typeof(T) == typeof(uint)
               || typeof(T) == typeof(long) || typeof(T) == typeof(ulong) || typeof(T) == typeof(char);

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static unsafe void TypedSumKernel<TIn, TAccum>(void** dataptrs, long* strides, long count, void* auxdata)
            where TIn : unmanaged, INumberBase<TIn>
            where TAccum : unmanaged, INumber<TAccum>
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];
            if (outS == 0)
            {
                TAccum acc = *(TAccum*)outp;
                for (long k = 0; k < count; k++) acc += ConvIn<TIn, TAccum>(inp + k * inS);
                *(TAccum*)outp = acc;
            }
            else
            {
                for (long k = 0; k < count; k++)
                {
                    TAccum* o = (TAccum*)(outp + k * outS);
                    *o += ConvIn<TIn, TAccum>(inp + k * inS);
                }
            }
        }

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static unsafe void TypedProdKernel<TIn, TAccum>(void** dataptrs, long* strides, long count, void* auxdata)
            where TIn : unmanaged, INumberBase<TIn>
            where TAccum : unmanaged, INumber<TAccum>
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];
            if (outS == 0)
            {
                TAccum acc = *(TAccum*)outp;
                for (long k = 0; k < count; k++) acc *= ConvIn<TIn, TAccum>(inp + k * inS);
                *(TAccum*)outp = acc;
            }
            else
            {
                for (long k = 0; k < count; k++)
                {
                    TAccum* o = (TAccum*)(outp + k * outS);
                    *o *= ConvIn<TIn, TAccum>(inp + k * inS);
                }
            }
        }

        private static unsafe void TypedMinKernel<TIn, TAccum>(void** dataptrs, long* strides, long count, void* auxdata)
            where TIn : unmanaged, INumberBase<TIn>
            where TAccum : unmanaged, INumber<TAccum>
            => TypedMinMax<TIn, TAccum>(dataptrs, strides, count, pickGreater: false);

        private static unsafe void TypedMaxKernel<TIn, TAccum>(void** dataptrs, long* strides, long count, void* auxdata)
            where TIn : unmanaged, INumberBase<TIn>
            where TAccum : unmanaged, INumber<TAccum>
            => TypedMinMax<TIn, TAccum>(dataptrs, strides, count, pickGreater: true);

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static unsafe void TypedMinMax<TIn, TAccum>(void** dataptrs, long* strides, long count, bool pickGreater)
            where TIn : unmanaged, INumberBase<TIn>
            where TAccum : unmanaged, INumber<TAccum>
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];
            if (outS == 0)
            {
                TAccum acc = *(TAccum*)outp;
                for (long k = 0; k < count; k++)
                    acc = MinMaxFold(acc, ConvIn<TIn, TAccum>(inp + k * inS), pickGreater);
                *(TAccum*)outp = acc;
            }
            else
            {
                for (long k = 0; k < count; k++)
                {
                    TAccum* o = (TAccum*)(outp + k * outS);
                    *o = MinMaxFold(*o, ConvIn<TIn, TAccum>(inp + k * inS), pickGreater);
                }
            }
        }

        /// <summary>
        /// NaN-propagating min/max fold (NumPy parity): a NaN accumulator wins (stays); an
        /// incoming NaN replaces a finite accumulator. For types without NaN (Decimal) the
        /// <c>IsNaN</c> calls fold to false and this is a plain comparison.
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static TAccum MinMaxFold<TAccum>(TAccum acc, TAccum v, bool pickGreater)
            where TAccum : unmanaged, INumber<TAccum>
        {
            if (TAccum.IsNaN(acc)) return acc;
            if (TAccum.IsNaN(v)) return v;
            bool takeV = pickGreater ? (v > acc) : (v < acc);
            return takeV ? v : acc;
        }

        // =====================================================================
        // SIMD same-type Min/Max (Phase 6 — double/float). Dual-mode, NaN-propagating
        // (NumPy parity). Mirrors PairwiseSumSameType so f64/f32 Min/Max ride the
        // stride-ordered NDIter reduce path at SIMD speed on EVERY layout
        // (transpose / F / strided / negcol / broadcast) instead of falling to the
        // C-contiguity-gated Direct axis kernel (which collapses to a cache-hostile
        // coordinate walk on negative-stride / broadcast inputs — measured 6–17×
        // slower than NumPy). The scalar tail / strided fallback reuse MinMaxFold so
        // they are bit-identical to the Half/Decimal/Complex TypedMinMax path.
        //
        // NaN handling: Vector256.Min/Max do NOT propagate NaN (x86 MINPD/MAXPD return
        // the SECOND operand when either input is NaN), so each step blends with
        // ConditionalSelect to keep the actual NaN OPERAND (accumulator first, then
        // input) — reproducing np.minimum/np.maximum's "NaN wins" AND preserving the
        // input NaN's bit pattern/sign exactly like the scalar MinMaxFold (NumPy
        // propagates the input NaN, not .NET's negative double.NaN). The non-NaN result
        // is order-independent (min/max associative+commutative), so the lane-parallel
        // fold matches the sequential scalar fold for all finite/inf inputs.
        // =====================================================================

        private static unsafe void SimdMinKernel<T>(void** dataptrs, long* strides, long count, void* auxdata)
            where T : unmanaged, IFloatingPointIeee754<T>
            => SimdMinMaxSameType<T>(dataptrs, strides, count, pickGreater: false);

        private static unsafe void SimdMaxKernel<T>(void** dataptrs, long* strides, long count, void* auxdata)
            where T : unmanaged, IFloatingPointIeee754<T>
            => SimdMinMaxSameType<T>(dataptrs, strides, count, pickGreater: true);

        // Per-lane NaN-propagating min/max: keep acc's NaN if acc is NaN, else input's
        // NaN if input is NaN, else the plain Vector256 min/max. Mirrors MinMaxFold lane-wise
        // (and therefore preserves the input NaN bits NumPy propagates).
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static Vector256<T> NaNFoldVec<T>(Vector256<T> acc, Vector256<T> v, bool pickGreater)
            where T : unmanaged, IFloatingPointIeee754<T>
        {
            // m is consumed only in lanes where NEITHER operand is NaN (the ConditionalSelect
            // chain overrides NaN lanes), so the raw VMINPD/VMAXPD (no NaN fixup) is exact here.
            var m = pickGreater ? RawMax256(acc, v) : RawMin256(acc, v);
            var accNaN = ~Vector256.Equals(acc, acc);   // acc != acc ⇒ NaN lane
            var vNaN = ~Vector256.Equals(v, v);
            return Vector256.ConditionalSelect(accNaN, acc, Vector256.ConditionalSelect(vNaN, v, m));
        }

        // Raw x86 Avx.Min/Max for float/double — a single VMINPD/VMAXPD WITHOUT the net9+ JIT
        // NaN-propagation fixup (an extra compare+blend) that Vector256.Min/Max carry (~2x the
        // raw instruction). SAFE in this file because every caller either tracks a separate
        // finite mask + cold NaN scan (SimdMinMaxSameType) or overrides NaN lanes via
        // ConditionalSelect (NaNFoldVec), so the raw op's NaN-dropping is intended. Non-x86
        // (ARM) falls back to Vector256.Min/Max (no x86 fixup there). typeof(T)/IsSupported
        // are JIT constants, so each specialization compiles to the single intrinsic.
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static Vector256<T> RawMin256<T>(Vector256<T> a, Vector256<T> b) where T : unmanaged
        {
            if (Avx.IsSupported)
            {
                if (typeof(T) == typeof(double)) return Avx.Min(a.AsDouble(), b.AsDouble()).As<double, T>();
                if (typeof(T) == typeof(float)) return Avx.Min(a.AsSingle(), b.AsSingle()).As<float, T>();
            }
            return Vector256.Min(a, b);
        }
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static Vector256<T> RawMax256<T>(Vector256<T> a, Vector256<T> b) where T : unmanaged
        {
            if (Avx.IsSupported)
            {
                if (typeof(T) == typeof(double)) return Avx.Max(a.AsDouble(), b.AsDouble()).As<double, T>();
                if (typeof(T) == typeof(float)) return Avx.Max(a.AsSingle(), b.AsSingle()).As<float, T>();
            }
            return Vector256.Max(a, b);
        }

        [MethodImpl(MethodImplOptions.AggressiveOptimization)]
        private static unsafe void SimdMinMaxSameType<T>(void** dataptrs, long* strides, long count, bool pickGreater)
            where T : unmanaged, IFloatingPointIeee754<T>
        {
            byte* inp = (byte*)dataptrs[0]; long inS = strides[0];
            byte* outp = (byte*)dataptrs[1]; long outS = strides[1];
            int sz = sizeof(T);
            int W = Vector256<T>.Count;

            if (outS == 0)
            {
                // PINNED: fold the contiguous stripe into the single slot. The slot is
                // seeded ±inf (SeedReduceIdentity) and accumulates across buffered chunks.
                T cur = *(T*)outp;
                if (T.IsNaN(cur)) return;            // NaN already won — stays NaN
                T* id = (T*)inp;

                if (inS == sz && Vector256.IsHardwareAccelerated && count >= W)
                {
                    // Hot loop: 4×-unrolled PLAIN min/max (1 op/lane, like the sum SLAB) plus a
                    // cheap per-lane finite mask. Vector256.Min/Max don't propagate NaN, but a
                    // NaN in a min/max stripe is rare; the finite mask flags it and a cold scalar
                    // scan then recomputes the exact result (the first NaN's bits — NumPy/scalar
                    // parity). This keeps the common (no-NaN) path at sum-class throughput; the
                    // earlier NaN-correct-every-step fold (NaNFoldVec) was ~5 ops/lane → 2–3×
                    // slower than the Direct whole-array kernel on contiguous reduced axes.
                    var a0 = Vector256.Create(cur); var a1 = a0; var a2 = a0; var a3 = a0;
                    var fin = Vector256<T>.AllBitsSet;
                    long i = 0, lim = count - count % (W * 4);
                    for (; i < lim; i += W * 4)
                    {
                        var v0 = Vector256.Load(id + i);
                        var v1 = Vector256.Load(id + i + W);
                        var v2 = Vector256.Load(id + i + 2 * W);
                        var v3 = Vector256.Load(id + i + 3 * W);
                        if (pickGreater)
                        { a0 = RawMax256(a0, v0); a1 = RawMax256(a1, v1); a2 = RawMax256(a2, v2); a3 = RawMax256(a3, v3); }
                        else
                        { a0 = RawMin256(a0, v0); a1 = RawMin256(a1, v1); a2 = RawMin256(a2, v2); a3 = RawMin256(a3, v3); }
                        fin &= Vector256.Equals(v0, v0) & Vector256.Equals(v1, v1) & Vector256.Equals(v2, v2) & Vector256.Equals(v3, v3);
                    }
                    var va = pickGreater
                        ? RawMax256(RawMax256(a0, a1), RawMax256(a2, a3))
                        : RawMin256(RawMin256(a0, a1), RawMin256(a2, a3));
                    for (; i + W <= count; i += W)
                    {
                        var v = Vector256.Load(id + i);
                        va = pickGreater ? RawMax256(va, v) : RawMin256(va, v);
                        fin &= Vector256.Equals(v, v);
                    }
                    // fin lane = all-ones ⇒ that lane was finite throughout; MSB-extract → bit set.
                    bool anyNaN = Vector256.ExtractMostSignificantBits(fin) != (uint)((1 << W) - 1);
                    T acc = va.GetElement(0);
                    for (int l = 1; l < W; l++) acc = MinMaxFold(acc, va.GetElement(l), pickGreater);
                    for (; i < count; i++)
                    {
                        T x = id[i];
                        if (T.IsNaN(x)) anyNaN = true;
                        acc = MinMaxFold(acc, x, pickGreater);
                    }
                    // Cold path: a NaN was present ⇒ the reduction is NaN. Propagate the FIRST
                    // input NaN's exact bits, matching the sequential scalar fold and np.minimum.
                    if (anyNaN)
                        for (long k = 0; k < count; k++) { T x = id[k]; if (T.IsNaN(x)) { acc = x; break; } }
                    *(T*)outp = acc;
                }
                else
                {
                    T acc = cur;
                    for (long k = 0; k < count; k++)
                        acc = MinMaxFold(acc, *(T*)(inp + k * inS), pickGreater);
                    *(T*)outp = acc;
                }
            }
            else
            {
                // SLAB: out[c] = minmax(out[c], in[c]); the slab is revisited across outer
                // iterations (seeded ±inf), so fold in place.
                if (inS == sz && outS == sz && Vector256.IsHardwareAccelerated)
                {
                    T* id = (T*)inp; T* od = (T*)outp;
                    long i = 0, lim = count - count % W;
                    for (; i < lim; i += W)
                        Vector256.Store(NaNFoldVec(Vector256.Load(od + i), Vector256.Load(id + i), pickGreater), od + i);
                    for (; i < count; i++) od[i] = MinMaxFold(od[i], id[i], pickGreater);
                }
                else
                {
                    for (long k = 0; k < count; k++)
                    {
                        T* o = (T*)(outp + k * outS);
                        *o = MinMaxFold(*o, *(T*)(inp + k * inS), pickGreater);
                    }
                }
            }
        }
    }
}