[Major Rewrite] NumPy nditer port, NpyExpr DSL with 3-tier custom-op API, C/F/A/K memory layout support, stride-native matmul

.agents/skills/np-function/SKILL.md

@@ -0,0 +1 @@
+../../../.claude/skills/np-function/SKILL.md

.agents/skills/np-tests/SKILL.md

@@ -0,0 +1 @@
+../../../.claude/skills/np-tests/SKILL.md

.claude/CLAUDE.md

@@ -27,8 +27,8 @@ Every np.* function and DefaultEngine operation MUST satisfy these criteria:
 - **Sliced views**: Correctly handles Shape.offset for base address calculation
 
 ### Dtype Support
-All 12 NumSharp types must be handled (or explicitly documented as unsupported):
-Boolean, Byte, Int16, UInt16, Int32, UInt32, Int64, UInt64, Char, Single, Double, Decimal
+All 15 NumSharp types must be handled (or explicitly documented as unsupported):
+Boolean, Byte, SByte, Int16, UInt16, Int32, UInt32, Int64, UInt64, Char, Half, Single, Double, Decimal, Complex
 
 ### NumPy API Parity
 - Function signature matches NumPy (parameter names, order, defaults)
@@ -43,18 +43,26 @@ Boolean, Byte, Int16, UInt16, Int32, UInt32, Int64, UInt64, Char, Single, Double
 
 **Full audit tracking:** See `docs/KERNEL_API_AUDIT.md`
 
-## Supported Types (12)
+## Supported Types (15)
 
 | NPTypeCode | C# Type | NPTypeCode | C# Type |
 |------------|---------|------------|---------|
-| Boolean | bool | Int64 | long |
-| Byte | byte | UInt64 | ulong |
-| Int16 | short | Char | char |
-| UInt16 | ushort | Single | float |
-| Int32 | int | Double | double |
-| UInt32 | uint | Decimal | decimal |
-
-All operations must handle all 12 types via type switch pattern.
+| Boolean | bool | UInt64 | ulong |
+| Byte | byte | Char | char |
+| SByte | sbyte | Half | System.Half |
+| Int16 | short | Single | float |
+| UInt16 | ushort | Double | double |
+| Int32 | int | Decimal | decimal |
+| UInt32 | uint | Complex | System.Numerics.Complex |
+| Int64 | long |  |  |
+
+All operations must handle all 15 types via type switch pattern.
+
+**Perf notes:**
+- SByte / Byte / Int*/UInt* / Single / Double — full SIMD via `MixedTypeKernel.SimdFull` (V128/V256/V512 detected at startup).
+- Half — scalar path (no `Vector<Half>` arithmetic in .NET BCL). Routes through `Half→double→Math.Pow→Half` for `np.power`; ~2× slower than NumPy.
+- Complex — scalar path via `System.Numerics.Complex` operators / `Complex.Pow`. ~2× slower than NumPy.
+- Decimal — scalar path via `DecimalMath.Pow`. Highest precision, slowest.
 
 ## Architecture
 
@@ -74,7 +82,7 @@ np                Static API class (like `import numpy as np`)
 | Decision | Rationale |
 |----------|-----------|
 | Unmanaged memory | Benchmarked fastest; optimized for performance |
-| C-order only | Only row-major (C-order) memory layout. Uses `ArrayFlags.C_CONTIGUOUS` flag. No F-order/column-major support. The `order` parameter on `ravel`, `flatten`, `copy`, `reshape` is accepted but ignored. |
+| Order-aware layout | Row-major (C-order) remains the default. `Shape` also tracks F-contiguity, and APIs with an `order` parameter resolve NumPy `C`/`F`/`A`/`K` modes through `OrderResolver`. |
 | Regen templating | Type-specific code generation (legacy, mostly replaced by ILKernel) |
 | TensorEngine abstract | Future GPU/SIMD backends possible |
 | View semantics | Slicing returns views (shared memory), not copies |
@@ -146,14 +154,15 @@ public readonly partial struct Shape
 | Flag | Value | Meaning |
 |------|-------|---------|
 | `C_CONTIGUOUS` | 0x0001 | Data is row-major contiguous |
-| `F_CONTIGUOUS` | 0x0002 | Reserved (always false for NumSharp) |
+| `F_CONTIGUOUS` | 0x0002 | Data is column-major contiguous |
 | `OWNDATA` | 0x0004 | Array owns its data buffer |
 | `ALIGNED` | 0x0100 | Always true for managed allocations |
 | `WRITEABLE` | 0x0400 | False for broadcast views |
 | `BROADCASTED` | 0x1000 | Has stride=0 with dim > 1 |
 
 **Key Shape properties:**
 - `IsContiguous` — O(1) check via `C_CONTIGUOUS` flag
+- `IsFContiguous` — O(1) check via `F_CONTIGUOUS` flag
 - `IsBroadcasted` — O(1) check via `BROADCASTED` flag
 - `IsWriteable` — False for broadcast views (prevents corruption)
 - `IsSliced` — True if offset != 0, different size, or non-contiguous
@@ -182,15 +191,14 @@ nd["..., -1"]     // Ellipsis fills dimensions
 
 ---
 
-## Missing Functions (20)
+## Missing Functions (18)
 
 These NumPy functions are **not implemented**:
 
 | Category | Functions |
 |----------|-----------|
 | Sorting | `np.sort` |
-| Selection | `np.where` |
-| Manipulation | `np.flip`, `np.fliplr`, `np.flipud`, `np.rot90`, `np.tile`, `np.pad` |
+| Manipulation | `np.flip`, `np.fliplr`, `np.flipud`, `np.rot90`, `np.pad` |
 | Splitting | `np.split`, `np.array_split`, `np.hsplit`, `np.vsplit`, `np.dsplit` |
 | Diagonal | `np.diag`, `np.diagonal`, `np.trace` |
 | Cumulative | `np.diff`, `np.gradient`, `np.ediff1d` |
@@ -206,7 +214,7 @@ Tested against NumPy 2.x.
 `arange`, `array`, `asanyarray`, `asarray`, `copy`, `empty`, `empty_like`, `eye`, `frombuffer`, `full`, `full_like`, `identity`, `linspace`, `meshgrid`, `mgrid`, `ones`, `ones_like`, `zeros`, `zeros_like`
 
 ### Shape Manipulation
-`atleast_1d`, `atleast_2d`, `atleast_3d`, `concatenate`, `dstack`, `expand_dims`, `flatten`, `hstack`, `moveaxis`, `ravel`, `repeat`, `reshape`, `roll`, `rollaxis`, `squeeze`, `stack`, `swapaxes`, `transpose`, `unique`, `vstack`
+`atleast_1d`, `atleast_2d`, `atleast_3d`, `concatenate`, `dstack`, `expand_dims`, `flatten`, `hstack`, `moveaxis`, `ravel`, `repeat`, `reshape`, `roll`, `rollaxis`, `squeeze`, `stack`, `swapaxes`, `tile`, `transpose`, `unique`, `vstack`
 
 ### Broadcasting
 `are_broadcastable`, `broadcast`, `broadcast_arrays`, `broadcast_to`
@@ -226,6 +234,9 @@ Tested against NumPy 2.x.
 ### Comparison & Logic
 `all`, `allclose`, `any`, `array_equal`, `find_common_type`, `isclose`, `isfinite`, `isinf`, `isnan`, `isscalar`, `maximum`, `minimum`
 
+### Selection
+`where`
+
 ### Sorting & Searching
 `argmax`, `argmin`, `argsort`, `nonzero`, `searchsorted`
 
@@ -264,7 +275,7 @@ Tested against NumPy 2.x.
 | TensorEngine | `Backends/TensorEngine.cs` |
 | DefaultEngine | `Backends/Default/DefaultEngine.*.cs` |
 | np API | `APIs/np.cs` |
-| Iterators | `Backends/Iterators/NDIterator.cs`, `MultiIterator.cs` |
+| Iterators | `Backends/Iterators/NDIterator.cs`, `NpyIter.cs`, `NpyExpr.cs` |
 | Type info | `Utilities/InfoOf.cs` |
 | Generic NDArray | `Generics/NDArray\`1.cs` |
 
@@ -574,10 +585,10 @@ A: The `Slice` class parses Python notation (e.g., "1:5:2") into `Start`, `Stop`
 A: `Slice.All` (`:` - all elements), `Slice.Ellipsis` (`...` - fill dimensions), `Slice.NewAxis` (insert dimension), `Slice.Index(n)` (single element, reduces dimensionality).
 
 **Q: What is NDIterator used for?**
-A: Traversing arrays with different memory layouts. Handles contiguous (fast pointer increment) and sliced (uses GetOffset) arrays. Has `MoveNext()`, `HasNext()`, `Reset()`. AutoReset mode for broadcasting smaller arrays.
+A: Legacy typed traversal surface over `NpyIter`. It keeps the existing `MoveNext()`, `HasNext()`, and `Reset()` API while delegating stride, broadcast, and view traversal to the NpyIter state machinery.
 
-**Q: What is MultiIterator?**
-A: Handles paired iteration for broadcasting. `MultiIterator.Assign(lhs, rhs)` copies with broadcasting. `GetIterators(lhs, rhs, broadcast)` creates synchronized iterators.
+**Q: What is NpyIter?**
+A: The NumPy-aligned multi-operand iterator. It handles C/F/A/K order, broadcasting, external loops, buffering, casting, masks, reductions, and synchronized traversal for copy and elementwise kernels. `MultiIterator` was removed in favor of `NpyIter.Copy` and multi-operand iterator execution.
 
 **Q: How does broadcasting work?**
 A: Shapes align from the right. Dimensions must be equal OR one must be 1. Dimension of 1 "stretches" to match. Implemented via `DefaultEngine.Broadcast()` which resolves compatible shapes.

AGENTS.md

@@ -0,0 +1 @@
+.claude/CLAUDE.md

ARCHITECTURE.md

@@ -480,6 +480,8 @@ MoveNext = () => *((T*)Address + shape.GetOffset(index++));
 
 ## Code Generation
 
+For the practical implementation rules used by `DefaultEngine` and `ILKernelGenerator`, see `docs/DEFAULTENGINE_ILKERNEL_PLAYBOOK.md`. That guide captures the recurring engine patterns, optimization conventions, and test expectations that are only implicit in the source code.
+
 ### Regen Templating
 
 NumSharp uses Regen (a custom templating engine) to generate type-specific code. This results in approximately **200,000 lines of generated code**.

benchmark/NumSharp.Benchmark.Exploration/Program.cs

@@ -393,7 +393,7 @@ private class Options
             Dtypes = Dtypes,
             Sizes = Sizes,
             OutputPath = OutputPath,
-            RemainingArgs = RemainingArgs
+            RemainingArgs = (string[])RemainingArgs.Clone()
         };
     }
 }