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));

            /// <summary><see cref="MathF.Atan2(float, float)"/>: the CRT's <c>atan2f</c>, NumPy's float32 arctan2 loop.</summary>
            public static readonly MethodInfo MathFAtan2 = typeof(MathF).GetMethod(nameof(MathF.Atan2), new[] { typeof(float), typeof(float) })
                ?? throw new MissingMethodException(typeof(MathF).FullName, nameof(MathF.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));

            // np.spacing scalar kernels — the UNARY IEEE-step sibling of nextafter (one loop-dtype
            // param, so resolved by a single-type GetMethod rather than the binary LA above). All four
            // overloads are named Spacing; the correct one is bound by its parameter type. Backs the
            // SIMD kernel's scalar tail (float32/float64) and the strided / Half / Decimal paths.
            private static MethodInfo SP(Type t) =>
                typeof(Utilities.NDSpacingMath).GetMethod(nameof(Utilities.NDSpacingMath.Spacing), new[] { t })
                ?? throw new MissingMethodException(typeof(Utilities.NDSpacingMath).FullName, nameof(Utilities.NDSpacingMath.Spacing));
            public static readonly MethodInfo SpacingD = SP(typeof(double));
            public static readonly MethodInfo SpacingF = SP(typeof(float));
            public static readonly MethodInfo SpacingH = SP(typeof(Half));
            public static readonly MethodInfo SpacingDec = SP(typeof(decimal));

            // np.hypot scalar kernels — same shape/family as the LogAddNext helpers above, resolved by
            // GetLogAddNextMethod. UCRT's sqrt(x*x + y*y) hypot (NumPy's npy_hypot) lives in NDHypotMath.
            private static MethodInfo HY(string name, Type t) =>
                typeof(Utilities.NDHypotMath).GetMethod(name, new[] { t, t })
                ?? throw new MissingMethodException(typeof(Utilities.NDHypotMath).FullName, name);
            public static readonly MethodInfo HypotD = HY(nameof(Utilities.NDHypotMath.Hypot), typeof(double));
            public static readonly MethodInfo HypotF = HY(nameof(Utilities.NDHypotMath.HypotF), typeof(float));
            public static readonly MethodInfo HypotH = HY(nameof(Utilities.NDHypotMath.HypotHalf), typeof(Half));
            public static readonly MethodInfo HypotDec = HY(nameof(Utilities.NDHypotMath.HypotDecimal), typeof(decimal));

            // np.heaviside scalar kernels — same (op, loop dtype) -> one 2-arg Call family shape as the
            // LogAddNext/hypot helpers above, but backed by NDHeavisideMath (the step function). Used by
            // the strided / mixed-dtype / f16 / decimal IL path; the contiguous & scalar-broadcast f32/f64
            // cases take NDHeavisideMath's SIMD driver in Default.LogAddExp.cs instead.
            private static MethodInfo HV(string name, Type t) =>
                typeof(Utilities.NDHeavisideMath).GetMethod(name, new[] { t, t })
                ?? throw new MissingMethodException(typeof(Utilities.NDHeavisideMath).FullName, name);
            public static readonly MethodInfo HeavisideD = HV(nameof(Utilities.NDHeavisideMath.Heaviside), typeof(double));
            public static readonly MethodInfo HeavisideF = HV(nameof(Utilities.NDHeavisideMath.HeavisideF), typeof(float));
            public static readonly MethodInfo HeavisideH = HV(nameof(Utilities.NDHeavisideMath.HeavisideHalf), typeof(Half));
            public static readonly MethodInfo HeavisideDec = HV(nameof(Utilities.NDHeavisideMath.HeavisideDecimal), 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));

            // C-style fmod helpers (NDDivision, np.fmod) — truncated remainder, sign of dividend.
            // Same 2-arg signature as the Rem* family, resolved via NDDiv.
            public static readonly MethodInfo FmodSByte = NDDiv(nameof(Utilities.NDDivision.FmodSByte), typeof(sbyte));
            public static readonly MethodInfo FmodByte = NDDiv(nameof(Utilities.NDDivision.FmodByte), typeof(byte));
            public static readonly MethodInfo FmodInt16 = NDDiv(nameof(Utilities.NDDivision.FmodInt16), typeof(short));
            public static readonly MethodInfo FmodUInt16 = NDDiv(nameof(Utilities.NDDivision.FmodUInt16), typeof(ushort));
            public static readonly MethodInfo FmodChar = NDDiv(nameof(Utilities.NDDivision.FmodChar), typeof(char));
            public static readonly MethodInfo FmodInt32 = NDDiv(nameof(Utilities.NDDivision.FmodInt32), typeof(int));
            public static readonly MethodInfo FmodUInt32 = NDDiv(nameof(Utilities.NDDivision.FmodUInt32), typeof(uint));
            public static readonly MethodInfo FmodInt64 = NDDiv(nameof(Utilities.NDDivision.FmodInt64), typeof(long));
            public static readonly MethodInfo FmodUInt64 = NDDiv(nameof(Utilities.NDDivision.FmodUInt64), typeof(ulong));
            public static readonly MethodInfo FmodSingle = NDDiv(nameof(Utilities.NDDivision.FmodSingle), typeof(float));
            public static readonly MethodInfo FmodDouble = NDDiv(nameof(Utilities.NDDivision.FmodDouble), typeof(double));
            public static readonly MethodInfo FmodDecimal = NDDiv(nameof(Utilities.NDDivision.FmodDecimal), typeof(decimal));

            // gcd/lcm helpers (NDGcdLcm, np.gcd / np.lcm) — integer-only, same 2-arg (T,T)->T shape as
            // the Rem*/Fmod* families but resolved from NDGcdLcm. NO float/decimal/complex entries: those
            // dtypes have no gcd/lcm loop and are rejected at the np.* boundary before a kernel is built.
            private static MethodInfo NDGcd(string name, Type t) =>
                typeof(Utilities.NDGcdLcm).GetMethod(name, new[] { t, t })
                ?? throw new MissingMethodException(typeof(Utilities.NDGcdLcm).FullName, name);

            public static readonly MethodInfo GcdSByte = NDGcd(nameof(Utilities.NDGcdLcm.GcdSByte), typeof(sbyte));
            public static readonly MethodInfo GcdByte = NDGcd(nameof(Utilities.NDGcdLcm.GcdByte), typeof(byte));
            public static readonly MethodInfo GcdInt16 = NDGcd(nameof(Utilities.NDGcdLcm.GcdInt16), typeof(short));
            public static readonly MethodInfo GcdUInt16 = NDGcd(nameof(Utilities.NDGcdLcm.GcdUInt16), typeof(ushort));
            public static readonly MethodInfo GcdChar = NDGcd(nameof(Utilities.NDGcdLcm.GcdChar), typeof(char));
            public static readonly MethodInfo GcdInt32 = NDGcd(nameof(Utilities.NDGcdLcm.GcdInt32), typeof(int));
            public static readonly MethodInfo GcdUInt32 = NDGcd(nameof(Utilities.NDGcdLcm.GcdUInt32), typeof(uint));
            public static readonly MethodInfo GcdInt64 = NDGcd(nameof(Utilities.NDGcdLcm.GcdInt64), typeof(long));
            public static readonly MethodInfo GcdUInt64 = NDGcd(nameof(Utilities.NDGcdLcm.GcdUInt64), typeof(ulong));

            public static readonly MethodInfo LcmSByte = NDGcd(nameof(Utilities.NDGcdLcm.LcmSByte), typeof(sbyte));
            public static readonly MethodInfo LcmByte = NDGcd(nameof(Utilities.NDGcdLcm.LcmByte), typeof(byte));
            public static readonly MethodInfo LcmInt16 = NDGcd(nameof(Utilities.NDGcdLcm.LcmInt16), typeof(short));
            public static readonly MethodInfo LcmUInt16 = NDGcd(nameof(Utilities.NDGcdLcm.LcmUInt16), typeof(ushort));
            public static readonly MethodInfo LcmChar = NDGcd(nameof(Utilities.NDGcdLcm.LcmChar), typeof(char));
            public static readonly MethodInfo LcmInt32 = NDGcd(nameof(Utilities.NDGcdLcm.LcmInt32), typeof(int));
            public static readonly MethodInfo LcmUInt32 = NDGcd(nameof(Utilities.NDGcdLcm.LcmUInt32), typeof(uint));
            public static readonly MethodInfo LcmInt64 = NDGcd(nameof(Utilities.NDGcdLcm.LcmInt64), typeof(long));
            public static readonly MethodInfo LcmUInt64 = NDGcd(nameof(Utilities.NDGcdLcm.LcmUInt64), typeof(ulong));

            // Fused divmod helpers (NDDivision, np.divmod) — floored quotient (return) + floored
            // remainder (out). Signature is (T, T, out T) -> T, so resolve with the by-ref third arg.
            private static MethodInfo NDDiv3(string name, Type t) =>
                typeof(Utilities.NDDivision).GetMethod(name, new[] { t, t, t.MakeByRefType() })
                ?? throw new MissingMethodException(typeof(Utilities.NDDivision).FullName, name);

            public static readonly MethodInfo DivmodSByte = NDDiv3(nameof(Utilities.NDDivision.DivmodSByte), typeof(sbyte));
            public static readonly MethodInfo DivmodByte = NDDiv3(nameof(Utilities.NDDivision.DivmodByte), typeof(byte));
            public static readonly MethodInfo DivmodInt16 = NDDiv3(nameof(Utilities.NDDivision.DivmodInt16), typeof(short));
            public static readonly MethodInfo DivmodUInt16 = NDDiv3(nameof(Utilities.NDDivision.DivmodUInt16), typeof(ushort));
            public static readonly MethodInfo DivmodChar = NDDiv3(nameof(Utilities.NDDivision.DivmodChar), typeof(char));
            public static readonly MethodInfo DivmodInt32 = NDDiv3(nameof(Utilities.NDDivision.DivmodInt32), typeof(int));
            public static readonly MethodInfo DivmodUInt32 = NDDiv3(nameof(Utilities.NDDivision.DivmodUInt32), typeof(uint));
            public static readonly MethodInfo DivmodInt64 = NDDiv3(nameof(Utilities.NDDivision.DivmodInt64), typeof(long));
            public static readonly MethodInfo DivmodUInt64 = NDDiv3(nameof(Utilities.NDDivision.DivmodUInt64), typeof(ulong));
            public static readonly MethodInfo DivmodHalf = NDDiv3(nameof(Utilities.NDDivision.DivmodHalf), typeof(Half));
            public static readonly MethodInfo DivmodSingle = NDDiv3(nameof(Utilities.NDDivision.DivmodSingle), typeof(float));
            public static readonly MethodInfo DivmodDouble = NDDiv3(nameof(Utilities.NDDivision.DivmodDouble), typeof(double));
            public static readonly MethodInfo DivmodDecimal = NDDiv3(nameof(Utilities.NDDivision.DivmodDecimal), typeof(decimal));

            // 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));

            // NumPy-exact float16 / uint64 conversions (Converts.NumPyFloatRules.cs) for the conversion emitter
            // (EmitConvertTo / EmitHalfOrComplexConversion) — which every IL kernel uses to convert an operand to its
            // loop dtype (the scalar cast kernels, mixed-dtype binary ops, np.evaluate, polynomial evaluation,
            // reductions) — and for the uint64 -> float64 operand widenings. NumPy converts float16 bit by bit — a
            // signalling NaN stays signalling, where the BCL casts above set the quiet bit, so maximum(f16, f32) returned
            // a quieted NaN — and rounds uint64 -> float32 / float64 once, where conv.r.un + conv.r4 / conv.r8 round
            // twice (.NET 8). The float16 ARITHMETIC paths (widen, compute, narrow) keep the BCL casts above: the
            // arithmetic quiets a NaN anyway.
            public static readonly MethodInfo ConvertsHalfToSingle = ConvertsOverload(nameof(Converts.ToSingle), typeof(Half));
            public static readonly MethodInfo ConvertsHalfToDouble = ConvertsOverload(nameof(Converts.ToDouble), typeof(Half));
            public static readonly MethodInfo ConvertsSingleToHalf = ConvertsOverload(nameof(Converts.ToHalf), typeof(float));
            public static readonly MethodInfo ConvertsDoubleToHalf = ConvertsOverload(nameof(Converts.ToHalf), typeof(double));
            public static readonly MethodInfo ConvertsUInt64ToSingle = ConvertsOverload(nameof(Converts.ToSingle), typeof(ulong));
            public static readonly MethodInfo ConvertsUInt64ToDouble = ConvertsOverload(nameof(Converts.ToDouble), typeof(ulong));

            // complex128 -> float64 as NumPy's unsafe cast does it (the real part; NumPy only warns ComplexWarning). A
            // STATIC method that takes the struct BY VALUE, which is exactly what the conversion emitter has on the IL
            // stack — the instance getter Complex.Real needs a managed reference as 'this' (ComplexGetReal above is
            // called after an ldloca), and calling it on the value made the JIT read the real part's bits as an
            // address: every emitted complex -> real / decimal conversion crashed the process (0x80131506).
            public static readonly MethodInfo ConvertsComplexToDouble = ConvertsOverload(nameof(Converts.ToDouble), typeof(System.Numerics.Complex));

            // The PACKED float32 <-> float64 conversions for a scalar value (see EmitConvertTo's Single/Double arms):
            // cvtps2pd / cvtpd2ps write their whole destination register, where the scalar cvtss2sd / cvtsd2ss that
            // conv.r8 / conv.r4 compile to MERGE into it — a false dependency on the register's previous value that
            // the JIT does not break, which serializes an element loop on the conversion's latency (a 100K
            // float32 -> float64 add loop: 128 us scalar, 37 us packed). Same IEEE conversion, same NaN quieting.
            public static readonly MethodInfo Vector128CreateScalarUnsafeSingle =
                typeof(Vector128).GetMethod(nameof(Vector128.CreateScalarUnsafe), new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(Vector128).FullName, "CreateScalarUnsafe(float)");
            public static readonly MethodInfo Vector128CreateScalarUnsafeDouble =
                typeof(Vector128).GetMethod(nameof(Vector128.CreateScalarUnsafe), new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(Vector128).FullName, "CreateScalarUnsafe(double)");
            public static readonly MethodInfo Sse2ConvertSingleToDouble =
                typeof(System.Runtime.Intrinsics.X86.Sse2).GetMethod(nameof(System.Runtime.Intrinsics.X86.Sse2.ConvertToVector128Double), new[] { typeof(Vector128<float>) })
                ?? throw new MissingMethodException("Sse2", "ConvertToVector128Double(Vector128<float>)");
            public static readonly MethodInfo Sse2ConvertDoubleToSingle =
                typeof(System.Runtime.Intrinsics.X86.Sse2).GetMethod(nameof(System.Runtime.Intrinsics.X86.Sse2.ConvertToVector128Single), new[] { typeof(Vector128<double>) })
                ?? throw new MissingMethodException("Sse2", "ConvertToVector128Single(Vector128<double>)");
            public static readonly MethodInfo Vector128ToScalarDouble = Vector128ToScalar(typeof(double));
            public static readonly MethodInfo Vector128ToScalarSingle = Vector128ToScalar(typeof(float));

            /// <summary>Resolves <c>Vector128.ToScalar&lt;T&gt;(Vector128&lt;T&gt;)</c> closed over <paramref name="element"/>.</summary>
            /// <param name="element">The element type.</param>
            /// <returns>The closed generic method.</returns>
            /// <exception cref="MissingMethodException">The BCL has no such method (fails the type's initialization, i.e. fast).</exception>
            private static MethodInfo Vector128ToScalar(Type element) =>
                (typeof(Vector128).GetMethods(BindingFlags.Public | BindingFlags.Static)
                    .FirstOrDefault(m => m.Name == nameof(Vector128.ToScalar) && m.IsGenericMethodDefinition && m.GetParameters().Length == 1)
                 ?? throw new MissingMethodException(typeof(Vector128).FullName, "ToScalar<T>"))
                .MakeGenericMethod(element);

            /// <summary>Resolves the public one-parameter <c>Converts.{name}({parameter})</c> overload.</summary>
            /// <param name="name">The method name (ToSingle, ToDouble, ToHalf, ...).</param>
            /// <param name="parameter">The single parameter's type.</param>
            /// <returns>The overload.</returns>
            /// <exception cref="MissingMethodException">No such overload exists (fails the type's initialization, i.e. fast).</exception>
            private static MethodInfo ConvertsOverload(string name, Type parameter) =>
                typeof(Converts).GetMethod(name, BindingFlags.Public | BindingFlags.Static, null, new[] { parameter }, null)
                ?? throw new MissingMethodException(typeof(Converts).FullName, $"{name}({parameter.Name})");
            // 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 rather than Complex.Abs: the helper is a port
            // of NumPy's SIMD simd_cabsolute, whose sqrt(FMA(r, r, 1))·larger is FUSED on every current
            // dispatch target — Complex.Abs's unfused form differed from numpy.abs on 35.5% of random
            // inputs — and it keeps C99 hypot's inf-beats-NaN rule, where the BCL's private Hypot
            // returns NaN for abs(NaN+inf*i) on net8.0 and NumPy returns +inf.
            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");
            // float32 / float64 log1p and expm1: NDLog1pExpm1Math, the models of the UCRT routines NumPy's win-amd64 loops
            // call (log1p / log1pf bit-exact; expm1 / expm1f bit-exact outside [-ln 2, ln 1.5], Sun's s_expm1 inside). The
            // Exp(x) - 1 / Log(1 + x) compositions they replace cancel near zero. Scalar only, as NumPy's loops are there.

            /// <summary>
            /// <see cref="NumSharp.Utilities.NDLog1pExpm1Math.Log1p(double)"/>: float64 <c>np.log1p</c>'s element function,
            /// NumPy's bits (the UCRT model); emitted by the scalar unary path, which every float64 log1p kernel takes.
            /// </summary>
            public static readonly MethodInfo DoubleLog1p = typeof(NumSharp.Utilities.NDLog1pExpm1Math).GetMethod("Log1p", BindingFlags.Public | BindingFlags.Static, new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDLog1pExpm1Math).FullName, "Log1p");

            /// <summary>
            /// <see cref="NumSharp.Utilities.NDLog1pExpm1Math.Log1pF(float)"/>: float32 <c>np.log1p</c>'s element function,
            /// the UCRT's <c>log1pf</c> bit for bit (float32 arithmetic throughout, not a narrowed double).
            /// </summary>
            public static readonly MethodInfo SingleLog1p = typeof(NumSharp.Utilities.NDLog1pExpm1Math).GetMethod("Log1pF", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDLog1pExpm1Math).FullName, "Log1pF");

            /// <summary>
            /// <see cref="NumSharp.Utilities.NDLog1pExpm1Math.Expm1(double)"/>: float64 <c>np.expm1</c>'s element function —
            /// NumPy's bits outside <c>[-ln 2, ln 1.5]</c>, within 2 ULP of them inside (Sun's <c>s_expm1</c>).
            /// </summary>
            public static readonly MethodInfo DoubleExpm1 = typeof(NumSharp.Utilities.NDLog1pExpm1Math).GetMethod("Expm1", BindingFlags.Public | BindingFlags.Static, new[] { typeof(double) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDLog1pExpm1Math).FullName, "Expm1");

            /// <summary>
            /// <see cref="NumSharp.Utilities.NDLog1pExpm1Math.Expm1F(float)"/>: float32 <c>np.expm1</c>'s element function —
            /// the UCRT's <c>expf(x) - 1f</c> outside the float32 band, the rounded float64 Sun port (within 3 ULP of
            /// NumPy's undocumented polynomial) inside it.
            /// </summary>
            public static readonly MethodInfo SingleExpm1 = typeof(NumSharp.Utilities.NDLog1pExpm1Math).GetMethod("Expm1F", BindingFlags.Public | BindingFlags.Static, new[] { typeof(float) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDLog1pExpm1Math).FullName, "Expm1F");
            // 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");
            // Multiply routes through NDComplexMath (FMA-contracted a*b == NumPy's simd_cmul): the BCL
            // Complex.op_Multiply is the NAIVE a_re*b_re - a_im*b_im, which diverges from NumPy's fused
            // vfmaddsub on ~14% of operands (up to ~2840 ULP in the cancellation regime).
            public static readonly MethodInfo ComplexMultiply = typeof(NumSharp.Utilities.NDComplexMath).GetMethod("Multiply", BindingFlags.Public | BindingFlags.Static, new[] { typeof(System.Numerics.Complex), typeof(System.Numerics.Complex) })
                ?? throw new MissingMethodException(typeof(NumSharp.Utilities.NDComplexMath).FullName, "Multiply");
            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:
                    // uint64 -> float32 rounds ONCE, like NumPy's C cast: conv.r.un + conv.r4 go through float64 and
                    // round twice for many values past 2^53 (.NET 8's own ulong -> float does the same).
                    if (from == NPTypeCode.UInt64)
                    {
                        il.EmitCall(OpCodes.Call, CachedMethods.ConvertsUInt64ToSingle, null);
                        break;
                    }
                    // float64 -> float32 as the packed cvtpd2ps (see CachedMethods: no false dependency on the target
                    // register, same round-to-nearest-even and NaN quieting as the scalar cvtsd2ss conv.r4 compiles to).
                    if (from == NPTypeCode.Double && System.Runtime.Intrinsics.X86.Sse2.IsSupported)
                    {
                        il.EmitCall(OpCodes.Call, CachedMethods.Vector128CreateScalarUnsafeDouble, null);
                        il.EmitCall(OpCodes.Call, CachedMethods.Sse2ConvertDoubleToSingle, null);
                        il.EmitCall(OpCodes.Call, CachedMethods.Vector128ToScalarSingle, null);
                        break;
                    }
                    if (IsUnsigned(from))
                        il.Emit(OpCodes.Conv_R_Un);
                    il.Emit(OpCodes.Conv_R4);
                    break;
                case NPTypeCode.Double:
                    // uint64 -> float64 rounds ONCE, like NumPy's C cast: .NET 8 converts a value >= 2^63 as signed
                    // and adds 2^64, which rounds twice.
                    if (from == NPTypeCode.UInt64)
                    {
                        il.EmitCall(OpCodes.Call, CachedMethods.ConvertsUInt64ToDouble, null);
                        break;
                    }
                    // float32 -> float64 as the packed cvtps2pd (see CachedMethods: the scalar cvtss2sd that conv.r8
                    // compiles to merges into its target register, and that false dependency made a float32 -> float64
                    // element loop latency-bound — 3.5x slower). Exact either way; a signalling NaN is quieted either way.
                    if (from == NPTypeCode.Single && System.Runtime.Intrinsics.X86.Sse2.IsSupported)
                    {
                        il.EmitCall(OpCodes.Call, CachedMethods.Vector128CreateScalarUnsafeSingle, null);
                        il.EmitCall(OpCodes.Call, CachedMethods.Sse2ConvertSingleToDouble, null);
                        il.EmitCall(OpCodes.Call, CachedMethods.Vector128ToScalarDouble, null);
                        break;
                    }
                    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). The struct is on the
                    // stack BY VALUE, so the real part comes from the static by-value helper, not the instance getter.
                    il.EmitCall(OpCodes.Call, CachedMethods.ConvertsComplexToDouble, 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>
        /// <param name="il">The IL stream; the source value is on the evaluation stack.</param>
        /// <param name="from">The source type (Half or Complex, or any type when <paramref name="to"/> is one).</param>
        /// <param name="to">The target type.</param>
        /// <exception cref="NotSupportedException">Neither side is Half nor Complex.</exception>
        /// <remarks>
        /// A float16 crosses with NumPy's bit rules (npy_halfbits_to_floatbits / _to_doublebits,
        /// npy_floatbits_to_halfbits / npy_doublebits_to_halfbits), through the NumPy-exact <see cref="Converts"/>
        /// overloads: a signalling NaN stays signalling. float16 &lt;-&gt; float32 converts DIRECTLY — going through a
        /// double would quiet the NaN in the float32 &lt;-&gt; float64 hop.
        /// </remarks>
        private static void EmitHalfOrComplexConversion(ILGenerator il, NPTypeCode from, NPTypeCode to)
        {
            // Half -> other: convert Half to double first, then to target
            if (from == NPTypeCode.Half)
            {
                // float16 -> float32 DIRECTLY (npy_halfbits_to_floatbits): the double route's narrowing would quiet a
                // signalling NaN.
                if (to == NPTypeCode.Single)
                {
                    il.EmitCall(OpCodes.Call, CachedMethods.ConvertsHalfToSingle, null);
                    return;
                }

                // float16 -> float64 with NumPy's NaN rule (npy_halfbits_to_doublebits; the BCL cast quiets). Every
                // float16 is exact in a double, so it is a lossless intermediate for every other target.
                il.EmitCall(OpCodes.Call, CachedMethods.ConvertsHalfToDouble, 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)
            {
                // The struct is on the stack BY VALUE: the static by-value helper reads its real part (the instance
                // getter Complex.Real would need a managed reference — see CachedMethods.ConvertsComplexToDouble).
                il.EmitCall(OpCodes.Call, CachedMethods.ConvertsComplexToDouble, 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)
            {
                // float32 -> float16 DIRECTLY (npy_floatbits_to_halfbits): widening to float64 first would quiet a
                // signalling NaN (and keep a different slice of its payload).
                if (from == NPTypeCode.Single)
                {
                    il.EmitCall(OpCodes.Call, CachedMethods.ConvertsSingleToHalf, null);
                    return;
                }

                // First convert source to double
                if (from != NPTypeCode.Double)
                    EmitConvertTo(il, from, NPTypeCode.Double);

                // double -> Half, rounded once to nearest-even, with NumPy's NaN rule (npy_doublebits_to_halfbits;
                // the BCL cast quiets a signalling NaN and turns a payload lost in the narrowing into 0x7E00, not 0x7C01)
                il.EmitCall(OpCodes.Call, CachedMethods.ConvertsDoubleToHalf, 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 / np.copysign / np.hypot — 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 || op == BinaryOp.Hypot || op == BinaryOp.Heaviside)
            {
                EmitLogAddNextOperation(il, op, resultType);
                return;
            }

            // np.gcd / np.lcm — one NDGcdLcm helper per integer dtype, signatured (resultType, resultType)
            // -> resultType, so the whole scalar op is a single Call. Intercepted BEFORE the decimal/half/
            // complex routing because gcd/lcm have integer loops ONLY: those dtypes are rejected at the
            // np.* boundary (Default.Gcd/Lcm) and can never reach here as a resultType, so routing them
            // through EmitDecimal/Half/ComplexOperation (which have no gcd/lcm branch) would be wrong.
            if (op == BinaryOp.Gcd || op == BinaryOp.Lcm)
            {
                EmitGcdLcmOperation(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 Fmod - C-style truncated remainder (sign of dividend),
            // unlike Mod's floored (sign of divisor) convention.
            if (op == BinaryOp.Fmod)
            {
                EmitFmodOperation(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 Fmod (np.fmod) via the <see cref="Utilities.NDDivision"/> Fmod* helpers, matching
        /// NumPy's integer <c>@TYPE@_fmod</c> (÷0 -> 0, truncated sign-of-dividend) and the float
        /// <c>npy_fmod</c> (C# <c>%</c> == C fmod). Stack: [dividend, divisor] -> [result].
        /// </summary>
        private static void EmitFmodOperation(ILGenerator il, NPTypeCode resultType)
        {
            var m = GetFmodMethod(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>
        /// Return the <see cref="Utilities.NDDivision"/> fmod helper for <paramref name="resultType"/>,
        /// or null if the dtype routes elsewhere (Half via <see cref="EmitHalfOperation"/>, Decimal via
        /// <see cref="EmitDecimalOperation"/>, Complex unsupported).
        /// </summary>
        private static MethodInfo? GetFmodMethod(NPTypeCode resultType)
        {
            return resultType switch
            {
                NPTypeCode.SByte => CachedMethods.FmodSByte,
                NPTypeCode.Byte => CachedMethods.FmodByte,
                NPTypeCode.Int16 => CachedMethods.FmodInt16,
                NPTypeCode.UInt16 => CachedMethods.FmodUInt16,
                NPTypeCode.Char => CachedMethods.FmodChar,
                NPTypeCode.Int32 => CachedMethods.FmodInt32,
                NPTypeCode.UInt32 => CachedMethods.FmodUInt32,
                NPTypeCode.Int64 => CachedMethods.FmodInt64,
                NPTypeCode.UInt64 => CachedMethods.FmodUInt64,
                NPTypeCode.Single => CachedMethods.FmodSingle,
                NPTypeCode.Double => CachedMethods.FmodDouble,
                _ => null
            };
        }

        /// <summary>
        /// Emit np.gcd / np.lcm via the <see cref="Utilities.NDGcdLcm"/> integer helpers. Stack:
        /// [x1, x2] (both already in the integer <paramref name="resultType"/> — the generic binary loops
        /// convert both operands) -> [result]. A single <c>call</c> to the per-(op,dtype) helper, mirroring
        /// <see cref="EmitFmodOperation"/>. The <paramref name="resultType"/> is always an integer dtype
        /// here (bool/float/complex/decimal are rejected in <c>Default.Gcd</c>/<c>Default.Lcm</c> before any
        /// kernel is generated), so an unresolved dtype is a real bug and throws rather than silently
        /// producing wrong bytes.
        /// </summary>
        /// <param name="il">The kernel's IL stream.</param>
        /// <param name="op"><see cref="BinaryOp.Gcd"/> or <see cref="BinaryOp.Lcm"/>.</param>
        /// <param name="resultType">The integer loop dtype both operands were converted to.</param>
        /// <exception cref="System.NotSupportedException">The (op, dtype) pair has no helper — only reachable
        /// if a non-integer dtype bypassed the np.* validation (a defect).</exception>
        private static void EmitGcdLcmOperation(ILGenerator il, BinaryOp op, NPTypeCode resultType)
        {
            var m = GetGcdLcmMethod(op, resultType)
                ?? throw new NotSupportedException(
                    $"np.{(op == BinaryOp.Gcd ? "gcd" : "lcm")} has no loop for dtype {resultType} " +
                    "(only the integer dtypes are supported).");
            il.EmitCall(OpCodes.Call, m, null);
        }

        /// <summary>
        /// Return the <see cref="Utilities.NDGcdLcm"/> helper for (<paramref name="op"/>,
        /// <paramref name="resultType"/>), or null when the dtype has no gcd/lcm loop (every non-integer
        /// dtype — the np.* boundary raises the no-loop error for those, so null here is a defect signal).
        /// Char rides along as the NumSharp unsigned-16-bit integer extension.
        /// </summary>
        /// <param name="op"><see cref="BinaryOp.Gcd"/> or <see cref="BinaryOp.Lcm"/>.</param>
        /// <param name="resultType">The integer loop dtype.</param>
        /// <returns>The cached 2-arg helper <see cref="MethodInfo"/>, or null if unsupported.</returns>
        private static MethodInfo? GetGcdLcmMethod(BinaryOp op, NPTypeCode resultType)
        {
            if (op == BinaryOp.Gcd)
                return resultType switch
                {
                    NPTypeCode.SByte => CachedMethods.GcdSByte,
                    NPTypeCode.Byte => CachedMethods.GcdByte,
                    NPTypeCode.Int16 => CachedMethods.GcdInt16,
                    NPTypeCode.UInt16 => CachedMethods.GcdUInt16,
                    NPTypeCode.Char => CachedMethods.GcdChar,
                    NPTypeCode.Int32 => CachedMethods.GcdInt32,
                    NPTypeCode.UInt32 => CachedMethods.GcdUInt32,
                    NPTypeCode.Int64 => CachedMethods.GcdInt64,
                    NPTypeCode.UInt64 => CachedMethods.GcdUInt64,
                    _ => null
                };
            return resultType switch
            {
                NPTypeCode.SByte => CachedMethods.LcmSByte,
                NPTypeCode.Byte => CachedMethods.LcmByte,
                NPTypeCode.Int16 => CachedMethods.LcmInt16,
                NPTypeCode.UInt16 => CachedMethods.LcmUInt16,
                NPTypeCode.Char => CachedMethods.LcmChar,
                NPTypeCode.Int32 => CachedMethods.LcmInt32,
                NPTypeCode.UInt32 => CachedMethods.LcmUInt32,
                NPTypeCode.Int64 => CachedMethods.LcmInt64,
                NPTypeCode.UInt64 => CachedMethods.LcmUInt64,
                _ => null
            };
        }

        /// <summary>
        /// Emit np.logaddexp / np.logaddexp2 / np.nextafter / np.copysign / np.hypot / np.heaviside via the
        /// <see cref="Utilities.NDLogAddExpMath"/> (hypot -> <see cref="Utilities.NDHypotMath"/>, heaviside
        /// -> <see cref="Utilities.NDHeavisideMath"/>) 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 scalar helper (<see cref="Utilities.NDLogAddExpMath"/>, or
        /// <see cref="Utilities.NDHypotMath"/> / <see cref="Utilities.NDHeavisideMath"/> for those ops) 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 these 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,
                    };
                case BinaryOp.Hypot:
                    return resultType switch
                    {
                        NPTypeCode.Half => CachedMethods.HypotH,
                        NPTypeCode.Single => CachedMethods.HypotF,
                        NPTypeCode.Decimal => CachedMethods.HypotDec,
                        _ => CachedMethods.HypotD,
                    };
                case BinaryOp.Heaviside:
                    return resultType switch
                    {
                        NPTypeCode.Half => CachedMethods.HeavisideH,
                        NPTypeCode.Single => CachedMethods.HeavisideF,
                        NPTypeCode.Decimal => CachedMethods.HeavisideDec,
                        _ => CachedMethods.HeavisideD,
                    };
                default: // BinaryOp.CopySign
                    return resultType switch
                    {
                        NPTypeCode.Half => CachedMethods.CopySignH,
                        NPTypeCode.Single => CachedMethods.CopySignF,
                        NPTypeCode.Decimal => CachedMethods.CopySignDec,
                        _ => CachedMethods.CopySignD,
                    };
            }
        }

        /// <summary>
        /// Emit ATan2 operation: the angle in radians between the positive x-axis and the point (x, y), with the correct
        /// quadrant. float32 is the CRT's <c>atan2f</c> (<see cref="MathF.Atan2(float, float)"/>, NumPy's float32 loop);
        /// every other loop dtype reaching here widens to double and calls <see cref="Math.Atan2(double, double)"/>
        /// (float16 takes <see cref="ATan2HalfNumPy"/> through <c>EmitHalfOperation</c>, decimal its own helper).
        /// Stack: [y, x] -> [result] (range [-pi, pi]).
        /// </summary>
        /// <param name="il">The IL stream; both operands are already converted to <paramref name="resultType"/>.</param>
        /// <param name="resultType">The loop dtype.</param>
        private static void EmitATan2Operation(ILGenerator il, NPTypeCode resultType)
        {
            // float32: NumPy's FLOAT_arctan2 is the CRT's atan2f, and MathF.Atan2 calls that same UCRT function (parity
            // QA-24). Widening to double, Math.Atan2 and rounding back is the correctly rounded angle, one float32 ULP from
            // atan2f wherever atan2f is not (116 of 2M probed pairs). Stack: [float y, float x] -> [float].
            if (resultType == NPTypeCode.Single)
            {
                il.EmitCall(OpCodes.Call, CachedMethods.MathFAtan2, null);
                return;
            }

            // 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;
            }

            // Fmod for decimal: C-style truncated remainder (decimal '%' truncates toward zero).
            if (op == BinaryOp.Fmod)
            {
                // Stack: [dividend, divisor] -> NDDivision.FmodDecimal(a, b) == a % b
                il.EmitCall(OpCodes.Call, CachedMethods.FmodDecimal, 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");

            // FloorDivide / Mod / Fmod: NumPy's HALF loops compute in FLOAT32 (astype 'e'->'f' —
            // HALF_floor_divide/remainder/divmod call npy_floor_dividef/npy_remainderf/npy_fmodf),
            // so route through the NDDivision *Single helpers rather than the generic double path.
            // The old double `a - floor(a/b)*b` / `floor(a/b) with inf->NaN` diverged from NumPy on
            // ÷0 (NaN vs the floored quotient / IEEE ±inf) — the W1-A known bug. A Half is exact in
            // both float and double, so HalfToDouble+Conv_R4 is an exact (float) widen and
            // Conv_R8+DoubleToHalf an exact narrow (== (Half)floatResult), keeping the arithmetic in
            // float32 as NumPy does. Stack in: [half1, half2].
            if (op == BinaryOp.FloorDivide || op == BinaryOp.Mod || op == BinaryOp.Fmod)
            {
                var locRhsHalf = il.DeclareLocal(typeof(Half));
                il.Emit(OpCodes.Stloc, locRhsHalf);                 // [half1]
                il.EmitCall(OpCodes.Call, CachedMethods.HalfToDouble, null);
                il.Emit(OpCodes.Conv_R4);                            // [floatA]
                il.Emit(OpCodes.Ldloc, locRhsHalf);
                il.EmitCall(OpCodes.Call, CachedMethods.HalfToDouble, null);
                il.Emit(OpCodes.Conv_R4);                            // [floatA, floatB]
                var single = op switch
                {
                    BinaryOp.FloorDivide => CachedMethods.FloorDivSingle,
                    BinaryOp.Mod => CachedMethods.RemSingle,
                    _ => CachedMethods.FmodSingle,
                };
                il.EmitCall(OpCodes.Call, single, null);             // [floatResult]
                il.Emit(OpCodes.Conv_R8);
                il.EmitCall(OpCodes.Call, CachedMethods.DoubleToHalf, null); // [halfResult]
                return;
            }

            // Add / Subtract / Multiply / Divide: NumPy's HALF loops compute in FLOAT32
            // (astype 'e'->'f'), NOT double. The double bridge below double-rounds where NumPy
            // single-rounds (exponent-gap sums — a real 1-ULP divergence) and BCL-quiets sNaN, so
            // route them through HalfArithScalarStruct (== HalfArithBits, the float32 + RTNE-narrow +
            // operand-order NaN pin the engine's own Half kernel and np.evaluate's vector path use).
            // Stack in: [Half a, Half b] -> push op -> [Half a, Half b, op] -> [Half result].
            if (op == BinaryOp.Add || op == BinaryOp.Subtract || op == BinaryOp.Multiply || op == BinaryOp.Divide)
            {
                il.Emit(OpCodes.Ldc_I4, (int)op);
                il.EmitCall(OpCodes.Call, GetHelper(nameof(HalfArithScalarStruct)), null);
                return;
            }

            // Power: NumPy's HALF_power is powf on the halves widened bit for bit (a signalling NaN stays signalling)
            // and narrowed with npy_float_to_half — not the double bridge below, which quiets the NaN before pow
            // (sNaN ** 0 became 1), loses a negative NaN's sign on integer exponents and rounds from double.
            // Stack in: [Half a, Half b] -> [Half result].
            if (op == BinaryOp.Power)
            {
                il.EmitCall(OpCodes.Call, GetHelper(nameof(PowerHalfNumPy)), null);
                return;
            }

            // ATan2: NumPy's HALF_arctan2 is atan2f on the halves widened bit for bit, narrowed with npy_float_to_half
            // (ATan2HalfNumPy, parity QA-24) — not the double bridge below, whose single rounding from Math.Atan2 lands on
            // the other float16 neighbour where atan2f's float32 result is a float16 midpoint. Stack in: [Half y, Half x].
            if (op == BinaryOp.ATan2)
            {
                il.EmitCall(OpCodes.Call, GetHelper(nameof(ATan2HalfNumPy)), null);
                return;
            }

            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;
                // Mod / FloorDivide / Fmod are handled above via the float32 NDDivision path
                // (NumPy's HALF loops compute in float32), so they never reach this double switch.
                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;
            }

            // Power goes through a NumPy-compatible helper rather than the BCL's Complex.Pow: NumPy's
            // npy_cpow special-cases integer exponents with EXACT repeated multiplication, whereas the
            // BCL always takes the polar exp(b*log a) form — up to ~1e14 ULP off on z**3/z**5/z**-1
            // (effectively wrong), and even z**2 diverges from np.square. See ComplexPowNumPy.
            if (op == BinaryOp.Power)
            {
                il.EmitCall(OpCodes.Call, GetHelper(nameof(ComplexPowNumPy)), null);
                return;
            }

            // Multiply goes through NDComplexMath.Multiply (FMA-contracted, == NumPy's simd_cmul)
            // rather than the BCL's naive op_Multiply — see the ComplexMultiply cache comment.
            if (op == BinaryOp.Multiply)
            {
                il.EmitCall(OpCodes.Call, CachedMethods.ComplexMultiply, null);
                return;
            }

            var method = op switch
            {
                BinaryOp.Add => complexType.GetMethod("op_Addition", new[] { complexType, complexType }),
                BinaryOp.Subtract => complexType.GetMethod("op_Subtraction", 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 — a byte-for-byte transcription of the CDOUBLE_divide
        /// scalar loop in NumPy 2.4.2 (<c>numpy/_core/src/umath/loops.c.src</c>), Smith's algorithm.
        /// </summary>
        /// <remarks>
        /// This must NOT defer to the BCL's <see cref="System.Numerics.Complex"/> <c>op_Division</c>:
        /// both use Smith's algorithm, but they differ in the last bit for finite inputs on ~33% of
        /// random operands. Two differences, each observable:
        /// <list type="bullet">
        /// <item>NumPy computes the reciprocal <c>scl = 1/denom</c> ONCE and MULTIPLIES; the BCL DIVIDES
        /// each component by <c>denom</c>. So <c>z / 3.0</c> is <c>z*(1.0/3.0)</c> in NumPy but <c>z/3.0</c>
        /// in the BCL — a 1-ULP gap wherever <c>a*(1/c) != a/c</c> (the reported footgun; trapezoid dodged
        /// it only because <c>/2.0</c> is a power of two, hence exact both ways).</item>
        /// <item>The branch pivot is <c>|c| &gt;= |d|</c> (NumPy) vs the BCL's <c>|d| &lt; |c|</c>, which pick
        /// different formulas exactly at <c>|c| == |d|</c>.</item>
        /// </list>
        /// The multiplications must stay un-fused (plain <c>*</c>/<c>+</c>, never <see cref="System.Math.FusedMultiplyAdd"/>):
        /// NumPy's win-amd64 wheel builds this loop under MSVC <c>/fp:precise</c>, which does not contract,
        /// and RyuJIT never introduces an FMA on its own — so the separate ops match bit-for-bit.
        /// The divide-by-zero sub-branch divides by <c>|c|</c>/<c>|d|</c> (both <c>+0.0</c> via <see cref="System.Math.Abs(double)"/>),
        /// so <c>z/(0+0j)</c> yields NumPy's component-wise IEEE result (±inf / NaN) rather than the BCL's (NaN, NaN).
        /// </remarks>
        // 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)
        {
            double in1r = a.Real, in1i = a.Imaginary;
            double in2r = b.Real, in2i = b.Imaginary;
            double in2r_abs = System.Math.Abs(in2r);
            double in2i_abs = System.Math.Abs(in2i);
            if (in2r_abs >= in2i_abs)
            {
                // Divisor real part dominates (covers the real-scalar case, in2i == 0).
                if (in2r_abs == 0.0 && in2i_abs == 0.0)
                    // Divide by zero → component-wise IEEE division by +0.0 (NumPy's complex inf/nan).
                    return new System.Numerics.Complex(in1r / in2r_abs, in1i / in2i_abs);

                double rat = in2i / in2r;
                double scl = 1.0 / (in2r + in2i * rat);
                return new System.Numerics.Complex((in1r + in1i * rat) * scl, (in1i - in1r * rat) * scl);
            }
            else
            {
                // Divisor imaginary part dominates.
                double rat = in2r / in2i;
                double scl = 1.0 / (in2i + in2r * rat);
                return new System.Numerics.Complex((in1r * rat + in1i) * scl, (in1i * rat - in1r) * scl);
            }
        }

        /// <summary>
        /// Naive complex multiply — a byte-for-byte match of NumPy's internal <c>cmul</c>
        /// (<c>npy_math_complex.c.src</c>: <c>(ar*br - ai*bi, ar*bi + ai*br)</c>, un-fused). Used by
        /// <see cref="ComplexPowNumPy"/> so that integer-power repeated multiplication rounds exactly
        /// as NumPy's does. Identical to the BCL <c>Complex.op_Multiply</c> for this scalar formula,
        /// but written out so the power port does not depend on the BCL operator's internal shape.
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        private static System.Numerics.Complex ComplexMulNumPy(System.Numerics.Complex a, System.Numerics.Complex b)
        {
            double ar = a.Real, ai = a.Imaginary, br = b.Real, bi = b.Imaginary;
            return new System.Numerics.Complex(ar * br - ai * bi, ar * bi + ai * br);
        }

        /// <summary>
        /// NumPy-compatible complex power — a transcription of NumPy 2.4.2's <c>npy_cpow</c>
        /// (<c>numpy/_core/src/npymath/npy_math_complex.c.src</c>).
        /// </summary>
        /// <remarks>
        /// This must NOT defer to the BCL's <see cref="System.Numerics.Complex"/> <c>Pow</c> for the
        /// common cases: <c>Complex.Pow</c> ALWAYS evaluates the polar form <c>exp(b·log a)</c>, which
        /// for an INTEGER exponent is catastrophically imprecise — measured up to ~1e14 raw ULP on
        /// <c>z**3</c>/<c>z**5</c>/<c>z**-1</c> (effectively the wrong value), and even <c>z**2</c>
        /// diverges from <see cref="System.Numerics.Complex"/>-free <c>np.square</c>. NumPy instead:
        /// <list type="bullet">
        /// <item><c>b == 0</c> → <c>1+0j</c> (incl. <c>0**0</c>).</item>
        /// <item><c>a == 0</c>, <c>b != 0</c> → <c>0+0j</c> when <c>Re(b) &gt; 0</c>, else complex NaN
        /// (<c>0</c> to a negative/imaginary power; NumPy raises the invalid flag via <c>inf-inf</c>).</item>
        /// <item>Integer <c>Re(b)</c> in (-100, 100) with <c>Im(b) == 0</c> → EXACT repeated
        /// multiplication (<c>z*z</c>, binary exponentiation via <see cref="ComplexMulNumPy"/>), with
        /// <c>1/result</c> through the byte-exact <see cref="ComplexDivideNumPy"/> for a negative exponent.</item>
        /// <item>Everything else (non-integer real / complex exponent) → the platform path
        /// <see cref="System.Numerics.Complex.Pow(System.Numerics.Complex,System.Numerics.Complex)"/>;
        /// NumPy uses the host <c>cpow</c> there, so a bounded ~ULP / libm divergence remains for that
        /// branch only (documented in the oracle).</item>
        /// </list>
        /// The <c>a == 0</c> branch is checked BEFORE the integer branch, exactly as NumPy orders them —
        /// otherwise <c>0**-1</c> would take the repeated-multiply path (<c>1/0 → inf/nan</c>) instead of
        /// yielding NumPy's complex NaN. The complex-NaN result uses the POSITIVE quiet-NaN payload
        /// NumPy's <c>NPY_NAN</c> carries (<c>0x7ff8…</c>, not .NET's negative <c>double.NaN</c>) so the
        /// bytes match should complex-power NaN ever become contractual.
        /// </remarks>
        // Per-element op called from the binary IL kernel's inner loop.
        [MethodImpl(MethodImplOptions.AggressiveInlining | MethodImplOptions.AggressiveOptimization)]
        internal static System.Numerics.Complex ComplexPowNumPy(System.Numerics.Complex a, System.Numerics.Complex b)
        {
            double ar = a.Real, ai = a.Imaginary;
            double br = b.Real, bi = b.Imaginary;

            // a**0 == 1 for any a (including 0**0), by the definition of the logarithm.
            if (br == 0.0 && bi == 0.0)
                return new System.Numerics.Complex(1.0, 0.0);

            // 0**b (b != 0): 0 for Re(b) > 0, otherwise complex NaN (negative / imaginary power of 0).
            if (ar == 0.0 && ai == 0.0)
            {
                if (br > 0.0)
                    return new System.Numerics.Complex(0.0, 0.0);
                // Raise the IEEE invalid flag the way NumPy does (inf - inf), then return complex NaN.
                double invalid = double.PositiveInfinity - double.PositiveInfinity;
                _ = invalid;
                return new System.Numerics.Complex(PositiveNaN, PositiveNaN);
            }

            // Integer real exponent strictly within (-100, 100): exact repeated multiplication.
            long n;
            if (bi == 0.0 && br > -100.0 && br < 100.0 && (n = (long)br) == br)
            {
                if (n == 1) return new System.Numerics.Complex(ar, ai);           // unroll (inf-friendly)
                if (n == 2) return ComplexMulNumPy(a, a);                          // == np.square
                if (n == 3) return ComplexMulNumPy(a, ComplexMulNumPy(a, a));

                // General |n| < 100 via binary exponentiation over ComplexMulNumPy.
                long nn = n < 0 ? -n : n;
                var aa = new System.Numerics.Complex(1.0, 0.0);
                var p = new System.Numerics.Complex(ar, ai);
                long mask = 1;
                while (true)
                {
                    if ((nn & mask) != 0)
                        aa = ComplexMulNumPy(aa, p);
                    mask <<= 1;
                    if (nn < mask || mask <= 0)   // mask <= 0 guards the (theoretical) shift overflow
                        break;
                    p = ComplexMulNumPy(p, p);
                }
                // Negative exponent → reciprocal through the byte-exact complex divide (NumPy's cdiv(1, r)).
                return br < 0.0 ? ComplexDivideNumPy(new System.Numerics.Complex(1.0, 0.0), aa) : aa;
            }

            // Non-integer real or complex exponent: NumPy uses the host cpow; the BCL's polar
            // exp(b*log a) is the closest managed equivalent (documented bounded libm/~ULP divergence).
            return System.Numerics.Complex.Pow(a, b);
        }

        /// <summary>NumPy's positive quiet NaN (<c>NPY_NAN</c>, <c>0x7ff8000000000000</c>) — .NET's
        /// <see cref="double.NaN"/> is the NEGATIVE <c>0xfff8…</c>, so this constant is used where a
        /// complex result must carry NumPy's NaN bytes.</summary>
        private static readonly double PositiveNaN =
            System.BitConverter.Int64BitsToDouble(unchecked((long)0x7ff8000000000000UL));

        /// <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>
        // Boolean rides byte lanes; Half rides its raw 16-bit (ushort) lanes — Vector{N}&lt;Half&gt;
        // throws (no BCL f16 vector arithmetic), so every SIMD Half kernel in the codebase works on
        // the bit patterns and widens to float32 on demand. No caller passes Half here except the
        // np.evaluate fused Half-arithmetic path (CanUseSimd(Half) is false, so the engine's own
        // per-op SIMD gates never reach Half), which is exactly the path that needs the ushort lane.
        internal static Type GetSimdLaneType(NPTypeCode type)
            => type == NPTypeCode.Boolean ? typeof(byte)
             : type == NPTypeCode.Half ? typeof(ushort)
             : 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>
        /// Width-explicit <see cref="EmitVectorLoad"/>: emit a vector load of <paramref name="bits"/>
        /// (128/256/512) at CLR lane type <paramref name="clrType"/>, regardless of the module-wide
        /// <see cref="VectorBits"/>. The fused MIXED-WIDTH path (ndexpr-evaluate.md P5.2) needs this
        /// because a half-lane operand (e.g. int32 under a float64 lane) loads a PARTIAL vector —
        /// <c>Vector128&lt;int&gt;</c> on a 256-bit host — of exactly <c>elemCount·sizeof(T)</c> bytes,
        /// so the load never over-reads past the operand's last element.
        /// </summary>
        /// <param name="il">The IL generator; stack [T*] → [Vector{bits}&lt;T&gt;].</param>
        /// <param name="bits">The container width in bits (128/256/512).</param>
        /// <param name="clrType">The CLR element type of the vector.</param>
        internal static void EmitVectorLoadAt(ILGenerator il, int bits, Type clrType)
        {
            // Same x86-preference rationale as EmitVectorLoad: the JIT generates ~1.8x faster code
            // for Avx/Sse loads than for the cross-platform Vector{N}.Load path.
            var x86 = VectorMethodCache.LoadX86(bits, clrType);
            if (x86 != null)
            {
                il.EmitCall(OpCodes.Call, x86, null);
                return;
            }
            il.EmitCall(OpCodes.Call, VectorMethodCache.Load(bits, clrType), null);
        }

        /// <summary>
        /// Width-explicit <see cref="EmitVectorCreate"/>: broadcast the scalar on the stack into a
        /// vector of <paramref name="bits"/> at <paramref name="type"/>'s SIMD lane type. Used by the
        /// fused mixed-width path, where a half-lane constant / broadcast operand / hoisted parameter
        /// builds a <c>Vector128</c> while the kernel's full lane rides <c>Vector256</c>.
        /// </summary>
        /// <param name="il">The IL generator; stack [scalar] → [Vector{bits}&lt;lane(type)&gt;].</param>
        /// <param name="type">The element dtype whose SIMD lane type the broadcast targets.</param>
        /// <param name="bits">The container width in bits (128/256/512).</param>
        internal static void EmitVectorCreateAt(ILGenerator il, NPTypeCode type, int bits)
            => il.EmitCall(OpCodes.Call, VectorMethodCache.CreateBroadcast(bits, GetSimdLaneType(type)), null);

        /// <summary>
        /// Width-explicit <see cref="EmitVectorOperation"/>: the binary vector op at an EXPLICIT
        /// container width. At <paramref name="bits"/> == <see cref="VectorBits"/> it forwards to
        /// <see cref="EmitVectorOperation"/> byte-for-byte (so uniform-lane callers see zero change);
        /// a NARROWER container serves the fused mixed-width path's half-lane interior nodes —
        /// e.g. the int32 multiply of <c>i4*2+f8</c> runs at <c>Vector128&lt;int&gt;</c> so it WRAPS at
        /// int32 exactly like NumPy's unfused sequence, before the edge widens to the f8 lane.
        /// </summary>
        /// <param name="il">The IL generator; stack [a, b] (both Vector{bits}&lt;T&gt;) → [result].</param>
        /// <param name="op">The binary op — the mixed plan admits only the arithmetic four and the
        /// bitwise three at a narrow container (min/max/other ops decline the whole tree to scalar).</param>
        /// <param name="type">The node's element dtype (never Boolean/Half at a narrow container —
        /// the mixed plan excludes bool-typed nodes and Half trees take the dedicated P5.1 path).</param>
        /// <param name="bits">The container width in bits.</param>
        /// <exception cref="NotSupportedException">The op has no vector body at a narrow container.</exception>
        internal static void EmitVectorOperationAt(ILGenerator il, BinaryOp op, NPTypeCode type, int bits)
        {
            // The full-width call must remain byte-identical to the pre-P5.2 emission (it carries the
            // Boolean byte-lane logical remap and the NaN-aware min/max wrappers the narrow path
            // never needs), so forward rather than duplicate.
            if (bits == VectorBits)
            {
                EmitVectorOperation(il, op, type);
                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(bits, methodName, clrType);
                if (x86Bit != null)
                {
                    il.EmitCall(OpCodes.Call, x86Bit, null);
                    return;
                }
                il.EmitCall(OpCodes.Call, VectorMethodCache.Generic(bits, methodName, clrType, paramCount: 2), null);
                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 a narrow SIMD container"),
            };
            var x86Arith = VectorMethodCache.BinaryX86(bits, arithName, clrType);
            if (x86Arith != null)
            {
                il.EmitCall(OpCodes.Call, x86Arith, null);
                return;
            }
            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 a narrow SIMD container"),
            };
            il.EmitCall(OpCodes.Call, VectorMethodCache.Operator(bits, clrType, operatorName), 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

    }
}