DirectILKernelGenerator.cs

using System;
using System.Linq;
using System.Reflection;
using System.Reflection.Emit;
using System.Runtime.CompilerServices;
using System.Runtime.Intrinsics;
using NumSharp.Utilities;

// =============================================================================
// DirectILKernelGenerator - IL-based SIMD kernel generation using DynamicMethod
// =============================================================================
//
// ARCHITECTURE OVERVIEW
// ---------------------
// This class generates high-performance kernels at runtime using IL emission.
// The JIT compiler can then optimize these kernels with full SIMD support (V128/V256/V512).
// Kernels are cached by operation key to avoid repeated IL generation.
//
// FLOW: Caller (DefaultEngine, np.*, NDArray ops)
//         -> Requests kernel via Get*Kernel() or *Helper() methods
//         -> DirectILKernelGenerator checks cache, generates IL if needed
//         -> Returns delegate that caller invokes with array pointers
//
// DESIGN: Static class - all kernel methods are static.
// Call them directly from DefaultEngine.
//
// EXCEPTION HANDLING
// ------------------
// All TryGet*Kernel() methods use catch-all exception handling that returns null.
// This is intentional graceful degradation: if IL generation fails for any reason
// (unsupported type, reflection error, invalid IL sequence), the caller receives
// null and falls back to an alternative code path (typically scalar loops or
// throwing NotSupportedException with a descriptive message).
//
// This pattern exists in 14 locations across the partial class files:
//   - Binary.cs: TryGenerateContiguousKernel
//   - MixedType.cs: TryGetMixedTypeKernel
//   - MatMul.cs: GenerateMatMulKernelIL
//   - Unary.cs: TryGetUnaryKernel
//   - Shift.cs: GetShiftScalarKernel, GetShiftArrayKernel
//   - Scan.cs: TryGetCumulativeKernel, TryGetCumulativeAxisKernel
//   - Reduction.cs: TryGetTypedElementReductionKernel
//   - Comparison.cs: TryGetComparisonKernel
//   - DirectILKernelGenerator.cs: Core kernel infrastructure
//
// =============================================================================
// PARTIAL CLASS FILES
// =============================================================================
//
// DirectILKernelGenerator.cs (THIS FILE)
//   OWNERSHIP: Core infrastructure - foundation for all other partial files
//   RESPONSIBILITY:
//     - Static kernel generation methods used by DefaultEngine
//     - Global state: Enabled flag, VectorBits/VectorBytes (detected at startup)
//     - Type mapping: NPTypeCode <-> CLR Type <-> Vector type conversions
//     - Shared IL emission primitives used by all other partials
//   DEPENDENCIES: None (other partials depend on this)
//   KEY MEMBERS:
//     - Enabled, VectorBits, VectorBytes - runtime SIMD capability
//     - GetVectorContainerType(), GetVectorType() - V128/V256/V512 type selection
//     - GetTypeSize(), GetClrType(), CanUseSimd(), IsUnsigned() - type utilities
//     - EmitLoadIndirect(), EmitStoreIndirect() - memory access IL
//     - EmitConvertTo(), EmitScalarOperation() - type conversion and scalar ops
//     - EmitVectorLoad/Store/Create/Operation() - SIMD operations
//
// DirectILKernelGenerator.Binary.cs
//   OWNERSHIP: Same-type binary operations on contiguous arrays (fast path)
//   RESPONSIBILITY:
//     - Optimized kernels when both operands have identical type and layout
//     - SIMD loop + scalar tail for Add, Sub, Mul, Div
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by DefaultEngine for same-type contiguous operations
//   KEY MEMBERS:
//     - ContiguousKernel<T> delegate, _contiguousKernelCache
//     - GetContiguousKernel<T>()
//     - Generic helpers: IsSimdSupported<T>(), EmitLoadIndirect<T>(), etc.
//
// DirectILKernelGenerator.MixedType.cs
//   OWNERSHIP: Mixed-type binary operations with type promotion
//   RESPONSIBILITY:
//     - Handles all binary ops where operand types may differ
//     - Generates path-specific kernels based on stride patterns
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by DefaultEngine for general binary operations
//   KEY MEMBERS:
//     - MixedTypeKernel delegate, _mixedTypeCache
//     - GetMixedTypeKernel(), TryGetMixedTypeKernel()
//     - Path generators: GenerateSimdFullKernel(), GenerateGeneralKernel(), etc.
//     - Loop emitters: EmitScalarFullLoop(), EmitSimdFullLoop(), EmitGeneralLoop()
//
// DirectILKernelGenerator.Unary.cs
//   OWNERSHIP: Unary element-wise operations
//   RESPONSIBILITY:
//     - Math functions: Negate, Abs, Sqrt, Sin, Cos, Exp, Log, Sign, Floor, Ceil, etc.
//     - Scalar delegate generation for single-value operations (Func<TIn,TOut>)
//     - Binary scalar delegates for mixed-type scalar operations
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by DefaultEngine for unary ops; scalar delegates used in broadcasting
//   KEY MEMBERS:
//     - UnaryKernel delegate, _unaryCache, _unaryScalarCache, _binaryScalarCache
//     - GetUnaryKernel(), GetUnaryScalarDelegate(), GetBinaryScalarDelegate()
//     - EmitUnaryScalarOperation(), EmitMathCall(), EmitSignCall()
//
// DirectILKernelGenerator.Comparison.cs
//   OWNERSHIP: Comparison operations returning boolean arrays
//   RESPONSIBILITY:
//     - Element-wise comparisons: ==, !=, <, >, <=, >=
//     - SIMD comparison with efficient mask-to-bool extraction
//     - Scalar comparison delegates for single-value operations
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by NDArray comparison operators (==, !=, <, >, etc.)
//   KEY MEMBERS:
//     - ComparisonKernel delegate, _comparisonCache, _comparisonScalarCache
//     - GetComparisonKernel(), GetComparisonScalarDelegate()
//     - EmitVectorComparison(), EmitMaskToBoolExtraction()
//
// DirectILKernelGenerator.Reduction.cs
//   OWNERSHIP: Reduction operations and specialized SIMD helpers
//   RESPONSIBILITY:
//     - Reductions: Sum, Prod, Min, Max, Mean, ArgMax, ArgMin, All, Any
//     - SIMD helpers called DIRECTLY by other NumSharp code (not just via kernels):
//       * All/Any with early-exit optimization
//       * ArgMax/ArgMin with SIMD two-pass (find value, then find index)
//       * Boolean masking: CountTrue, CopyMaskedElements
//   (np.nonzero / np.argwhere live in DirectILKernelGenerator.Argwhere.cs +
//    DirectILKernelGenerator.NonZero.cs as per-dtype IL kernels.)
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Kernels called by DefaultEngine; helpers called directly by np.all/any/nonzero/masking
//   KEY MEMBERS:
//     - TypedElementReductionKernel<T> delegate, _elementReductionCache
//     - GetTypedElementReductionKernel<T>()
//     - AllSimdHelper<T>(), AnySimdHelper<T>() - early-exit boolean reductions
//     - ArgMaxSimdHelper<T>(), ArgMinSimdHelper<T>() - index-tracking reductions
//     - CountTrueSimdHelper(), CopyMaskedElementsHelper<T>()
//     - EmitTreeReduction(), EmitVectorHorizontalReduction()
//
// =============================================================================

// =============================================================================
// DirectILKernelGenerator - IL-based SIMD kernel generation using DynamicMethod
// =============================================================================
//
// ARCHITECTURE OVERVIEW
// ---------------------
// This partial class generates high-performance kernels at runtime using IL emission.
// The JIT compiler can then optimize these kernels with full SIMD support (V128/V256/V512).
// Kernels are cached by operation key to avoid repeated IL generation.
//
// FLOW: Caller (DefaultEngine, np.*, NDArray ops)
//         -> Requests kernel via Get*Kernel() or *Helper() methods
//         -> DirectILKernelGenerator checks cache, generates IL if needed
//         -> Returns delegate that caller invokes with array pointers
//
// =============================================================================
// PARTIAL CLASS FILES
// =============================================================================
//
// DirectILKernelGenerator.cs (THIS FILE)
//   OWNERSHIP: Core infrastructure - foundation for all other partial files
//   RESPONSIBILITY:
//     - Global state: Enabled flag, VectorBits/VectorBytes (detected at startup)
//     - Type mapping: NPTypeCode <-> CLR Type <-> Vector type conversions
//     - Shared IL emission primitives used by all other partials
//   DEPENDENCIES: None (other partials depend on this)
//   KEY MEMBERS:
//     - Enabled, VectorBits, VectorBytes - runtime SIMD capability
//     - GetVectorContainerType(), GetVectorType() - V128/V256/V512 type selection
//     - GetTypeSize(), GetClrType(), CanUseSimd(), IsUnsigned() - type utilities
//     - EmitLoadIndirect(), EmitStoreIndirect() - memory access IL
//     - EmitConvertTo(), EmitScalarOperation() - type conversion and scalar ops
//     - EmitVectorLoad/Store/Create/Operation() - SIMD operations
//
// DirectILKernelGenerator.Binary.cs
//   OWNERSHIP: Same-type binary operations on contiguous arrays (fast path)
//   RESPONSIBILITY:
//     - Optimized kernels when both operands have identical type and layout
//     - SIMD loop + scalar tail for Add, Sub, Mul, Div
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by DefaultEngine for same-type contiguous operations
//   KEY MEMBERS:
//     - ContiguousKernel<T> delegate, _contiguousKernelCache
//     - GetContiguousKernel<T>(), GenerateUnifiedKernel<T>()
//     - Generic helpers: IsSimdSupported<T>(), EmitLoadIndirect<T>(), etc.
//
// DirectILKernelGenerator.MixedType.cs
//   OWNERSHIP: Mixed-type binary operations with type promotion
//   RESPONSIBILITY:
//     - Handles all binary ops where operand types may differ
//     - Generates path-specific kernels based on stride patterns
//     - Owns ClearAll() which clears ALL caches across all partials
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by DefaultEngine for general binary operations
//   KEY MEMBERS:
//     - MixedTypeKernel delegate, _mixedTypeCache
//     - GetMixedTypeKernel(), TryGetMixedTypeKernel(), ClearAll()
//     - Path generators: GenerateSimdFullKernel(), GenerateGeneralKernel(), etc.
//     - Loop emitters: EmitScalarFullLoop(), EmitSimdFullLoop(), EmitGeneralLoop()
//
// DirectILKernelGenerator.Unary.cs
//   OWNERSHIP: Unary element-wise operations
//   RESPONSIBILITY:
//     - Math functions: Negate, Abs, Sqrt, Sin, Cos, Exp, Log, Sign, Floor, Ceil, etc.
//     - Scalar delegate generation for single-value operations (Func<TIn,TOut>)
//     - Binary scalar delegates for mixed-type scalar operations
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by DefaultEngine for unary ops; scalar delegates used in broadcasting
//   KEY MEMBERS:
//     - UnaryKernel delegate, _unaryCache, _unaryScalarCache, _binaryScalarCache
//     - GetUnaryKernel(), GetUnaryScalarDelegate(), GetBinaryScalarDelegate()
//     - EmitUnaryScalarOperation(), EmitMathCall(), EmitSignCall()
//
// DirectILKernelGenerator.Comparison.cs
//   OWNERSHIP: Comparison operations returning boolean arrays
//   RESPONSIBILITY:
//     - Element-wise comparisons: ==, !=, <, >, <=, >=
//     - SIMD comparison with efficient mask-to-bool extraction
//     - Scalar comparison delegates for single-value operations
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Called by NDArray comparison operators (==, !=, <, >, etc.)
//   KEY MEMBERS:
//     - ComparisonKernel delegate, _comparisonCache, _comparisonScalarCache
//     - GetComparisonKernel(), GetComparisonScalarDelegate()
//     - EmitVectorComparison(), EmitMaskToBoolExtraction()
//
// DirectILKernelGenerator.Reduction.cs
//   OWNERSHIP: Reduction operations and specialized SIMD helpers
//   RESPONSIBILITY:
//     - Reductions: Sum, Prod, Min, Max, Mean, ArgMax, ArgMin, All, Any
//     - SIMD helpers called DIRECTLY by other NumSharp code (not just via kernels):
//       * All/Any with early-exit optimization
//       * ArgMax/ArgMin with SIMD two-pass (find value, then find index)
//       * Boolean masking: CountTrue, CopyMaskedElements
//   (np.nonzero / np.argwhere live in DirectILKernelGenerator.Argwhere.cs +
//    DirectILKernelGenerator.NonZero.cs as per-dtype IL kernels.)
//   DEPENDENCIES: Uses core emit helpers from DirectILKernelGenerator.cs
//   FLOW: Kernels called by DefaultEngine; helpers called directly by np.all/any/nonzero/masking
//   KEY MEMBERS:
//     - TypedElementReductionKernel<T> delegate, _elementReductionCache
//     - GetTypedElementReductionKernel<T>(), ClearReduction()
//     - AllSimdHelper<T>(), AnySimdHelper<T>() - early-exit boolean reductions
//     - ArgMaxSimdHelper<T>(), ArgMinSimdHelper<T>() - index-tracking reductions
//     - CountTrueSimdHelper(), CopyMaskedElementsHelper<T>()
//     - EmitTreeReduction(), EmitVectorHorizontalReduction()
//
// =============================================================================

namespace NumSharp.Backends.Kernels
{
    /// <summary>
    /// Generates IL-based SIMD kernels using DynamicMethod.
    /// These kernels provide ~10-15% speedup over the C# reference implementations
    /// by allowing the JIT to inline Vector256 operations more aggressively.
    ///
    /// This class is internal to NumSharp.Backends - all kernel access should go
    /// through TensorEngine/DefaultEngine, not directly from np.* or NDArray.
    /// </summary>
    public static partial class DirectILKernelGenerator
    {
        #region Static Configuration

        /// <summary>
        /// Provider name for diagnostics.
        /// </summary>
        public static string Name => "IL";

        /// <summary>
        /// Whether IL generation is enabled. Can be disabled for debugging.
        /// </summary>
        public static bool Enabled { get; set; } = true;

        /// <summary>
        /// Detected vector width at startup: 512, 256, 128, or 0 (no SIMD).
        /// </summary>
        public static readonly int VectorBits =
            Vector512.IsHardwareAccelerated ? 512 :
            Vector256.IsHardwareAccelerated ? 256 :
            Vector128.IsHardwareAccelerated ? 128 : 0;

        /// <summary>
        /// Number of bytes per vector register.
        /// </summary>
        public static readonly int VectorBytes = VectorBits / 8;

        #endregion

        #region Cached MethodInfo Lookups

        /// <summary>
        /// Pre-cached MethodInfo references for frequently used reflection calls.
        /// Caching these avoids repeated GetMethod() lookups during kernel generation.
        /// All fields use ?? throw to fail fast at type load if a method is not found.
        /// </summary>
        private static partial class CachedMethods
        {
            // Math methods (double versions)
            public static readonly MethodInfo MathPow = typeof(Math).GetMethod(nameof(Math.Pow), new[] { typeof(double), typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, nameof(Math.Pow));
            public static readonly MethodInfo MathFPow = typeof(MathF).GetMethod(nameof(MathF.Pow), new[] { typeof(float), typeof(float) })
                ?? throw new MissingMethodException(typeof(MathF).FullName, nameof(MathF.Pow));
            public static readonly MethodInfo MathFloor = typeof(Math).GetMethod(nameof(Math.Floor), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, nameof(Math.Floor));
            public static readonly MethodInfo MathAtan2 = typeof(Math).GetMethod(nameof(Math.Atan2), new[] { typeof(double), typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, nameof(Math.Atan2));

            // np.logaddexp / np.logaddexp2 / np.nextafter / np.copysign scalar kernels — one static
            // helper per (op, loop dtype), resolved through GetLogAddNextMethod. Operands reach
            // EmitScalarOperation already in the loop dtype, so the kernel is a single Call to the
            // matching signature (double/float/Half/decimal). Same shape as the NDDiv helpers below.
            private static MethodInfo LA(string name, Type t) =>
                typeof(Utilities.NDLogAddExpMath).GetMethod(name, new[] { t, t })
                ?? throw new MissingMethodException(typeof(Utilities.NDLogAddExpMath).FullName, name);
            public static readonly MethodInfo LogAddExpD = LA(nameof(Utilities.NDLogAddExpMath.LogAddExp), typeof(double));
            public static readonly MethodInfo LogAddExpF = LA(nameof(Utilities.NDLogAddExpMath.LogAddExpF), typeof(float));
            public static readonly MethodInfo LogAddExpH = LA(nameof(Utilities.NDLogAddExpMath.LogAddExpHalf), typeof(Half));
            public static readonly MethodInfo LogAddExpDec = LA(nameof(Utilities.NDLogAddExpMath.LogAddExpDecimal), typeof(decimal));
            public static readonly MethodInfo LogAddExp2D = LA(nameof(Utilities.NDLogAddExpMath.LogAddExp2), typeof(double));
            public static readonly MethodInfo LogAddExp2F = LA(nameof(Utilities.NDLogAddExpMath.LogAddExp2F), typeof(float));
            public static readonly MethodInfo LogAddExp2H = LA(nameof(Utilities.NDLogAddExpMath.LogAddExp2Half), typeof(Half));
            public static readonly MethodInfo LogAddExp2Dec = LA(nameof(Utilities.NDLogAddExpMath.LogAddExp2Decimal), typeof(decimal));
            public static readonly MethodInfo NextAfterD = LA(nameof(Utilities.NDLogAddExpMath.NextAfter), typeof(double));
            public static readonly MethodInfo NextAfterF = LA(nameof(Utilities.NDLogAddExpMath.NextAfterF), typeof(float));
            public static readonly MethodInfo NextAfterH = LA(nameof(Utilities.NDLogAddExpMath.NextAfterHalf), typeof(Half));
            public static readonly MethodInfo NextAfterDec = LA(nameof(Utilities.NDLogAddExpMath.NextAfterDecimal), typeof(decimal));
            public static readonly MethodInfo CopySignD = LA(nameof(Utilities.NDLogAddExpMath.CopySign), typeof(double));
            public static readonly MethodInfo CopySignF = LA(nameof(Utilities.NDLogAddExpMath.CopySignF), typeof(float));
            public static readonly MethodInfo CopySignH = LA(nameof(Utilities.NDLogAddExpMath.CopySignHalf), typeof(Half));
            public static readonly MethodInfo CopySignDec = LA(nameof(Utilities.NDLogAddExpMath.CopySignDecimal), typeof(decimal));

            // Integer power helpers (squared-exponentiation with native wrapping).
            // Used by EmitPowerOperation when result type is integer to preserve
            // NumPy's exact-wrap semantics that Math.Pow's double round-trip loses.
            public static readonly MethodInfo IntPowSByte = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowSByte), new[] { typeof(sbyte), typeof(sbyte) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowSByte));
            public static readonly MethodInfo IntPowByte = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowByte), new[] { typeof(byte), typeof(byte) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowByte));
            public static readonly MethodInfo IntPowInt16 = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowInt16), new[] { typeof(short), typeof(short) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowInt16));
            public static readonly MethodInfo IntPowUInt16 = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowUInt16), new[] { typeof(ushort), typeof(ushort) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowUInt16));
            public static readonly MethodInfo IntPowChar = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowChar), new[] { typeof(char), typeof(char) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowChar));
            public static readonly MethodInfo IntPowInt32 = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowInt32), new[] { typeof(int), typeof(int) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowInt32));
            public static readonly MethodInfo IntPowUInt32 = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowUInt32), new[] { typeof(uint), typeof(uint) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowUInt32));
            public static readonly MethodInfo IntPowInt64 = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowInt64), new[] { typeof(long), typeof(long) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowInt64));
            public static readonly MethodInfo IntPowUInt64 = typeof(Utilities.NDIntegerPower).GetMethod(
                nameof(Utilities.NDIntegerPower.PowUInt64), new[] { typeof(ulong), typeof(ulong) })
                ?? throw new MissingMethodException(typeof(Utilities.NDIntegerPower).FullName, nameof(Utilities.NDIntegerPower.PowUInt64));

            // floor-division / remainder helpers (NDDivision) — NumPy-exact divide-by-zero (-> 0
            // for integers, ±inf/nan for floats) and floored-sign semantics. Used by
            // EmitFloorDivideOperation / EmitModOperation in place of the old double round-trip.
            private static MethodInfo NDDiv(string name, Type t) =>
                typeof(Utilities.NDDivision).GetMethod(name, new[] { t, t })
                ?? throw new MissingMethodException(typeof(Utilities.NDDivision).FullName, name);

            public static readonly MethodInfo FloorDivSByte = NDDiv(nameof(Utilities.NDDivision.FloorDivSByte), typeof(sbyte));
            public static readonly MethodInfo FloorDivByte = NDDiv(nameof(Utilities.NDDivision.FloorDivByte), typeof(byte));
            public static readonly MethodInfo FloorDivInt16 = NDDiv(nameof(Utilities.NDDivision.FloorDivInt16), typeof(short));
            public static readonly MethodInfo FloorDivUInt16 = NDDiv(nameof(Utilities.NDDivision.FloorDivUInt16), typeof(ushort));
            public static readonly MethodInfo FloorDivChar = NDDiv(nameof(Utilities.NDDivision.FloorDivChar), typeof(char));
            public static readonly MethodInfo FloorDivInt32 = NDDiv(nameof(Utilities.NDDivision.FloorDivInt32), typeof(int));
            public static readonly MethodInfo FloorDivUInt32 = NDDiv(nameof(Utilities.NDDivision.FloorDivUInt32), typeof(uint));
            public static readonly MethodInfo FloorDivInt64 = NDDiv(nameof(Utilities.NDDivision.FloorDivInt64), typeof(long));
            public static readonly MethodInfo FloorDivUInt64 = NDDiv(nameof(Utilities.NDDivision.FloorDivUInt64), typeof(ulong));
            public static readonly MethodInfo FloorDivSingle = NDDiv(nameof(Utilities.NDDivision.FloorDivSingle), typeof(float));
            public static readonly MethodInfo FloorDivDouble = NDDiv(nameof(Utilities.NDDivision.FloorDivDouble), typeof(double));

            public static readonly MethodInfo RemSByte = NDDiv(nameof(Utilities.NDDivision.RemSByte), typeof(sbyte));
            public static readonly MethodInfo RemByte = NDDiv(nameof(Utilities.NDDivision.RemByte), typeof(byte));
            public static readonly MethodInfo RemInt16 = NDDiv(nameof(Utilities.NDDivision.RemInt16), typeof(short));
            public static readonly MethodInfo RemUInt16 = NDDiv(nameof(Utilities.NDDivision.RemUInt16), typeof(ushort));
            public static readonly MethodInfo RemChar = NDDiv(nameof(Utilities.NDDivision.RemChar), typeof(char));
            public static readonly MethodInfo RemInt32 = NDDiv(nameof(Utilities.NDDivision.RemInt32), typeof(int));
            public static readonly MethodInfo RemUInt32 = NDDiv(nameof(Utilities.NDDivision.RemUInt32), typeof(uint));
            public static readonly MethodInfo RemInt64 = NDDiv(nameof(Utilities.NDDivision.RemInt64), typeof(long));
            public static readonly MethodInfo RemUInt64 = NDDiv(nameof(Utilities.NDDivision.RemUInt64), typeof(ulong));
            public static readonly MethodInfo RemSingle = NDDiv(nameof(Utilities.NDDivision.RemSingle), typeof(float));
            public static readonly MethodInfo RemDouble = NDDiv(nameof(Utilities.NDDivision.RemDouble), typeof(double));

            // Decimal conversion methods (to decimal)
            public static readonly MethodInfo DecimalImplicitFromInt = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(int) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(int)");
            public static readonly MethodInfo DecimalImplicitFromByte = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(byte) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(byte)");
            public static readonly MethodInfo DecimalImplicitFromSByte = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(sbyte) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(sbyte)");
            public static readonly MethodInfo DecimalImplicitFromShort = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(short) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(short)");
            public static readonly MethodInfo DecimalImplicitFromUShort = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(ushort) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(ushort)");
            public static readonly MethodInfo DecimalImplicitFromUInt = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(uint) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(uint)");
            public static readonly MethodInfo DecimalImplicitFromLong = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(long) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(long)");
            public static readonly MethodInfo DecimalImplicitFromULong = typeof(decimal).GetMethod("op_Implicit", new[] { typeof(ulong) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Implicit(ulong)");
            public static readonly MethodInfo DecimalExplicitFromFloat = typeof(decimal).GetMethod("op_Explicit", new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Explicit(float)");
            public static readonly MethodInfo DecimalExplicitFromDouble = typeof(decimal).GetMethod("op_Explicit", new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Explicit(double)");

            // Decimal conversion methods (from decimal)
            public static readonly MethodInfo DecimalToByte = typeof(decimal).GetMethod("ToByte", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToByte");
            public static readonly MethodInfo DecimalToSByte = typeof(decimal).GetMethod("ToSByte", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToSByte");
            public static readonly MethodInfo DecimalToInt16 = typeof(decimal).GetMethod("ToInt16", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToInt16");
            public static readonly MethodInfo DecimalToUInt16 = typeof(decimal).GetMethod("ToUInt16", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToUInt16");
            public static readonly MethodInfo DecimalToInt32 = typeof(decimal).GetMethod("ToInt32", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToInt32");
            public static readonly MethodInfo DecimalToUInt32 = typeof(decimal).GetMethod("ToUInt32", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToUInt32");
            public static readonly MethodInfo DecimalToInt64 = typeof(decimal).GetMethod("ToInt64", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToInt64");
            public static readonly MethodInfo DecimalToUInt64 = typeof(decimal).GetMethod("ToUInt64", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToUInt64");
            public static readonly MethodInfo DecimalToSingle = typeof(decimal).GetMethod("ToSingle", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToSingle");
            public static readonly MethodInfo DecimalToDouble = typeof(decimal).GetMethod("ToDouble", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "ToDouble");

            // Decimal operator methods
            public static readonly MethodInfo DecimalOpAddition = typeof(decimal).GetMethod("op_Addition",
                BindingFlags.Public | BindingFlags.Static, null, new[] { typeof(decimal), typeof(decimal) }, null)
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Addition");
            public static readonly MethodInfo DecimalOpSubtraction = typeof(decimal).GetMethod("op_Subtraction",
                BindingFlags.Public | BindingFlags.Static, null, new[] { typeof(decimal), typeof(decimal) }, null)
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Subtraction");
            public static readonly MethodInfo DecimalOpMultiply = typeof(decimal).GetMethod("op_Multiply",
                BindingFlags.Public | BindingFlags.Static, null, new[] { typeof(decimal), typeof(decimal) }, null)
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Multiply");
            public static readonly MethodInfo DecimalOpDivision = typeof(decimal).GetMethod("op_Division",
                BindingFlags.Public | BindingFlags.Static, null, new[] { typeof(decimal), typeof(decimal) }, null)
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Division");
            public static readonly MethodInfo DecimalFloor = typeof(decimal).GetMethod(nameof(decimal.Floor),
                BindingFlags.Public | BindingFlags.Static, null, new[] { typeof(decimal) }, null)
                ?? throw new MissingMethodException(typeof(decimal).FullName, nameof(decimal.Floor));

            // NumSharp.Utilities.DecimalMath methods
            public static readonly MethodInfo DecimalMathPow = typeof(Utilities.DecimalMath).GetMethod(
                nameof(Utilities.DecimalMath.Pow), BindingFlags.Public | BindingFlags.Static, null,
                new[] { typeof(decimal), typeof(decimal) }, null)
                ?? throw new MissingMethodException(typeof(Utilities.DecimalMath).FullName, nameof(Utilities.DecimalMath.Pow));
            public static readonly MethodInfo DecimalMathATan2 = typeof(Utilities.DecimalMath).GetMethod(
                nameof(Utilities.DecimalMath.ATan2), BindingFlags.Public | BindingFlags.Static, null,
                new[] { typeof(decimal), typeof(decimal) }, null)
                ?? throw new MissingMethodException(typeof(Utilities.DecimalMath).FullName, nameof(Utilities.DecimalMath.ATan2));

            // Decimal fields
            public static readonly FieldInfo DecimalZero = typeof(decimal).GetField(nameof(decimal.Zero))
                ?? throw new MissingFieldException(typeof(decimal).FullName, nameof(decimal.Zero));
            public static readonly FieldInfo DecimalOne = typeof(decimal).GetField(nameof(decimal.One))
                ?? throw new MissingFieldException(typeof(decimal).FullName, nameof(decimal.One));
            public static readonly FieldInfo DecimalMinValue = typeof(decimal).GetField(nameof(decimal.MinValue))
                ?? throw new MissingFieldException(typeof(decimal).FullName, nameof(decimal.MinValue));
            public static readonly FieldInfo DecimalMaxValue = typeof(decimal).GetField(nameof(decimal.MaxValue))
                ?? throw new MissingFieldException(typeof(decimal).FullName, nameof(decimal.MaxValue));

            // Additional decimal operator methods
            public static readonly MethodInfo DecimalOpUnaryNegation = typeof(decimal).GetMethod("op_UnaryNegation", new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_UnaryNegation");
            public static readonly MethodInfo DecimalOpEquality = typeof(decimal).GetMethod("op_Equality", new[] { typeof(decimal), typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, "op_Equality");
            public static readonly MethodInfo DecimalTruncate = typeof(decimal).GetMethod(nameof(decimal.Truncate), new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(decimal).FullName, nameof(decimal.Truncate));

            // Math methods for decimal
            public static readonly MethodInfo MathAbsDecimal = typeof(Math).GetMethod(nameof(Math.Abs), new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(Math).FullName, "Abs(decimal)");
            public static readonly MethodInfo MathSignDecimal = typeof(Math).GetMethod(nameof(Math.Sign), new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(Math).FullName, "Sign(decimal)");
            public static readonly MethodInfo MathCeilingDecimal = typeof(Math).GetMethod(nameof(Math.Ceiling), new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(Math).FullName, "Ceiling(decimal)");
            public static readonly MethodInfo MathFloorDecimal = typeof(Math).GetMethod(nameof(Math.Floor), new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(Math).FullName, "Floor(decimal)");
            public static readonly MethodInfo MathRoundDecimal = typeof(Math).GetMethod(nameof(Math.Round), new[] { typeof(decimal) })
                ?? throw new MissingMethodException(typeof(Math).FullName, "Round(decimal)");

            // Math methods for double
            public static readonly MethodInfo MathAbsDouble = typeof(Math).GetMethod(nameof(Math.Abs), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, "Abs(double)");
            public static readonly MethodInfo MathExp = typeof(Math).GetMethod(nameof(Math.Exp), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, nameof(Math.Exp));
            public static readonly MethodInfo MathLog = typeof(Math).GetMethod(nameof(Math.Log), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, nameof(Math.Log));
            public static readonly MethodInfo MathCbrt = typeof(Math).GetMethod(nameof(Math.Cbrt), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, nameof(Math.Cbrt));

            // MathF methods for float
            public static readonly MethodInfo MathFAbsFloat = typeof(MathF).GetMethod(nameof(MathF.Abs), new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(MathF).FullName, nameof(MathF.Abs));
            public static readonly MethodInfo MathFSign = typeof(MathF).GetMethod(nameof(MathF.Sign), new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(MathF).FullName, nameof(MathF.Sign));

            // Math.Sign methods
            public static readonly MethodInfo MathSignDouble = typeof(Math).GetMethod(nameof(Math.Sign), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, "Sign(double)");

            // IsNaN / IsInfinity / IsFinite methods
            public static readonly MethodInfo FloatIsNaN = typeof(float).GetMethod(nameof(float.IsNaN), new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(float).FullName, nameof(float.IsNaN));
            public static readonly MethodInfo DoubleIsNaN = typeof(double).GetMethod(nameof(double.IsNaN), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(double).FullName, nameof(double.IsNaN));
            public static readonly MethodInfo DoubleIsInfinity = typeof(double).GetMethod(nameof(double.IsInfinity), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(double).FullName, nameof(double.IsInfinity));
            public static readonly MethodInfo DoubleIsFinite = typeof(double).GetMethod(nameof(double.IsFinite), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(double).FullName, nameof(double.IsFinite));
            public static readonly MethodInfo MathCopySign = typeof(Math).GetMethod(nameof(Math.CopySign), new[] { typeof(double), typeof(double) })
                ?? throw new MissingMethodException(typeof(Math).FullName, nameof(Math.CopySign));

            // Unsafe methods
            public static readonly MethodInfo UnsafeInitBlockUnaligned = typeof(Unsafe).GetMethod(nameof(Unsafe.InitBlockUnaligned),
                new[] { typeof(void*), typeof(byte), typeof(uint) })
                ?? throw new MissingMethodException(typeof(Unsafe).FullName, nameof(Unsafe.InitBlockUnaligned));

            // (Vector256 operator methods used by MatMul moved to VectorMethodCache.Operator.)

            // Half conversion methods (Half is a struct with operator methods, not IConvertible)
            public static readonly MethodInfo HalfToDouble = typeof(Half).GetMethods(BindingFlags.Public | BindingFlags.Static)
                .First(m => m.Name == "op_Explicit" && m.ReturnType == typeof(double) && m.GetParameters().Length == 1 && m.GetParameters()[0].ParameterType == typeof(Half));
            public static readonly MethodInfo DoubleToHalf = typeof(Half).GetMethods(BindingFlags.Public | BindingFlags.Static)
                .First(m => m.Name == "op_Explicit" && m.ReturnType == typeof(Half) && m.GetParameters().Length == 1 && m.GetParameters()[0].ParameterType == typeof(double));
            // Half <-> float (F16C-lowered by the JIT: [Intrinsic] op_Explicit → vcvtph2ps / vcvtps2ph).
            // These are the fast, NumPy-HALF-loop-matching conversions: NumPy's npy_half loops widen
            // to float32 (npy_half_to_float), compute the CRT float function, and round back
            // (npy_float_to_half) — see EmitUnaryHalfViaFloat. The double bridge above is 2 conversions
            // deep (Half→float→double) and drags in the slower double math functions.
            public static readonly MethodInfo HalfToFloat = typeof(Half).GetMethods(BindingFlags.Public | BindingFlags.Static)
                .First(m => m.Name == "op_Explicit" && m.ReturnType == typeof(float) && m.GetParameters().Length == 1 && m.GetParameters()[0].ParameterType == typeof(Half));
            public static readonly MethodInfo FloatToHalf = typeof(Half).GetMethods(BindingFlags.Public | BindingFlags.Static)
                .First(m => m.Name == "op_Explicit" && m.ReturnType == typeof(Half) && m.GetParameters().Length == 1 && m.GetParameters()[0].ParameterType == typeof(float));
            // float.Exp2 == the CRT exp2f NumPy's HALF_exp2 loop calls (npy_exp2f). Used ONLY on the
            // Half path: it is the correctly-rounded software 2^x, byte-identical to NumPy's half loop
            // (finite AND NaN payload). The float32-array path uses NDFloatMath.Exp2 (the fast SIMD
            // ≤1-ULP kernel) instead — but that kernel's scalar entry is far slower here than float.Exp2
            // and blanks the NaN payload, so the half loop deliberately keeps the CRT.
            public static readonly MethodInfo SingleExp2Crt = typeof(float).GetMethod("Exp2", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(float).FullName, "Exp2(float)");
            public static readonly MethodInfo HalfIsNaN = typeof(Half).GetMethod("IsNaN", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "IsNaN");

            // Half static properties (NaN, Zero, PositiveInfinity, NegativeInfinity are properties, not fields)
            public static readonly MethodInfo HalfNaN = typeof(Half).GetProperty("NaN", BindingFlags.Public | BindingFlags.Static)!.GetGetMethod()
                ?? throw new MissingMethodException(typeof(Half).FullName, "NaN");
            public static readonly MethodInfo HalfZero = typeof(Half).GetProperty("Zero", BindingFlags.Public | BindingFlags.Static)!.GetGetMethod()
                ?? throw new MissingMethodException(typeof(Half).FullName, "Zero");
            public static readonly MethodInfo HalfPositiveInfinity = typeof(Half).GetProperty("PositiveInfinity", BindingFlags.Public | BindingFlags.Static)!.GetGetMethod()
                ?? throw new MissingMethodException(typeof(Half).FullName, "PositiveInfinity");
            public static readonly MethodInfo HalfNegativeInfinity = typeof(Half).GetProperty("NegativeInfinity", BindingFlags.Public | BindingFlags.Static)!.GetGetMethod()
                ?? throw new MissingMethodException(typeof(Half).FullName, "NegativeInfinity");

            // Complex methods and fields (Complex uses static fields, not properties).
            // ComplexAbs routes through NDComplexMath.Abs (npy_cabs / C99 hypot semantics) rather
            // than Complex.Abs directly: the BCL's private Hypot returns NaN for abs(NaN+inf*i) on
            // net8.0, where NumPy returns +inf. The helper defers to Complex.Abs for every
            // finite/NaN-only input, so magnitudes that already match NumPy stay bit-identical.
            public static readonly MethodInfo ComplexAbs = typeof(Utilities.NDComplexMath).GetMethod("Abs", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(Utilities.NDComplexMath).FullName, "Abs");
            public static readonly MethodInfo ComplexDivisionByDouble = typeof(System.Numerics.Complex).GetMethod("op_Division", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex), typeof(double) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "op_Division(Complex, double)");
            public static readonly FieldInfo ComplexZero = typeof(System.Numerics.Complex).GetField("Zero", BindingFlags.Public | BindingFlags.Static)
                ?? throw new MissingFieldException(typeof(System.Numerics.Complex).FullName, "Zero");
            public static readonly FieldInfo ComplexOne = typeof(System.Numerics.Complex).GetField("One", BindingFlags.Public | BindingFlags.Static)
                ?? throw new MissingFieldException(typeof(System.Numerics.Complex).FullName, "One");
            public static readonly ConstructorInfo ComplexCtor = typeof(System.Numerics.Complex).GetConstructor(new[] { typeof(double), typeof(double) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, ".ctor(double, double)");

            // Complex binary operator methods
            public static readonly MethodInfo ComplexOpAddition = typeof(System.Numerics.Complex).GetMethod("op_Addition", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex), typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "op_Addition");
            public static readonly MethodInfo ComplexOpMultiply = typeof(System.Numerics.Complex).GetMethod("op_Multiply", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex), typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "op_Multiply");

            // Complex unary operator methods
            public static readonly MethodInfo ComplexNegate = typeof(System.Numerics.Complex).GetMethod("op_UnaryNegation", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "op_UnaryNegation");
            // np.conjugate: Complex.Conjugate(z) flips the sign of the imaginary part (re - i*im).
            public static readonly MethodInfo ComplexConjugate = typeof(System.Numerics.Complex).GetMethod("Conjugate", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "Conjugate");
            // Sqrt/Exp/Sin/Cos/Tan route through NDComplexMath (not Complex.* directly): the BCL
            // matches NumPy on finite interiors but diverges at the C99 edges (non-finite, branch-cut
            // signs, signed zeros). NDComplexMath delegates the interior to the BCL and adds the fixups.
            public static readonly MethodInfo ComplexSqrt = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Sqrt", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Sqrt");
            public static readonly MethodInfo ComplexExp = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Exp", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Exp");
            // float32 exp routes through NDFloatMath.Exp (port of NumPy's simd_exp_FLOAT), NOT
            // MathF.Exp: the platform libm is ~correctly rounded while NumPy's Remez/Cody-Waite
            // kernel carries up to 2.52 ULP of its own error, so "more accurate" reads as a
            // divergence. The port is bit-exact against NumPy 2.4.2 over all 2^32 inputs.
            public static readonly MethodInfo SingleExp = typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Exp", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDFloatMath).FullName, "Exp");
            // The same kernel a whole register at a time, resolved for the width detected at
            // startup. NULLABLE by design (unlike the entries above): a host with no SIMD has no
            // width to bind, and the eligibility gate (ExpVectorSimdAvailable) reads this field to
            // decide, so an absent overload must be a false rather than a type-load failure.
            public static readonly MethodInfo SingleExpVector =
                VectorBits == 0 ? null
                : typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Exp", BindingFlags.Public | BindingFlags.Static,
                    new[] { VectorMethodCache.V(VectorBits, typeof(float)) });
            // float32 exp2 (2^x): NDFloatMath.Exp2, a fast replacement for (float)Math.Pow(2, x). NOT
            // a NumPy-kernel port — the win-amd64 wheel runs a scalar exp2f (its SVML vector kernel is
            // AVX-512/Linux-gated) — but ≤1 ULP and ~2.4x faster than that scalar loop. Same scalar/
            // vector split as SingleExp; the vector overload is nullable for the no-SIMD host.
            public static readonly MethodInfo SingleExp2 = typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Exp2", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDFloatMath).FullName, "Exp2");
            public static readonly MethodInfo SingleExp2Vector =
                VectorBits == 0 ? null
                : typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Exp2", BindingFlags.Public | BindingFlags.Static,
                    new[] { VectorMethodCache.V(VectorBits, typeof(float)) });
            // float32 log: the sibling port (NumPy's simd_log_FLOAT), same reasoning as SingleExp.
            public static readonly MethodInfo SingleLog = typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Log", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDFloatMath).FullName, "Log");
            public static readonly MethodInfo SingleLogVector =
                VectorBits == 0 ? null
                : typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Log", BindingFlags.Public | BindingFlags.Static,
                    new[] { VectorMethodCache.V(VectorBits, typeof(float)) });
            // float32 sin: port of NumPy's simd_sincos_f32 (loops_trigonometric).
            public static readonly MethodInfo SingleSin = typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Sin", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDFloatMath).FullName, "Sin");
            public static readonly MethodInfo SingleSinVector =
                VectorBits == 0 ? null
                : typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Sin", BindingFlags.Public | BindingFlags.Static,
                    new[] { VectorMethodCache.V(VectorBits, typeof(float)) });
            // float32 cos: port of NumPy's simd_sincos_f32 (loops_trigonometric).
            public static readonly MethodInfo SingleCos = typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Cos", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDFloatMath).FullName, "Cos");
            public static readonly MethodInfo SingleCosVector =
                VectorBits == 0 ? null
                : typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Cos", BindingFlags.Public | BindingFlags.Static,
                    new[] { VectorMethodCache.V(VectorBits, typeof(float)) });
            // tanh is the FIRST of these ports to cover float64 as well: NumPy ships its own
            // kernel (loops_hyperbolic) at both widths, so unlike exp/log/sin/cos — where the
            // platform libm already agrees with NumPy at f8 — BOTH the float and double BCL calls
            // diverged and both are replaced. Non-nullable: the scalar entry points always exist.
            public static readonly MethodInfo SingleTanh = typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Tanh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDFloatMath).FullName, "Tanh");
            public static readonly MethodInfo DoubleTanh = typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Tanh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDFloatMath).FullName, "Tanh");
            // Vector form is float32-only (the f64 lookup would cost 18 gathers per 4 lanes, and the
            // scalar f64 kernel already outruns NumPy). Nullable for the usual reason — a host with
            // no SIMD has no width to bind.
            public static readonly MethodInfo SingleTanhVector =
                VectorBits == 0 ? null
                : typeof(NumSharp.Utilities.NDFloatMath).GetMethod("Tanh", BindingFlags.Public | BindingFlags.Static,
                    new[] { VectorMethodCache.V(VectorBits, typeof(float)) });
            // ComplexLog routes through NDComplexMath.Log (full npy_clog port): Complex.Log drops the
            // real part to 0 near |z|=1 (it lacks clog's log1p path). Reused by the Log2 composition
            // and by NDComplexMath.Log10/Log1p.
            public static readonly MethodInfo ComplexLog = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Log", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Log");
            public static readonly MethodInfo ComplexSin = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Sin", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Sin");
            public static readonly MethodInfo ComplexCos = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Cos", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Cos");
            public static readonly MethodInfo ComplexTan = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Tan", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Tan");
            // Hyperbolic and inverse-trig route through NDComplexMath (not Complex.* directly): the BCL
            // matches NumPy only on finite interiors; NDComplexMath adds the C99 Annex G non-finite
            // tables and signed-zero/branch-cut fixups so every input matches NumPy.
            public static readonly MethodInfo ComplexSinh = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Sinh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Sinh");
            public static readonly MethodInfo ComplexCosh = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Cosh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Cosh");
            public static readonly MethodInfo ComplexTanh = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Tanh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Tanh");
            public static readonly MethodInfo ComplexAsin = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Asin", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Asin");
            public static readonly MethodInfo ComplexAcos = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Acos", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Acos");
            public static readonly MethodInfo ComplexAtan = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Atan", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Atan");
            // Inverse hyperbolic: NumSharp derives casinh/cacosh from the byte-exact Asin/Acos via the
            // exact msun involution, and exposes the already-ported catanh (which drives Atan).
            public static readonly MethodInfo ComplexAsinh = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Asinh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Asinh");
            public static readonly MethodInfo ComplexAcosh = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Acosh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Acosh");
            public static readonly MethodInfo ComplexAtanh = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Atanh", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Atanh");
            public static readonly MethodInfo ComplexPow = typeof(System.Numerics.Complex).GetMethod("Pow", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex), typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "Pow");
            // Log10/Reciprocal/Log1p route through NDComplexMath (Complex.Log10 drifts past 1 ULP from
            // NumPy; Complex.op_Division / Complex.One+z drop NumPy's signed zeros).
            public static readonly MethodInfo ComplexLog10 = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Log10", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Log10");
            public static readonly MethodInfo ComplexReciprocal = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Reciprocal", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Reciprocal");
            // Square routes through NDComplexMath (FMA-contracted z*z): Complex.op_Multiply lacks FMA,
            // so it loses NumPy's square(1e-10+1e-10i).real = -2.275e-37 and turns the 1e300 overflow
            // into NaN instead of NumPy's -inf.
            public static readonly MethodInfo ComplexSquare = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Square", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Square");
            public static readonly MethodInfo ComplexLog1p = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Log1p", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Log1p");
            public static readonly MethodInfo ComplexExp2 = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Exp2", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Exp2");
            // Expm1 routes through NDComplexMath (nc_expm1 formula); Complex.Exp(z)-1 mis-handles the
            // non-finite imaginary parts (e.g. expm1(+Inf+0i).imag must be NaN = exp(+Inf)*sin(0)).
            public static readonly MethodInfo ComplexExpm1 = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Expm1", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Expm1");
            // Complex doesn't have Log2/Exp2/Log1p directly — composed via Log(z, 2), Pow(2, z),
            // NDComplexMath.Log1p in EmitUnaryComplexOperation.
            public static readonly MethodInfo ComplexLogBase = typeof(System.Numerics.Complex).GetMethod("Log", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex), typeof(double) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "Log(Complex, double)");
            public static readonly MethodInfo ComplexOpSubtraction = typeof(System.Numerics.Complex).GetMethod("op_Subtraction", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex), typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "op_Subtraction");

            // Complex instance property getters — called via Ldloca + Call (struct instance method
            // requires a managed reference for 'this').
            public static readonly MethodInfo ComplexGetReal = typeof(System.Numerics.Complex)
                .GetProperty("Real", BindingFlags.Public | BindingFlags.Instance)!.GetGetMethod()
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "get_Real");
            public static readonly MethodInfo ComplexGetImaginary = typeof(System.Numerics.Complex)
                .GetProperty("Imaginary", BindingFlags.Public | BindingFlags.Instance)!.GetGetMethod()
                ?? throw new MissingMethodException(typeof(System.Numerics.Complex).FullName, "get_Imaginary");

            // Field handle for the runtime-computed 1/ln(2) constant used by Complex log2 inline IL.
            public static readonly FieldInfo LogE_Inv_Ln2Field = typeof(DirectILKernelGenerator)
                .GetField(nameof(DirectILKernelGenerator.LogE_Inv_Ln2), BindingFlags.NonPublic | BindingFlags.Static)
                ?? throw new MissingFieldException(typeof(DirectILKernelGenerator).FullName, nameof(DirectILKernelGenerator.LogE_Inv_Ln2));

            // Half unary operator methods. Negate is OUR sign-bit-flip helper, not
            // Half.op_UnaryNegation (that operator does a (Half)(-(float)h) roundtrip measured
            // 7.3× slower — it made f16 negate the worst cell in the elementwise matrix, ~0.14×).
            public static readonly MethodInfo HalfNegate = typeof(DirectILKernelGenerator).GetMethod(nameof(DirectILKernelGenerator.NegateHalf), BindingFlags.NonPublic | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(DirectILKernelGenerator).FullName, "NegateHalf");
            public static readonly MethodInfo HalfSqrt = typeof(Half).GetMethod("Sqrt", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Sqrt");
            public static readonly MethodInfo HalfSin = typeof(Half).GetMethod("Sin", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Sin");
            public static readonly MethodInfo HalfCos = typeof(Half).GetMethod("Cos", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Cos");
            public static readonly MethodInfo HalfExp = typeof(Half).GetMethod("Exp", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Exp");
            public static readonly MethodInfo HalfLog = typeof(Half).GetMethod("Log", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Log");
            public static readonly MethodInfo HalfFloor = typeof(Half).GetMethod("Floor", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Floor");
            public static readonly MethodInfo HalfCeiling = typeof(Half).GetMethod("Ceiling", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Ceiling");
            public static readonly MethodInfo HalfTruncate = typeof(Half).GetMethod("Truncate", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Truncate");
            public static readonly MethodInfo HalfAbs = typeof(Half).GetMethod("Abs", BindingFlags.Public | BindingFlags.Static, new[] { typeof(Half) })
                ?? throw new MissingMethodException(typeof(Half).FullName, "Abs");
            // Sinh/Cosh/Tanh/ASin/ACos/ATan/Asinh/Acosh/Atanh/Tan/Cbrt/Log2/Log10/Exp2/Reciprocal/
            // Square on Half now take the float32 fast path (EmitUnaryHalfViaFloat: `(Half)Xf((float)h)`,
            // the F16C round-trip that mirrors NumPy's npy_half loop) rather than a BCL Half.* call or
            // the double bridge — so the corresponding Half.* MethodInfos (HalfTan/Log10/Log2/Cbrt/Exp2/
            // Asinh/Acosh/Atanh) were removed. Only the ops still on a dedicated Half path keep a cached
            // method (Sqrt/Sin/Cos/Exp/Log/Floor/Ceiling/Truncate/Abs/Negate above).
            // Note: .NET's Half exposes log1p as LogP1 and expm1 as ExpM1 (IFloatingPointIeee754<Half>).
            // Half.LogP1/ExpM1 lose subnormal precision because internally they compute (1 + x) in
            // Half, which rounds x < Half.Epsilon (≈ 2^-11) to 0. NumPy promotes to a higher-precision
            // intermediate before log1p/expm1, then casts back — we replicate that with double
            // (via existing HalfToDouble / DoubleToHalf op_Explicit helpers).
            public static readonly MethodInfo DoubleLogP1 = typeof(double)
                .GetMethod("LogP1", BindingFlags.Public | BindingFlags.Static, new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(double).FullName, "LogP1");
            public static readonly MethodInfo DoubleExpM1 = typeof(double)
                .GetMethod("ExpM1", BindingFlags.Public | BindingFlags.Static, new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(double).FullName, "ExpM1");
        }

        #endregion

        /// <summary>
        /// Get the Vector container type (Vector128, Vector256, or Vector512).
        /// Thin wrapper around <see cref="VectorMethodCache.Container"/> preserved
        /// so the existing call sites in the partial-class files keep their shape.
        /// </summary>
        internal static Type GetVectorContainerType() => VectorMethodCache.Container(VectorBits);

        /// <summary>
        /// Get the Vector{Width}&lt;T&gt; generic type.
        /// Thin wrapper around <see cref="VectorMethodCache.V"/>.
        /// </summary>
        internal static Type GetVectorType(Type elementType) => VectorMethodCache.V(VectorBits, elementType);

        /// <summary>
        /// Resolve a NonPublic|Static helper method on <see cref="DirectILKernelGenerator"/> by name.
        /// Throws <see cref="MissingMethodException"/> on miss rather than returning <c>null</c>,
        /// so renamed or removed helpers fail loudly at IL emit time rather than as NREs later.
        /// </summary>
        internal static MethodInfo GetHelper(string name)
            => typeof(DirectILKernelGenerator).GetMethod(name, BindingFlags.NonPublic | BindingFlags.Static)
               ?? throw new MissingMethodException(typeof(DirectILKernelGenerator).FullName, name);

        /// <summary>
        /// Same as <see cref="GetHelper"/> but immediately closes a single generic argument.
        /// </summary>
        internal static MethodInfo GetGenericHelper(string name, Type genericArg)
            => GetHelper(name).MakeGenericMethod(genericArg);

        #region NPTypeCode-Based IL Helpers

        /// <summary>
        /// Get size in bytes for NPTypeCode.
        /// </summary>
        internal static int GetTypeSize(NPTypeCode type)
        {
            return type switch
            {
                NPTypeCode.Boolean => 1,
                NPTypeCode.Byte => 1,
                NPTypeCode.SByte => 1,
                NPTypeCode.Int16 => 2,
                NPTypeCode.UInt16 => 2,
                NPTypeCode.Half => 2,
                NPTypeCode.Int32 => 4,
                NPTypeCode.UInt32 => 4,
                NPTypeCode.Int64 => 8,
                NPTypeCode.UInt64 => 8,
                NPTypeCode.Char => 2,
                NPTypeCode.Single => 4,
                NPTypeCode.Double => 8,
                NPTypeCode.Decimal => 16,
                NPTypeCode.Complex => 16,
                _ => throw new NotSupportedException($"Type {type} not supported")
            };
        }

        /// <summary>
        /// Get CLR Type for NPTypeCode.
        /// </summary>
        internal static Type GetClrType(NPTypeCode type)
        {
            return type switch
            {
                NPTypeCode.Boolean => typeof(bool),
                NPTypeCode.Byte => typeof(byte),
                NPTypeCode.SByte => typeof(sbyte),
                NPTypeCode.Int16 => typeof(short),
                NPTypeCode.UInt16 => typeof(ushort),
                NPTypeCode.Half => typeof(Half),
                NPTypeCode.Int32 => typeof(int),
                NPTypeCode.UInt32 => typeof(uint),
                NPTypeCode.Int64 => typeof(long),
                NPTypeCode.UInt64 => typeof(ulong),
                NPTypeCode.Char => typeof(char),
                NPTypeCode.Single => typeof(float),
                NPTypeCode.Double => typeof(double),
                NPTypeCode.Decimal => typeof(decimal),
                NPTypeCode.Complex => typeof(System.Numerics.Complex),
                _ => throw new NotSupportedException($"Type {type} not supported")
            };
        }

        /// <summary>
        /// Check if type supports SIMD operations (V128/V256/V512).
        /// </summary>
        internal static bool CanUseSimd(NPTypeCode type)
        {
            if (VectorBits == 0) return false;  // No SIMD hardware

            return type switch
            {
                NPTypeCode.Byte or NPTypeCode.SByte => true,
                NPTypeCode.Int16 or NPTypeCode.UInt16 => true,
                NPTypeCode.Int32 or NPTypeCode.UInt32 => true,
                NPTypeCode.Int64 or NPTypeCode.UInt64 => true,
                NPTypeCode.Single or NPTypeCode.Double => true,
                _ => false  // Boolean, Char, Decimal, Half, Complex
            };
        }

        /// <summary>
        /// Get vector element count for type (adapts to V128/V256/V512).
        /// </summary>
        internal static int GetVectorCount(NPTypeCode type)
        {
            if (VectorBits == 0) return 1;  // Scalar fallback
            return VectorBytes / GetTypeSize(type);
        }

        /// <summary>
        /// Emit load indirect for NPTypeCode.
        /// </summary>
        internal static void EmitLoadIndirect(ILGenerator il, NPTypeCode type)
        {
            switch (type)
            {
                case NPTypeCode.Boolean:
                case NPTypeCode.Byte:
                    il.Emit(OpCodes.Ldind_U1);
                    break;
                case NPTypeCode.SByte:
                    il.Emit(OpCodes.Ldind_I1);
                    break;
                case NPTypeCode.Int16:
                    il.Emit(OpCodes.Ldind_I2);
                    break;
                case NPTypeCode.UInt16:
                case NPTypeCode.Char:
                    il.Emit(OpCodes.Ldind_U2);
                    break;
                case NPTypeCode.Half:
                    il.Emit(OpCodes.Ldobj, typeof(Half));
                    break;
                case NPTypeCode.Int32:
                    il.Emit(OpCodes.Ldind_I4);
                    break;
                case NPTypeCode.UInt32:
                    il.Emit(OpCodes.Ldind_U4);
                    break;
                case NPTypeCode.Int64:
                case NPTypeCode.UInt64:
                    il.Emit(OpCodes.Ldind_I8);
                    break;
                case NPTypeCode.Single:
                    il.Emit(OpCodes.Ldind_R4);
                    break;
                case NPTypeCode.Double:
                    il.Emit(OpCodes.Ldind_R8);
                    break;
                case NPTypeCode.Decimal:
                    il.Emit(OpCodes.Ldobj, typeof(decimal));
                    break;
                case NPTypeCode.Complex:
                    il.Emit(OpCodes.Ldobj, typeof(System.Numerics.Complex));
                    break;
                default:
                    throw new NotSupportedException($"Type {type} not supported for ldind");
            }
        }

        /// <summary>
        /// Emit store indirect for NPTypeCode.
        /// </summary>
        internal static void EmitStoreIndirect(ILGenerator il, NPTypeCode type)
        {
            switch (type)
            {
                case NPTypeCode.Boolean:
                case NPTypeCode.Byte:
                case NPTypeCode.SByte:
                    il.Emit(OpCodes.Stind_I1);
                    break;
                case NPTypeCode.Int16:
                case NPTypeCode.UInt16:
                case NPTypeCode.Char:
                    il.Emit(OpCodes.Stind_I2);
                    break;
                case NPTypeCode.Half:
                    il.Emit(OpCodes.Stobj, typeof(Half));
                    break;
                case NPTypeCode.Int32:
                case NPTypeCode.UInt32:
                    il.Emit(OpCodes.Stind_I4);
                    break;
                case NPTypeCode.Int64:
                case NPTypeCode.UInt64:
                    il.Emit(OpCodes.Stind_I8);
                    break;
                case NPTypeCode.Single:
                    il.Emit(OpCodes.Stind_R4);
                    break;
                case NPTypeCode.Double:
                    il.Emit(OpCodes.Stind_R8);
                    break;
                case NPTypeCode.Decimal:
                    il.Emit(OpCodes.Stobj, typeof(decimal));
                    break;
                case NPTypeCode.Complex:
                    il.Emit(OpCodes.Stobj, typeof(System.Numerics.Complex));
                    break;
                default:
                    throw new NotSupportedException($"Type {type} not supported for stind");
            }
        }

        /// <summary>
        /// Emit type conversion from source to target type.
        /// </summary>
        internal static void EmitConvertTo(ILGenerator il, NPTypeCode from, NPTypeCode to)
        {
            if (from == to)
                return; // No conversion needed

            // A Boolean's numeric value is exactly 0 or 1 — never its raw storage byte. A bool buffer
            // can legally hold non-0/1 bytes (np.frombuffer is a zero-copy VIEW, like NumPy; interop
            // wraps a foreign buffer), so normalize nonzero->1 BEFORE widening to any numeric type.
            // Without this, sum/mean/var/std over such a buffer accumulate the raw bytes (e.g. byte 255
            // contributes 255) instead of counting True. Mirrors the to==Boolean '!= 0' below, applied
            // to the source side. Idempotent for proper 0/1 bools (0->0, 1->1).
            if (from == NPTypeCode.Boolean) // to != Boolean here (from == to already returned)
            {
                il.Emit(OpCodes.Ldc_I4_0);
                il.Emit(OpCodes.Cgt_Un);    // value = (value != 0) ? 1 : 0
                from = NPTypeCode.Byte;      // continue the conversion from a clean 0/1 byte
            }

            // Special case: decimal conversions require method calls
            if (from == NPTypeCode.Decimal || to == NPTypeCode.Decimal)
            {
                EmitDecimalConversion(il, from, to);
                return;
            }

            // Special case: Half and Complex require method calls
            if (from == NPTypeCode.Half || from == NPTypeCode.Complex || to == NPTypeCode.Half || to == NPTypeCode.Complex)
            {
                EmitHalfOrComplexConversion(il, from, to);
                return;
            }

            // float -> integer / char: the IL conv.* opcodes SATURATE on overflow and yield 0
            // for NaN, but NumPy truncates toward zero and yields the type-min sentinel
            // (int.MinValue & co.) on NaN/overflow, then wraps to narrower widths. Route through
            // the Converts.* table, which is bit-exact with NumPy across the full cast matrix.
            if ((from == NPTypeCode.Single || from == NPTypeCode.Double)
                && to != NPTypeCode.Boolean && to != NPTypeCode.Single && to != NPTypeCode.Double)
            {
                il.EmitCall(OpCodes.Call, ConvertsFloatToIntMethod(from, to), null);
                return;
            }

            // For numeric types, use conv.* opcodes
            switch (to)
            {
                case NPTypeCode.Boolean:
                    // Convert to bool: != 0
                    il.Emit(OpCodes.Ldc_I4_0);
                    il.Emit(OpCodes.Cgt_Un);
                    break;
                case NPTypeCode.Byte:
                    il.Emit(OpCodes.Conv_U1);
                    break;
                case NPTypeCode.SByte:
                    il.Emit(OpCodes.Conv_I1);
                    break;
                case NPTypeCode.Int16:
                    il.Emit(OpCodes.Conv_I2);
                    break;
                case NPTypeCode.UInt16:
                case NPTypeCode.Char:
                    il.Emit(OpCodes.Conv_U2);
                    break;
                case NPTypeCode.Int32:
                    il.Emit(OpCodes.Conv_I4);
                    break;
                case NPTypeCode.UInt32:
                    il.Emit(OpCodes.Conv_U4);
                    break;
                case NPTypeCode.Int64:
                    if (IsUnsigned(from))
                        il.Emit(OpCodes.Conv_U8);
                    else
                        il.Emit(OpCodes.Conv_I8);
                    break;
                case NPTypeCode.UInt64:
                    // signed source must SIGN-extend to 64 bits (NumPy: int16(-1)->uint64 =
                    // 0xFFFF...FFFF). Conv_U8 zero-extends, truncating the sign. (float->uint64
                    // was already routed through Converts above.)
                    il.Emit(IsUnsigned(from) ? OpCodes.Conv_U8 : OpCodes.Conv_I8);
                    break;
                case NPTypeCode.Single:
                    if (IsUnsigned(from))
                        il.Emit(OpCodes.Conv_R_Un);
                    il.Emit(OpCodes.Conv_R4);
                    break;
                case NPTypeCode.Double:
                    if (IsUnsigned(from))
                        il.Emit(OpCodes.Conv_R_Un);
                    il.Emit(OpCodes.Conv_R8);
                    break;
                default:
                    throw new NotSupportedException($"Conversion to {to} not supported");
            }
        }

        private static readonly System.Collections.Concurrent.ConcurrentDictionary<(NPTypeCode from, NPTypeCode to), MethodInfo> _convertsFloatToInt = new();

        /// <summary>
        /// Resolves the NumPy-faithful Converts.To{Int}({float|double}) method used to convert a
        /// float source to an integer/char destination (truncate toward zero, type-min sentinel on
        /// NaN/overflow, wrap to narrower widths). Cached per (from,to).
        /// </summary>
        private static MethodInfo ConvertsFloatToIntMethod(NPTypeCode from, NPTypeCode to)
            => _convertsFloatToInt.GetOrAdd((from, to), static k =>
            {
                Type srcType = k.from == NPTypeCode.Single ? typeof(float) : typeof(double);
                string name = "To" + k.to.ToString(); // ToByte, ToInt32, ..., ToChar
                return typeof(Converts).GetMethod(name, BindingFlags.Public | BindingFlags.Static, null, new[] { srcType }, null)
                    ?? throw new InvalidOperationException($"Converts.{name}({srcType.Name}) not found for {k.from}->{k.to}");
            });

        /// <summary>
        /// Emit decimal-specific conversions.
        /// </summary>
        private static void EmitDecimalConversion(ILGenerator il, NPTypeCode from, NPTypeCode to)
        {
            if (to == NPTypeCode.Decimal)
            {
                // Convert to decimal - need to handle bool/char first
                if (from == NPTypeCode.Boolean)
                {
                    // bool -> int -> decimal
                    il.Emit(OpCodes.Conv_I4);
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalImplicitFromInt, null);
                    return;
                }
                if (from == NPTypeCode.Char)
                {
                    // char -> int -> decimal
                    il.Emit(OpCodes.Conv_I4);
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalImplicitFromInt, null);
                    return;
                }
                if (from == NPTypeCode.Half)
                {
                    // Half -> double -> decimal
                    il.EmitCall(OpCodes.Call, CachedMethods.HalfToDouble, null);
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalExplicitFromDouble, null);
                    return;
                }
                if (from == NPTypeCode.Complex)
                {
                    // Complex -> Real (double) -> decimal (matches NumPy ComplexWarning truncation)
                    var realGetter = typeof(System.Numerics.Complex).GetProperty("Real")?.GetGetMethod()
                        ?? throw new InvalidOperationException("Complex.Real not found");
                    il.EmitCall(OpCodes.Call, realGetter, null);
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalExplicitFromDouble, null);
                    return;
                }

                var method = from switch
                {
                    NPTypeCode.Byte => CachedMethods.DecimalImplicitFromByte,
                    NPTypeCode.SByte => CachedMethods.DecimalImplicitFromSByte,
                    NPTypeCode.Int16 => CachedMethods.DecimalImplicitFromShort,
                    NPTypeCode.UInt16 => CachedMethods.DecimalImplicitFromUShort,
                    NPTypeCode.Int32 => CachedMethods.DecimalImplicitFromInt,
                    NPTypeCode.UInt32 => CachedMethods.DecimalImplicitFromUInt,
                    NPTypeCode.Int64 => CachedMethods.DecimalImplicitFromLong,
                    NPTypeCode.UInt64 => CachedMethods.DecimalImplicitFromULong,
                    NPTypeCode.Single => CachedMethods.DecimalExplicitFromFloat,
                    NPTypeCode.Double => CachedMethods.DecimalExplicitFromDouble,
                    _ => throw new NotSupportedException($"Cannot convert {from} to decimal")
                };
                il.EmitCall(OpCodes.Call, method, null);
            }
            else
            {
                // Convert from decimal - need to handle bool/char
                if (to == NPTypeCode.Boolean)
                {
                    // decimal -> int -> bool (compare with 0)
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalToInt32, null);
                    il.Emit(OpCodes.Ldc_I4_0);
                    il.Emit(OpCodes.Cgt_Un);
                    return;
                }
                if (to == NPTypeCode.Char)
                {
                    // decimal -> int -> char
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalToInt32, null);
                    il.Emit(OpCodes.Conv_U2);
                    return;
                }
                if (to == NPTypeCode.Half)
                {
                    // decimal -> double -> Half
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalToDouble, null);
                    il.EmitCall(OpCodes.Call, CachedMethods.DoubleToHalf, null);
                    return;
                }
                if (to == NPTypeCode.Complex)
                {
                    // decimal -> double -> Complex(d, 0)
                    il.EmitCall(OpCodes.Call, CachedMethods.DecimalToDouble, null);
                    il.Emit(OpCodes.Ldc_R8, 0.0);
                    var ctor = typeof(System.Numerics.Complex).GetConstructor(new[] { typeof(double), typeof(double) })
                        ?? throw new InvalidOperationException("Complex constructor not found");
                    il.Emit(OpCodes.Newobj, ctor);
                    return;
                }

                var method = to switch
                {
                    NPTypeCode.Byte => CachedMethods.DecimalToByte,
                    NPTypeCode.SByte => CachedMethods.DecimalToSByte,
                    NPTypeCode.Int16 => CachedMethods.DecimalToInt16,
                    NPTypeCode.UInt16 => CachedMethods.DecimalToUInt16,
                    NPTypeCode.Int32 => CachedMethods.DecimalToInt32,
                    NPTypeCode.UInt32 => CachedMethods.DecimalToUInt32,
                    NPTypeCode.Int64 => CachedMethods.DecimalToInt64,
                    NPTypeCode.UInt64 => CachedMethods.DecimalToUInt64,
                    NPTypeCode.Single => CachedMethods.DecimalToSingle,
                    NPTypeCode.Double => CachedMethods.DecimalToDouble,
                    _ => throw new NotSupportedException($"Cannot convert decimal to {to}")
                };
                il.EmitCall(OpCodes.Call, method, null);
            }
        }

        /// <summary>
        /// Emit Half or Complex type conversions (require method calls).
        /// </summary>
        private static void EmitHalfOrComplexConversion(ILGenerator il, NPTypeCode from, NPTypeCode to)
        {
            // Half -> other: convert Half to double first, then to target
            if (from == NPTypeCode.Half)
            {
                // Half.op_Explicit(Half) -> double (use cached method to avoid ambiguous match)
                il.EmitCall(OpCodes.Call, CachedMethods.HalfToDouble, null);

                if (to == NPTypeCode.Double)
                    return;  // Already double

                // Convert double to target type
                EmitConvertTo(il, NPTypeCode.Double, to);
                return;
            }

            // Complex -> other: get Real part as double, then convert
            if (from == NPTypeCode.Complex)
            {
                // Complex.Real property getter
                var realGetter = typeof(System.Numerics.Complex).GetProperty("Real")?.GetGetMethod()
                    ?? throw new InvalidOperationException("Complex.Real not found");
                il.EmitCall(OpCodes.Call, realGetter, null);

                if (to == NPTypeCode.Double)
                    return;  // Already double

                // Convert double to target type
                EmitConvertTo(il, NPTypeCode.Double, to);
                return;
            }

            // other -> Half: convert to double first, then to Half
            if (to == NPTypeCode.Half)
            {
                // First convert source to double
                if (from != NPTypeCode.Double && from != NPTypeCode.Single)
                    EmitConvertTo(il, from, NPTypeCode.Double);
                else if (from == NPTypeCode.Single)
                    il.Emit(OpCodes.Conv_R8);  // float to double

                // double -> Half via explicit cast (use cached method to avoid ambiguous match)
                il.EmitCall(OpCodes.Call, CachedMethods.DoubleToHalf, null);
                return;
            }

            // other -> Complex: convert to double, then create Complex with imaginary = 0
            if (to == NPTypeCode.Complex)
            {
                // First convert source to double
                if (from != NPTypeCode.Double && from != NPTypeCode.Single)
                    EmitConvertTo(il, from, NPTypeCode.Double);
                else if (from == NPTypeCode.Single)
                    il.Emit(OpCodes.Conv_R8);  // float to double

                // Load 0.0 for imaginary part
                il.Emit(OpCodes.Ldc_R8, 0.0);

                // new Complex(real, imaginary)
                var ctor = typeof(System.Numerics.Complex).GetConstructor(new[] { typeof(double), typeof(double) })
                    ?? throw new InvalidOperationException("Complex constructor not found");
                il.Emit(OpCodes.Newobj, ctor);
                return;
            }

            throw new NotSupportedException($"Conversion from {from} to {to} not supported");
        }

        /// <summary>
        /// Check if type is unsigned.
        /// </summary>
        internal static bool IsUnsigned(NPTypeCode type)
        {
            return type == NPTypeCode.Byte || type == NPTypeCode.UInt16 ||
                   type == NPTypeCode.UInt32 || type == NPTypeCode.UInt64 ||
                   type == NPTypeCode.Char;
        }

        /// <summary>
        /// Emit scalar operation for NPTypeCode.
        /// </summary>
        internal static void EmitScalarOperation(ILGenerator il, BinaryOp op, NPTypeCode resultType)
        {
            // Element-wise min/max (np.maximum / np.minimum / np.fmax / np.fmin): one scalar
            // clamp covers every dtype — the NaN-aware float/double/half/complex helpers,
            // Math.Max/Min for int + decimal, and a manual select for char. nanIgnore selects
            // the fmax/fmin (skip NaN) vs maximum/minimum (propagate NaN) policy. Intercept
            // BEFORE the decimal/half/complex routing so those dtypes also flow through the
            // shared clamp instead of EmitDecimal/Half/ComplexOperation (which have no min/max).
            if (op == BinaryOp.Maximum || op == BinaryOp.Minimum || op == BinaryOp.FMax || op == BinaryOp.FMin)
            {
                EmitScalarClamp(il, resultType,
                    isMax: op == BinaryOp.Maximum || op == BinaryOp.FMax,
                    nanIgnore: op == BinaryOp.FMax || op == BinaryOp.FMin);
                return;
            }

            // np.logaddexp / np.logaddexp2 / np.nextafter — one static helper per (op, loop dtype).
            // Operands arrive already in resultType (the generic binary loops converted both), and
            // each helper is signatured (resultType, resultType) -> resultType, so the whole scalar
            // op is a single Call. Intercept BEFORE the decimal/half routing (like min/max) so those
            // dtypes flow through the shared helpers rather than EmitDecimal/HalfOperation.
            if (op == BinaryOp.LogAddExp || op == BinaryOp.LogAddExp2 || op == BinaryOp.NextAfter || op == BinaryOp.CopySign)
            {
                EmitLogAddNextOperation(il, op, resultType);
                return;
            }

            // Special handling for decimal (uses operator methods)
            if (resultType == NPTypeCode.Decimal)
            {
                EmitDecimalOperation(il, op);
                return;
            }

            // Special handling for Half (uses operator methods)
            if (resultType == NPTypeCode.Half)
            {
                EmitHalfOperation(il, op);
                return;
            }

            // Special handling for Complex (uses operator methods)
            if (resultType == NPTypeCode.Complex)
            {
                EmitComplexOperation(il, op);
                return;
            }

            // Special handling for Power - requires Math.Pow call
            if (op == BinaryOp.Power)
            {
                EmitPowerOperation(il, resultType);
                return;
            }

            // Special handling for FloorDivide - division followed by floor for floats
            if (op == BinaryOp.FloorDivide)
            {
                EmitFloorDivideOperation(il, resultType);
                return;
            }

            // Special handling for Mod - NumPy uses floored division semantics (Python %)
            // Result sign matches divisor sign, not dividend sign (unlike C# %)
            if (op == BinaryOp.Mod)
            {
                EmitModOperation(il, resultType);
                return;
            }

            // Special handling for ATan2 - requires Math.Atan2 call
            if (op == BinaryOp.ATan2)
            {
                EmitATan2Operation(il, resultType);
                return;
            }

            // Special handling for LeftShift / RightShift - NumPy bit-shift semantics.
            // Stack on entry is [value(resultType), count(resultType)] (the generic binary
            // loops converted BOTH operands to resultType). EmitShiftFromStack consumes the
            // count as the shift amount (truncated to int for the IL shl/shr opcode) and
            // applies NumPy's overflow rule (count < 0 || count >= bitWidth) at the
            // resultType width — see DirectILKernelGenerator.Shift.cs.
            if (op == BinaryOp.LeftShift || op == BinaryOp.RightShift)
            {
                EmitShiftFromStack(il, resultType, op == BinaryOp.LeftShift);
                return;
            }

            // Boolean and/or/xor use NumPy's LOGICAL loops (loops.h.src aliases the bool bitwise
            // ufuncs onto logical_and / logical_or / not_equal): any nonzero byte is True, output
            // always canonical 0/1. Must mirror EmitBoolVectorLogicalOp exactly, or the SIMD body
            // and the scalar tail of one array would disagree on a non-canonical (frombuffer) bool.
            // Stack in: [a(u1-as-i4), b(u1-as-i4)].
            if (resultType == NPTypeCode.Boolean &&
                (op == BinaryOp.BitwiseAnd || op == BinaryOp.BitwiseOr || op == BinaryOp.BitwiseXor))
            {
                if (op == BinaryOp.BitwiseOr)
                {
                    // (a | b) != 0
                    il.Emit(OpCodes.Or);
                    il.Emit(OpCodes.Ldc_I4_0);
                    il.Emit(OpCodes.Cgt_Un);
                    return;
                }

                // (a != 0) and/xor (b != 0)
                var locB = il.DeclareLocal(typeof(int));
                il.Emit(OpCodes.Stloc, locB);       // [a]
                il.Emit(OpCodes.Ldc_I4_0);
                il.Emit(OpCodes.Cgt_Un);            // [na]
                il.Emit(OpCodes.Ldloc, locB);
                il.Emit(OpCodes.Ldc_I4_0);
                il.Emit(OpCodes.Cgt_Un);            // [na, nb]
                il.Emit(op == BinaryOp.BitwiseAnd ? OpCodes.And : OpCodes.Xor);
                return;
            }

            var opcode = op switch
            {
                BinaryOp.Add => OpCodes.Add,
                BinaryOp.Subtract => OpCodes.Sub,
                BinaryOp.Multiply => OpCodes.Mul,
                BinaryOp.Divide => IsUnsigned(resultType) ? OpCodes.Div_Un : OpCodes.Div,
                BinaryOp.BitwiseAnd => OpCodes.And,
                BinaryOp.BitwiseOr => OpCodes.Or,
                BinaryOp.BitwiseXor => OpCodes.Xor,
                _ => throw new NotSupportedException($"Operation {op} not supported")
            };

            il.Emit(opcode);
        }

        /// <summary>
        /// Emit Power operation.
        /// Stack: [base, exponent] -> [result]
        ///
        /// For integer result types, calls <see cref="Utilities.NDIntegerPower"/> to
        /// preserve dtype-native wrapping (matches NumPy's per-dtype integer power loop).
        /// For float32 result, uses MathF.Pow (single-precision, no f64 round-trip).
        /// For float64 result, uses Math.Pow.
        /// For any other result type that reaches this emit (e.g. Boolean), falls back to
        /// double round-trip to keep the kernel verifiable.
        /// </summary>
        private static void EmitPowerOperation(ILGenerator il, NPTypeCode resultType)
        {
            // Stack has: base (resultType), exponent (resultType)
            //
            // Integer result types use the NDIntegerPower helpers (squared-exp with
            // native wrapping). The helpers expect non-negative exponents; the caller
            // (DefaultEngine.Power) is responsible for the neg-exp ValueError pre-check.
            var intPow = GetIntegerPowMethod(resultType);
            if (intPow != null)
            {
                il.EmitCall(OpCodes.Call, intPow, null);
                return;
            }

            // float32 result → MathF.Pow (single-precision, matches NumPy powf)
            if (resultType == NPTypeCode.Single)
            {
                il.EmitCall(OpCodes.Call, CachedMethods.MathFPow, null);
                return;
            }

            // float64 result → Math.Pow directly (operands already double on stack).
            if (resultType == NPTypeCode.Double)
            {
                il.EmitCall(OpCodes.Call, CachedMethods.MathPow, null);
                return;
            }

            // Fallback (e.g. Boolean as result type): convert to double, call Math.Pow,
            // convert back. Keeps IL verification happy for non-numeric promotion targets.
            var locExp = il.DeclareLocal(GetClrType(resultType));
            il.Emit(OpCodes.Stloc, locExp);
            if (IsUnsigned(resultType))
                il.Emit(OpCodes.Conv_R_Un);
            il.Emit(OpCodes.Conv_R8);
            il.Emit(OpCodes.Ldloc, locExp);
            if (IsUnsigned(resultType))
                il.Emit(OpCodes.Conv_R_Un);
            il.Emit(OpCodes.Conv_R8);
            il.EmitCall(OpCodes.Call, CachedMethods.MathPow, null);
            EmitConvertFromDouble(il, resultType);
        }

        /// <summary>
        /// Return the integer-power helper MethodInfo for <paramref name="resultType"/>,
        /// or null when the result dtype is not integer (caller falls back to float power).
        /// </summary>
        private static MethodInfo? GetIntegerPowMethod(NPTypeCode resultType)
        {
            return resultType switch
            {
                NPTypeCode.SByte => CachedMethods.IntPowSByte,
                NPTypeCode.Byte => CachedMethods.IntPowByte,
                NPTypeCode.Int16 => CachedMethods.IntPowInt16,
                NPTypeCode.UInt16 => CachedMethods.IntPowUInt16,
                NPTypeCode.Char => CachedMethods.IntPowChar,
                NPTypeCode.Int32 => CachedMethods.IntPowInt32,
                NPTypeCode.UInt32 => CachedMethods.IntPowUInt32,
                NPTypeCode.Int64 => CachedMethods.IntPowInt64,
                NPTypeCode.UInt64 => CachedMethods.IntPowUInt64,
                _ => null
            };
        }

        /// <summary>
        /// Return the <see cref="Utilities.NDDivision"/> floor-division helper for
        /// <paramref name="resultType"/>, or null if the dtype routes elsewhere
        /// (Decimal/Half/Complex are handled before this is reached).
        /// </summary>
        private static MethodInfo? GetFloorDivideMethod(NPTypeCode resultType)
        {
            return resultType switch
            {
                NPTypeCode.SByte => CachedMethods.FloorDivSByte,
                NPTypeCode.Byte => CachedMethods.FloorDivByte,
                NPTypeCode.Int16 => CachedMethods.FloorDivInt16,
                NPTypeCode.UInt16 => CachedMethods.FloorDivUInt16,
                NPTypeCode.Char => CachedMethods.FloorDivChar,
                NPTypeCode.Int32 => CachedMethods.FloorDivInt32,
                NPTypeCode.UInt32 => CachedMethods.FloorDivUInt32,
                NPTypeCode.Int64 => CachedMethods.FloorDivInt64,
                NPTypeCode.UInt64 => CachedMethods.FloorDivUInt64,
                NPTypeCode.Single => CachedMethods.FloorDivSingle,
                NPTypeCode.Double => CachedMethods.FloorDivDouble,
                _ => null
            };
        }

        /// <summary>
        /// Return the <see cref="Utilities.NDDivision"/> remainder helper for
        /// <paramref name="resultType"/>, or null if the dtype routes elsewhere.
        /// </summary>
        private static MethodInfo? GetRemainderMethod(NPTypeCode resultType)
        {
            return resultType switch
            {
                NPTypeCode.SByte => CachedMethods.RemSByte,
                NPTypeCode.Byte => CachedMethods.RemByte,
                NPTypeCode.Int16 => CachedMethods.RemInt16,
                NPTypeCode.UInt16 => CachedMethods.RemUInt16,
                NPTypeCode.Char => CachedMethods.RemChar,
                NPTypeCode.Int32 => CachedMethods.RemInt32,
                NPTypeCode.UInt32 => CachedMethods.RemUInt32,
                NPTypeCode.Int64 => CachedMethods.RemInt64,
                NPTypeCode.UInt64 => CachedMethods.RemUInt64,
                NPTypeCode.Single => CachedMethods.RemSingle,
                NPTypeCode.Double => CachedMethods.RemDouble,
                _ => null
            };
        }

        /// <summary>
        /// Emit FloorDivide via the <see cref="Utilities.NDDivision"/> helpers, matching NumPy's
        /// <c>floor_div_@TYPE@</c> (integer ÷0 -> 0, signed floor toward -inf, MIN/-1 wrap) and
        /// <c>npy_floor_divide</c> (CPython divmod port: float ÷0 -> ±inf/nan, snap-to-nearest).
        /// Stack: [dividend, divisor] -> [result]
        /// </summary>
        private static void EmitFloorDivideOperation(ILGenerator il, NPTypeCode resultType)
        {
            var m = GetFloorDivideMethod(resultType);
            if (m != null)
            {
                il.EmitCall(OpCodes.Call, m, null);
                return;
            }
            // Boolean (or any other) result reaching here: fall back to plain integer division.
            il.Emit(IsUnsigned(resultType) ? OpCodes.Div_Un : OpCodes.Div);
        }

        /// <summary>
        /// Emit Mod via the <see cref="Utilities.NDDivision"/> helpers, matching NumPy's integer
        /// remainder (÷0 -> 0, floored sign) and <c>npy_remainder</c> (fmod-based CPython divmod).
        /// Stack: [dividend, divisor] -> [result]
        /// </summary>
        private static void EmitModOperation(ILGenerator il, NPTypeCode resultType)
        {
            var m = GetRemainderMethod(resultType);
            if (m != null)
            {
                il.EmitCall(OpCodes.Call, m, null);
                return;
            }
            // Boolean (or any other) result reaching here: fall back to plain remainder.
            il.Emit(IsUnsigned(resultType) ? OpCodes.Rem_Un : OpCodes.Rem);
        }

        /// <summary>
        /// Emit np.logaddexp / np.logaddexp2 / np.nextafter / np.copysign via the
        /// <see cref="Utilities.NDLogAddExpMath"/> scalar helpers. Stack: [x1, x2] (already in the loop
        /// dtype) -> [result]. Mirrors <see cref="EmitFloorDivideOperation"/>: resolve the per-dtype
        /// helper, then a single Call.
        /// </summary>
        private static void EmitLogAddNextOperation(ILGenerator il, BinaryOp op, NPTypeCode resultType)
            => il.EmitCall(OpCodes.Call, GetLogAddNextMethod(op, resultType), null);

        /// <summary>
        /// Return the <see cref="Utilities.NDLogAddExpMath"/> scalar helper for
        /// (<paramref name="op"/>, <paramref name="resultType"/>). Every loop dtype resolves to
        /// Half / Single / Double / Decimal (float-tier promotion, same as ATan2), so exactly one of
        /// four signatures applies per op. Parallels <see cref="GetFloorDivideMethod"/> /
        /// <see cref="GetRemainderMethod"/>, but the four ops are one family (intercepted together like
        /// maximum/minimum), so they share a single op-keyed resolver.
        /// </summary>
        private static MethodInfo GetLogAddNextMethod(BinaryOp op, NPTypeCode resultType)
        {
            switch (op)
            {
                case BinaryOp.LogAddExp:
                    return resultType switch
                    {
                        NPTypeCode.Half => CachedMethods.LogAddExpH,
                        NPTypeCode.Single => CachedMethods.LogAddExpF,
                        NPTypeCode.Decimal => CachedMethods.LogAddExpDec,
                        _ => CachedMethods.LogAddExpD,
                    };
                case BinaryOp.LogAddExp2:
                    return resultType switch
                    {
                        NPTypeCode.Half => CachedMethods.LogAddExp2H,
                        NPTypeCode.Single => CachedMethods.LogAddExp2F,
                        NPTypeCode.Decimal => CachedMethods.LogAddExp2Dec,
                        _ => CachedMethods.LogAddExp2D,
                    };
                case BinaryOp.NextAfter:
                    return resultType switch
                    {
                        NPTypeCode.Half => CachedMethods.NextAfterH,
                        NPTypeCode.Single => CachedMethods.NextAfterF,
                        NPTypeCode.Decimal => CachedMethods.NextAfterDec,
                        _ => CachedMethods.NextAfterD,
                    };
                default: // BinaryOp.CopySign
                    return resultType switch
                    {
                        NPTypeCode.Half => CachedMethods.CopySignH,
                        NPTypeCode.Single => CachedMethods.CopySignF,
                        NPTypeCode.Decimal => CachedMethods.CopySignDec,
                        _ => CachedMethods.CopySignD,
                    };
            }
        }

        /// <summary>
        /// Emit ATan2 operation using Math.Atan2.
        /// NumPy arctan2(y, x) computes the angle in radians between the positive x-axis
        /// and the point (x, y), with the correct quadrant determination.
        /// Stack: [y, x] -> [result] (angle in radians, range [-pi, pi])
        /// </summary>
        private static void EmitATan2Operation(ILGenerator il, NPTypeCode resultType)
        {
            // Math.Atan2(double y, double x) -> double
            // Stack has: y (resultType), x (resultType)
            // We need to convert both to double for Math.Atan2

            // Store x temporarily
            var locX = il.DeclareLocal(GetClrType(resultType));
            il.Emit(OpCodes.Stloc, locX);

            // Convert y to double
            if (resultType != NPTypeCode.Double)
            {
                if (IsUnsigned(resultType))
                    il.Emit(OpCodes.Conv_R_Un);
                il.Emit(OpCodes.Conv_R8);
            }

            // Load and convert x to double
            il.Emit(OpCodes.Ldloc, locX);
            if (resultType != NPTypeCode.Double)
            {
                if (IsUnsigned(resultType))
                    il.Emit(OpCodes.Conv_R_Un);
                il.Emit(OpCodes.Conv_R8);
            }

            // Call Math.Atan2(double y, double x)
            il.EmitCall(OpCodes.Call, CachedMethods.MathAtan2, null);

            // Convert result back to target type
            if (resultType != NPTypeCode.Double)
            {
                EmitConvertFromDouble(il, resultType);
            }
        }

        /// <summary>
        /// Emit conversion from double to target type.
        /// </summary>
        private static void EmitConvertFromDouble(ILGenerator il, NPTypeCode targetType)
        {
            switch (targetType)
            {
                case NPTypeCode.Byte:
                    il.Emit(OpCodes.Conv_U1);
                    break;
                case NPTypeCode.SByte:
                    il.Emit(OpCodes.Conv_I1);
                    break;
                case NPTypeCode.Int16:
                    il.Emit(OpCodes.Conv_I2);
                    break;
                case NPTypeCode.UInt16:
                case NPTypeCode.Char:
                    il.Emit(OpCodes.Conv_U2);
                    break;
                case NPTypeCode.Int32:
                    il.Emit(OpCodes.Conv_I4);
                    break;
                case NPTypeCode.UInt32:
                    il.Emit(OpCodes.Conv_U4);
                    break;
                case NPTypeCode.Int64:
                    il.Emit(OpCodes.Conv_I8);
                    break;
                case NPTypeCode.UInt64:
                    il.Emit(OpCodes.Conv_U8);
                    break;
                case NPTypeCode.Single:
                    il.Emit(OpCodes.Conv_R4);
                    break;
                case NPTypeCode.Boolean:
                    // double -> bool: != 0
                    il.Emit(OpCodes.Ldc_R8, 0.0);
                    il.Emit(OpCodes.Ceq);
                    il.Emit(OpCodes.Ldc_I4_0);
                    il.Emit(OpCodes.Ceq);
                    break;
                // NPTypeCode.Double needs no conversion
            }
        }

        /// <summary>
        /// Emit decimal-specific operation using operator methods.
        /// </summary>
        private static void EmitDecimalOperation(ILGenerator il, BinaryOp op)
        {
            // Bitwise operations not supported for decimal
            if (op == BinaryOp.BitwiseAnd || op == BinaryOp.BitwiseOr || op == BinaryOp.BitwiseXor)
                throw new NotSupportedException($"Bitwise operation {op} not supported for decimal type");

            // Power for decimal uses DecimalEx.Pow
            if (op == BinaryOp.Power)
            {
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalMathPow, null);
                return;
            }

            // FloorDivide for decimal: divide then floor toward negative infinity
            if (op == BinaryOp.FloorDivide)
            {
                // Stack: [dividend, divisor]
                // For NumPy semantics: floor(a/b), not truncate
                var locDivisor = il.DeclareLocal(typeof(decimal));
                var locDividend = il.DeclareLocal(typeof(decimal));
                il.Emit(OpCodes.Stloc, locDivisor);
                il.Emit(OpCodes.Stloc, locDividend);

                // Compute a / b
                il.Emit(OpCodes.Ldloc, locDividend);
                il.Emit(OpCodes.Ldloc, locDivisor);
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalOpDivision, null);

                // Call decimal.Floor for floored division toward negative infinity
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalFloor, null);
                return;
            }

            // Mod for decimal: NumPy floored modulo semantics
            // result = a - floor(a / b) * b
            if (op == BinaryOp.Mod)
            {
                // Stack: [dividend, divisor]
                var locDivisor = il.DeclareLocal(typeof(decimal));
                var locDividend = il.DeclareLocal(typeof(decimal));
                il.Emit(OpCodes.Stloc, locDivisor);
                il.Emit(OpCodes.Stloc, locDividend);

                // Load a for final subtraction
                il.Emit(OpCodes.Ldloc, locDividend);

                // Compute floor(a / b)
                il.Emit(OpCodes.Ldloc, locDividend);
                il.Emit(OpCodes.Ldloc, locDivisor);
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalOpDivision, null);
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalFloor, null);

                // Multiply by b
                il.Emit(OpCodes.Ldloc, locDivisor);
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalOpMultiply, null);

                // Subtract: a - floor(a/b)*b
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalOpSubtraction, null);
                return;
            }

            // ATan2 for decimal uses DecimalEx.ATan2
            if (op == BinaryOp.ATan2)
            {
                il.EmitCall(OpCodes.Call, CachedMethods.DecimalMathATan2, null);
                return;
            }

            var method = op switch
            {
                BinaryOp.Add => CachedMethods.DecimalOpAddition,
                BinaryOp.Subtract => CachedMethods.DecimalOpSubtraction,
                BinaryOp.Multiply => CachedMethods.DecimalOpMultiply,
                BinaryOp.Divide => CachedMethods.DecimalOpDivision,
                _ => throw new NotSupportedException($"Operation {op} not supported for decimal")
            };

            il.EmitCall(OpCodes.Call, method, null);
        }

        /// <summary>
        /// Emit Half-specific operation using operator methods.
        /// </summary>
        private static void EmitHalfOperation(ILGenerator il, BinaryOp op)
        {
            // Bitwise operations not supported for Half
            if (op == BinaryOp.BitwiseAnd || op == BinaryOp.BitwiseOr || op == BinaryOp.BitwiseXor)
                throw new NotSupportedException($"Bitwise operation {op} not supported for Half type");

            var halfToDouble = CachedMethods.HalfToDouble;

            // For all other operations, convert to double, perform operation, convert back
            // Stack: [half1, half2]
            var locRight = il.DeclareLocal(typeof(Half));
            il.Emit(OpCodes.Stloc, locRight);

            // Convert left to double
            il.EmitCall(OpCodes.Call, halfToDouble, null);

            // Convert right to double
            il.Emit(OpCodes.Ldloc, locRight);
            il.EmitCall(OpCodes.Call, halfToDouble, null);

            // Perform the operation in double
            switch (op)
            {
                case BinaryOp.Add:
                    il.Emit(OpCodes.Add);
                    break;
                case BinaryOp.Subtract:
                    il.Emit(OpCodes.Sub);
                    break;
                case BinaryOp.Multiply:
                    il.Emit(OpCodes.Mul);
                    break;
                case BinaryOp.Divide:
                    il.Emit(OpCodes.Div);
                    break;
                case BinaryOp.Power:
                    il.EmitCall(OpCodes.Call, CachedMethods.MathPow, null);
                    break;
                case BinaryOp.Mod:
                    // NumPy floored modulo: a - floor(a/b) * b
                    var locB = il.DeclareLocal(typeof(double));
                    var locA = il.DeclareLocal(typeof(double));
                    il.Emit(OpCodes.Stloc, locB);
                    il.Emit(OpCodes.Stloc, locA);
                    il.Emit(OpCodes.Ldloc, locA);
                    il.Emit(OpCodes.Ldloc, locA);
                    il.Emit(OpCodes.Ldloc, locB);
                    il.Emit(OpCodes.Div);
                    il.EmitCall(OpCodes.Call, CachedMethods.MathFloor, null);
                    il.Emit(OpCodes.Ldloc, locB);
                    il.Emit(OpCodes.Mul);
                    il.Emit(OpCodes.Sub);
                    break;
                case BinaryOp.FloorDivide:
                    // NumPy rule: floor_divide returns NaN when a/b is non-finite (inf or -inf).
                    // This matches numpy/core/src/umath/loops_arithmetic's npy_floor_divide_@type@.
                    il.Emit(OpCodes.Div);
                    EmitFloorWithInfToNaN(il);
                    break;
                case BinaryOp.ATan2:
                    il.EmitCall(OpCodes.Call, CachedMethods.MathAtan2, null);
                    break;
                default:
                    throw new NotSupportedException($"Operation {op} not supported for Half");
            }

            // Convert result back to Half (double -> Half)
            il.EmitCall(OpCodes.Call, CachedMethods.DoubleToHalf, null);
        }

        /// <summary>
        /// Emit Complex-specific operation using operator methods.
        /// </summary>
        private static void EmitComplexOperation(ILGenerator il, BinaryOp op)
        {
            // Bitwise operations not supported for Complex
            if (op == BinaryOp.BitwiseAnd || op == BinaryOp.BitwiseOr || op == BinaryOp.BitwiseXor)
                throw new NotSupportedException($"Bitwise operation {op} not supported for Complex type");

            // Complex has operator overloads we can call
            var complexType = typeof(System.Numerics.Complex);

            // Divide goes through a NumPy-compatible helper rather than the BCL's
            // op_Division: BCL's Smith's algorithm returns (NaN, NaN) for a/(0+0j),
            // whereas NumPy returns IEEE component-wise division (e.g. 1+0j -> inf+nanj).
            if (op == BinaryOp.Divide)
            {
                il.EmitCall(OpCodes.Call, GetHelper(nameof(ComplexDivideNumPy)), null);
                return;
            }

            var method = op switch
            {
                BinaryOp.Add => complexType.GetMethod("op_Addition", new[] { complexType, complexType }),
                BinaryOp.Subtract => complexType.GetMethod("op_Subtraction", new[] { complexType, complexType }),
                BinaryOp.Multiply => complexType.GetMethod("op_Multiply", new[] { complexType, complexType }),
                BinaryOp.Power => complexType.GetMethod("Pow", new[] { complexType, complexType }),
                _ => throw new NotSupportedException($"Operation {op} not supported for Complex")
            };

            if (method == null)
                throw new InvalidOperationException($"Could not find method for {op} on Complex");

            il.EmitCall(OpCodes.Call, method, null);
        }

        /// <summary>
        /// NumPy-compatible complex division. The .NET BCL's Complex.op_Division uses
        /// Smith's algorithm, which returns (NaN, NaN) when the divisor is (0+0j).
        /// NumPy instead produces IEEE component-wise division: (a.real/0, a.imag/0),
        /// giving (±inf, NaN) / (±inf, ±inf) / (NaN, NaN) depending on a's components.
        /// For all other cases we defer to the BCL operator — it's ULP-identical to
        /// NumPy for finite inputs.
        /// </summary>
        // Per-element op called from the binary IL kernel's inner loop — inline where possible,
        // full tier-1 codegen when invoked standalone via the kernel's Call.
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static System.Numerics.Complex ComplexDivideNumPy(System.Numerics.Complex a, System.Numerics.Complex b)
        {
            if (b.Real == 0.0 && b.Imaginary == 0.0)
                return new System.Numerics.Complex(a.Real / 0.0, a.Imaginary / 0.0);
            return a / b;
        }

        /// <summary>
        /// SIMD lane type for NPTypeCode. Boolean maps to <c>byte</c> — there is no
        /// <c>Vector{N}&lt;bool&gt;</c> in the BCL, and a NumSharp bool is stored as one byte, so
        /// bool kernels ride the byte lanes (NumPy's own BOOL loops do the same: its
        /// <c>loops_logical</c> SIMD operates on <c>u8</c> vectors). Every other type is its CLR
        /// type. Only the vector load/store/create helpers consult this; scalar loads keep
        /// <see cref="GetClrType"/>.
        /// </summary>
        internal static Type GetSimdLaneType(NPTypeCode type)
            => type == NPTypeCode.Boolean ? typeof(byte) : GetClrType(type);

        /// <summary>
        /// Emit Vector.Load for NPTypeCode (adapts to V128/V256/V512).
        /// Prefers <c>Avx.LoadVector256</c> / <c>Sse.LoadVector128</c> / <c>Avx512F.LoadVector512</c>
        /// when running on x86 — the JIT generates ~1.8x faster code than the cross-platform
        /// <c>Vector{N}.Load</c> path.
        /// </summary>
        internal static void EmitVectorLoad(ILGenerator il, NPTypeCode type)
        {
            var clrType = GetSimdLaneType(type);
            var x86 = VectorMethodCache.LoadX86(VectorBits, clrType);
            if (x86 != null)
            {
                il.EmitCall(OpCodes.Call, x86, null);
                return;
            }
            il.EmitCall(OpCodes.Call, VectorMethodCache.Load(VectorBits, clrType), null);
        }

        /// <summary>
        /// Emit Vector.Create for NPTypeCode (broadcasts scalar to all vector elements).
        /// Stack must have scalar value on top; result is Vector on stack.
        /// (Boolean broadcasts through the byte lanes — the bool scalar sits on the IL stack as
        /// its 0/1 byte value, which is exactly the <c>Vector.Create(byte)</c> argument.)
        /// </summary>
        internal static void EmitVectorCreate(ILGenerator il, NPTypeCode type)
            => il.EmitCall(OpCodes.Call, VectorMethodCache.CreateBroadcast(VectorBits, GetSimdLaneType(type)), null);

        /// <summary>
        /// Emit Vector.Store for NPTypeCode (adapts to V128/V256/V512).
        /// Stack convention: <c>..., Vector, T*</c> (matches <c>Vector{N}.Store(V, T*)</c>).
        /// When the x86 fast path applies, we swap the stack args because <c>Avx.Store</c>
        /// has reversed parameter order (<c>T*, V</c>) and use the X86 intrinsic instead.
        /// </summary>
        internal static void EmitVectorStore(ILGenerator il, NPTypeCode type)
        {
            var clrType = GetSimdLaneType(type);
            var x86 = VectorMethodCache.StoreX86(VectorBits, clrType);
            if (x86 != null)
            {
                // Stack: ..., V<T>, T*
                // Need:  ..., T*, V<T>  — swap via two locals.
                var vecType = VectorMethodCache.V(VectorBits, clrType);
                var locPtr = il.DeclareLocal(typeof(void*));
                var locVec = il.DeclareLocal(vecType);
                il.Emit(OpCodes.Stloc, locPtr);  // pop T*  → stack: ..., V
                il.Emit(OpCodes.Stloc, locVec);  // pop V   → stack: ...
                il.Emit(OpCodes.Ldloc, locPtr);  // push T* → stack: ..., T*
                il.Emit(OpCodes.Ldloc, locVec);  // push V  → stack: ..., T*, V
                il.EmitCall(OpCodes.Call, x86, null);
                return;
            }
            il.EmitCall(OpCodes.Call, VectorMethodCache.Store(VectorBits, clrType), null);
        }

        /// <summary>
        /// Emit Vector min/max (width-adaptive). Stack must hold two vectors;
        /// result is one vector. Prefers Avx/Avx2/Sse2 intrinsics on x86; falls back
        /// to the cross-platform <c>Vector{N}.Min/Max</c> for unsupported (op, T) — e.g.
        /// int64 on Avx2.
        /// <para>
        /// When <paramref name="propagateNaN"/> is set (and the element type is a float),
        /// the result propagates NaN per NumPy <c>maximum</c>/<c>minimum</c>/<c>clip</c>
        /// semantics: a lane is NaN whenever <em>either</em> operand is NaN. The hardware
        /// MAXPS/MAXPD (and <c>Avx.Max</c>/<c>Sse.Max</c>) instructions return the SECOND
        /// operand on an unordered (NaN) compare — i.e. <c>Max(NaN, finite)</c> yields the
        /// finite value, silently dropping the NaN. See the wrapper below for the fix.
        /// </para>
        /// </summary>
        internal static void EmitVectorMinOrMax(ILGenerator il, bool isMax, NPTypeCode type, bool propagateNaN = false)
        {
            string methodName = isMax ? "Max" : "Min";
            var clrType = GetClrType(type);

            // Integer lanes have no NaN, and callers that don't opt in keep the raw
            // hardware semantics — emit the bare min/max and return.
            bool nanAware = propagateNaN && (type == NPTypeCode.Single || type == NPTypeCode.Double);
            if (!nanAware)
            {
                EmitRawVectorMinOrMax(il, methodName, clrType);
                return;
            }

            // NaN-propagating float min/max (NumPy maximum/minimum/clip semantics).
            //
            //   Stack in : [a, b]      a = value vector, b = bound vector
            //   Stack out: [result]    result lane = isNaN(a) ? a : hwMinMax(a, b)
            //
            // hwMinMax returns the SECOND operand (b) on any NaN lane, so:
            //   * b is NaN   -> hwMinMax = b = NaN            (already propagates)
            //   * a is NaN   -> hwMinMax = b (finite)          (WRONG — restored below)
            //   * neither    -> hwMinMax = true min/max        (correct)
            // Selecting `a` back into exactly the a-is-NaN lanes restores propagation.
            // Robust even if the underlying min/max already propagates (then a is NaN
            // anyway, so the select is a no-op). Verified against Avx.Max/Min and the
            // cross-platform Vector{N}.Max/Min fallback.
            var vecType = VectorMethodCache.V(VectorBits, clrType);
            var locA = il.DeclareLocal(vecType);
            var locB = il.DeclareLocal(vecType);
            var locHw = il.DeclareLocal(vecType);

            il.Emit(OpCodes.Stloc, locB);                       // [a]
            il.Emit(OpCodes.Stloc, locA);                       // []

            // hw = hwMinMax(a, b)
            il.Emit(OpCodes.Ldloc, locA);
            il.Emit(OpCodes.Ldloc, locB);
            EmitRawVectorMinOrMax(il, methodName, clrType);     // [hw]
            il.Emit(OpCodes.Stloc, locHw);                      // []

            // equalMask = Equals(a, a) -> all-ones in non-NaN lanes, zero in NaN lanes
            // (ordered compare: NaN == NaN is false). This IS the lane-wise "not NaN" mask.
            il.Emit(OpCodes.Ldloc, locA);
            il.Emit(OpCodes.Ldloc, locA);
            il.EmitCall(OpCodes.Call, VectorMethodCache.Equals(VectorBits, clrType), null);       // [equalMask]

            // result = ConditionalSelect(equalMask, hw, a) -> non-NaN ? hw : NaN(=a)
            il.Emit(OpCodes.Ldloc, locHw);                     // [equalMask, hw]
            il.Emit(OpCodes.Ldloc, locA);                      // [equalMask, hw, a]
            il.EmitCall(OpCodes.Call, VectorMethodCache.ConditionalSelect(VectorBits, clrType), null); // [result]
        }

        // Bare width-adaptive Vector min/max: x86 intrinsic when available, else the
        // cross-platform Vector{N}.Min/Max. Consumes two vectors, leaves one. Does NOT
        // propagate NaN (hardware MAXPS/MINPD return the second operand on NaN).
        private static void EmitRawVectorMinOrMax(ILGenerator il, string methodName, Type clrType)
        {
            var x86 = VectorMethodCache.BinaryX86(VectorBits, methodName, clrType);
            if (x86 != null)
            {
                il.EmitCall(OpCodes.Call, x86, null);
                return;
            }
            // Generic fallback (e.g. ARM64). For float/double, Vector{N}.Max/Min lower to ARM
            // FMAX/FMINNM, whose signed-zero tie-break is IEEE maxNum/minNum (order-independent,
            // +0 > -0) — diverging from x86 MAXPS/MINPS (and NumPy maximum/minimum), which return
            // the SECOND operand on a ±0 tie. Emit an explicit strict compare + select so every
            // platform agrees: max = a > b ? a : b, min = a < b ? a : b (tie -> b). NaN lanes
            // also select b here; the EmitVectorMinOrMax wrapper restores NaN-from-a afterwards,
            // so NaN still propagates per NumPy. Integer min/max keep the bare hardware op.
            if (clrType == typeof(float) || clrType == typeof(double))
            {
                EmitFloatTieMinOrMax(il, methodName == "Max", clrType);
                return;
            }
            il.EmitCall(OpCodes.Call,
                VectorMethodCache.Generic(VectorBits, methodName, clrType, paramCount: 2), null);
        }

        // Cross-platform float/double min/max with MAXPS/MINPS tie semantics (second operand on
        // a ±0 / equal tie). Stack in: [a, b]; stack out: [result].
        private static void EmitFloatTieMinOrMax(ILGenerator il, bool isMax, Type clrType)
        {
            var vecType = VectorMethodCache.V(VectorBits, clrType);
            var locA = il.DeclareLocal(vecType);
            var locB = il.DeclareLocal(vecType);
            il.Emit(OpCodes.Stloc, locB);                       // [a]
            il.Emit(OpCodes.Stloc, locA);                       // []
            // mask = a (>|<) b  — strict; ties and NaN compares are false, so they select b.
            il.Emit(OpCodes.Ldloc, locA);
            il.Emit(OpCodes.Ldloc, locB);
            il.EmitCall(OpCodes.Call,
                isMax ? VectorMethodCache.GreaterThan(VectorBits, clrType)
                      : VectorMethodCache.LessThan(VectorBits, clrType), null);   // [mask]
            // result = ConditionalSelect(mask, a, b)
            il.Emit(OpCodes.Ldloc, locA);                       // [mask, a]
            il.Emit(OpCodes.Ldloc, locB);                       // [mask, a, b]
            il.EmitCall(OpCodes.Call, VectorMethodCache.ConditionalSelect(VectorBits, clrType), null); // [result]
        }

        // NaN-IGNORING vector min/max — np.fmax / np.fmin (ARRAY semantics).
        //
        //   Stack in : [a, b]
        //   Stack out: [result]   result = isNaN(a) ? b : (isNaN(b) ? a : hwMinMax(a, b))
        //
        // hwMinMax (raw MAXPS/MINPS, or the EmitFloatTieMinOrMax fallback) returns the SECOND
        // operand on a ±0 / equal tie — which is exactly what np.fmax/fmin's ARRAY ufunc does
        // (probed 2.4.2: fmax([+0],[-0]) -> -0, like maximum; NumPy's scalar fmax() differs but
        // the array loop does not). Only the NaN lanes need overriding: a NaN -> b, b NaN -> a,
        // both NaN -> b. Integer lanes have no NaN, so the bare hardware op is fmax/fmin too.
        internal static void EmitVectorFMinOrFMax(ILGenerator il, bool isMax, NPTypeCode type)
        {
            string methodName = isMax ? "Max" : "Min";
            var clrType = GetClrType(type);

            if (type != NPTypeCode.Single && type != NPTypeCode.Double)
            {
                EmitRawVectorMinOrMax(il, methodName, clrType);
                return;
            }

            var vecType = VectorMethodCache.V(VectorBits, clrType);
            var locA = il.DeclareLocal(vecType);
            var locB = il.DeclareLocal(vecType);
            var locHw = il.DeclareLocal(vecType);
            var locInner = il.DeclareLocal(vecType);

            il.Emit(OpCodes.Stloc, locB);                       // [a]
            il.Emit(OpCodes.Stloc, locA);                       // []

            // hw = hwMinMax(a, b)
            il.Emit(OpCodes.Ldloc, locA);
            il.Emit(OpCodes.Ldloc, locB);
            EmitRawVectorMinOrMax(il, methodName, clrType);     // [hw]
            il.Emit(OpCodes.Stloc, locHw);                      // []

            // inner = ConditionalSelect(Equals(a, a), hw, b)  -> a NaN ? b : hw
            il.Emit(OpCodes.Ldloc, locA);
            il.Emit(OpCodes.Ldloc, locA);
            il.EmitCall(OpCodes.Call, VectorMethodCache.Equals(VectorBits, clrType), null);   // [aEq]
            il.Emit(OpCodes.Ldloc, locHw);                     // [aEq, hw]
            il.Emit(OpCodes.Ldloc, locB);                      // [aEq, hw, b]
            il.EmitCall(OpCodes.Call, VectorMethodCache.ConditionalSelect(VectorBits, clrType), null); // [inner]
            il.Emit(OpCodes.Stloc, locInner);                  // []

            // result = ConditionalSelect(Equals(b, b), inner, a)  -> b NaN ? a : inner
            // (b-first ordering makes both-NaN return a, the FIRST operand, per the array loop.)
            il.Emit(OpCodes.Ldloc, locB);
            il.Emit(OpCodes.Ldloc, locB);
            il.EmitCall(OpCodes.Call, VectorMethodCache.Equals(VectorBits, clrType), null);   // [bEq]
            il.Emit(OpCodes.Ldloc, locInner);                  // [bEq, inner]
            il.Emit(OpCodes.Ldloc, locA);                      // [bEq, inner, a]
            il.EmitCall(OpCodes.Call, VectorMethodCache.ConditionalSelect(VectorBits, clrType), null); // [result]
        }

        /// <summary>
        /// Emit Vector operation for NPTypeCode (adapts to V128/V256/V512).
        /// Prefers x86 Avx/Avx2 intrinsics where available; falls back to operator overload.
        /// </summary>
        internal static void EmitVectorOperation(ILGenerator il, BinaryOp op, NPTypeCode type)
        {
            // Boolean rides the byte lanes with NumPy's LOGICAL semantics (loops.h.src aliases
            // BOOL_bitwise_and -> BOOL_logical_and, _or -> logical_or, _xor -> not_equal): any
            // nonzero byte counts as True and the output is always canonical 0/1 — a plain
            // bytewise AND would call 0x01 & 0x80 False where NumPy says True. Normalization is
            // one unsigned-byte Min against 1 (min(v,1): 0 -> 0, nonzero -> 1), NumPy's own
            // byte_to_true reshaped for pminub.
            if (type == NPTypeCode.Boolean &&
                (op == BinaryOp.BitwiseAnd || op == BinaryOp.BitwiseOr || op == BinaryOp.BitwiseXor))
            {
                EmitBoolVectorLogicalOp(il, op);
                return;
            }

            var clrType = GetClrType(type);

            if (op == BinaryOp.BitwiseAnd || op == BinaryOp.BitwiseOr || op == BinaryOp.BitwiseXor)
            {
                string methodName = op switch
                {
                    BinaryOp.BitwiseAnd => "BitwiseAnd",
                    BinaryOp.BitwiseOr => "BitwiseOr",
                    BinaryOp.BitwiseXor => "Xor",
                    _ => throw new NotSupportedException()
                };
                var x86Bit = VectorMethodCache.BinaryX86(VectorBits, methodName, clrType);
                if (x86Bit != null)
                {
                    il.EmitCall(OpCodes.Call, x86Bit, null);
                    return;
                }
                il.EmitCall(OpCodes.Call,
                    VectorMethodCache.Generic(VectorBits, methodName, clrType, paramCount: 2), null);
                return;
            }

            // Element-wise min/max SIMD (np.maximum / np.minimum / np.fmax / np.fmin).
            // maximum/minimum reuse the NaN-PROPAGATING clamp (EmitVectorMinOrMax) — the same
            // width-adaptive body np.clip uses, fuzz-validated: float/double propagate a NaN
            // operand (first-operand-NaN preferred, ±0 tie -> second), integers get the bare
            // hardware min/max. fmax/fmin use the NaN-IGNORING variant (skip a NaN operand).
            if (op == BinaryOp.Maximum || op == BinaryOp.Minimum)
            {
                EmitVectorMinOrMax(il, isMax: op == BinaryOp.Maximum, type, propagateNaN: true);
                return;
            }
            if (op == BinaryOp.FMax || op == BinaryOp.FMin)
            {
                EmitVectorFMinOrFMax(il, isMax: op == BinaryOp.FMax, type);
                return;
            }

            string arithName = op switch
            {
                BinaryOp.Add => "Add",
                BinaryOp.Subtract => "Subtract",
                BinaryOp.Multiply => "Multiply",
                BinaryOp.Divide => "Divide",
                _ => throw new NotSupportedException($"Operation {op} not supported for SIMD")
            };
            var x86Arith = VectorMethodCache.BinaryX86(VectorBits, arithName, clrType);
            if (x86Arith != null)
            {
                il.EmitCall(OpCodes.Call, x86Arith, null);
                return;
            }
            // Arithmetic uses operator overloads on Vector256<T> / Vector128<T>.
            string operatorName = op switch
            {
                BinaryOp.Add => "op_Addition",
                BinaryOp.Subtract => "op_Subtraction",
                BinaryOp.Multiply => "op_Multiply",
                BinaryOp.Divide => "op_Division",
                _ => throw new NotSupportedException($"Operation {op} not supported for SIMD")
            };
            il.EmitCall(OpCodes.Call, VectorMethodCache.Operator(VectorBits, clrType, operatorName), null);
        }

        /// <summary>
        /// Boolean and/or/xor over byte lanes, NumPy-normalized (see <see cref="EmitVectorOperation"/>).
        /// Stack: [a(V&lt;byte&gt;), b(V&lt;byte&gt;)] -> [result(V&lt;byte&gt;)], every result lane 0 or 1.
        ///   and: min(min(a,b), 1)      — min(a,b) nonzero iff both nonzero (NumPy simd_logical_and_u8)
        ///   or : min(a|b, 1)           — a|b nonzero iff either nonzero (NumPy simd_logical_or_u8)
        ///   xor: min(a,1) ^ min(b,1)   — normalize each, then xor (NumPy BOOL_not_equal's cmpeq+xor)
        /// </summary>
        private static void EmitBoolVectorLogicalOp(ILGenerator il, BinaryOp op)
        {
            var laneType = typeof(byte);

            void PushOnes()
            {
                il.Emit(OpCodes.Ldc_I4_1);
                il.EmitCall(OpCodes.Call, VectorMethodCache.CreateBroadcast(VectorBits, laneType), null);
            }

            switch (op)
            {
                case BinaryOp.BitwiseAnd:
                    EmitRawVectorMinOrMax(il, "Min", laneType);   // [min(a,b)]
                    PushOnes();                                   // [min, 1]
                    EmitRawVectorMinOrMax(il, "Min", laneType);   // [r]
                    return;

                case BinaryOp.BitwiseOr:
                {
                    var or = VectorMethodCache.BinaryX86(VectorBits, "BitwiseOr", laneType)
                             ?? VectorMethodCache.Generic(VectorBits, "BitwiseOr", laneType, paramCount: 2);
                    il.EmitCall(OpCodes.Call, or, null);          // [a|b]
                    PushOnes();                                   // [a|b, 1]
                    EmitRawVectorMinOrMax(il, "Min", laneType);   // [r]
                    return;
                }

                case BinaryOp.BitwiseXor:
                {
                    var xor = VectorMethodCache.BinaryX86(VectorBits, "Xor", laneType)
                              ?? VectorMethodCache.Generic(VectorBits, "Xor", laneType, paramCount: 2);
                    var locNb = il.DeclareLocal(VectorMethodCache.V(VectorBits, laneType));
                    PushOnes();                                   // [a, b, 1]
                    EmitRawVectorMinOrMax(il, "Min", laneType);   // [a, nb]
                    il.Emit(OpCodes.Stloc, locNb);                // [a]
                    PushOnes();                                   // [a, 1]
                    EmitRawVectorMinOrMax(il, "Min", laneType);   // [na]
                    il.Emit(OpCodes.Ldloc, locNb);                // [na, nb]
                    il.EmitCall(OpCodes.Call, xor, null);         // [r]
                    return;
                }

                default:
                    throw new NotSupportedException($"Operation {op} not supported for Boolean SIMD");
            }
        }

        #endregion

    }
}