after merge from remote

src/Nitload.Mathematics/.discard/CoordinationNaming.md

@@ -0,0 +1,88 @@
+# Matrix Coordinate & Naming Conventions (StgSharp)
+
+This document defines the coordinate terminology and naming rules used across `Matrix<T>`, `MatrixKernel<T>`, and packed/triangular storage layouts.  
+Goal: eliminate ambiguity around `column`/`row` names (e.g., “index **of** a column” vs “index **along** a column”).
+
+---
+
+## 1) The key rule: encode prepositions into identifiers (`Of` vs `In`)
+
+### `...OfX` — “which X” (object identity / ordinal)
+Use **`Of`** when the value identifies **which object** (which column, which row, which kernel, etc.).
+
+Examples:
+- `columnOfElement`: the element is located in **which column**
+- `rowOfKernel`: the kernel block is located in **which row**
+
+### `...InX` / `...WithinX` — “inside X” (local coordinate)
+Use **`In`** (or `Within`) when the value is a **local coordinate inside** another object.
+
+Examples:
+- `columnInKernel`: local column coordinate **inside** a kernel block
+- `rowInKernel`: local row coordinate **inside** a kernel block
+
+> Quick mental model: **`Of` = identity**, **`In` = inside**.
+
+---
+
+## 2) Coordinate levels: `Element`, `KernelGrid`, `InKernel`
+
+StgSharp commonly uses three coordinate levels:
+
+1. **Element-level**: final matrix element coordinates (`Matrix<T>[col,row]`)
+2. **Kernel-grid level**: coordinates of a kernel block in the kernel grid
+3. **In-kernel level**: local element coordinates within a kernel (e.g., 4×4)
+
+Recommended names:
+
+| Level | Column name | Row name | Meaning |
+|---|---|---|---|
+| Element | `columnOfElement` | `rowOfElement` | element located in which column/row (0-based) |
+| Kernel grid | `columnOfKernel` | `rowOfKernel` | kernel block located in which kernel-column/kernel-row (0-based) |
+| In-kernel | `columnInKernel` | `rowInKernel` | local coordinate inside the kernel (0..KernelSize-1) |
+
+---
+
+## 3) Index base: 0-based (mandatory)
+
+Unless explicitly documented otherwise, all coordinates are **0-based**:
+
+- `columnOfElement ∈ [0, ColumnLength-1]`
+- `rowOfElement ∈ [0, RowLength-1]`
+- `columnOfKernel ∈ [0, KernelColumnLength-1]`
+- `rowOfKernel ∈ [0, KernelRowLength-1]`
+- `columnInKernel ∈ [0, KernelSize-1]`
+- `rowInKernel ∈ [0, KernelSize-1]`
+
+---
+
+## 4) Mapping element coordinates to kernel coordinates (KernelSize = 4 example)
+
+For a fixed kernel size `KernelSize = 4`:
+
+```csharp
+int columnOfKernel = columnOfElement / 4;
+int rowOfKernel    = rowOfElement    / 4;
+
+int columnInKernel = columnOfElement % 4;
+int rowInKernel    = rowOfElement    % 4;
+```
+
+If `columnOfElement` and `rowOfElement` are guaranteed non-negative, bit operations can be used:
+
+```csharp
+int columnOfKernel = columnOfElement >> 2;
+int rowOfKernel    = rowOfElement    >> 2;
+
+int columnInKernel = columnOfElement & 3;
+int rowInKernel    = rowOfElement    & 3;
+```
+
+---
+
+## 5) Packed storage: distinguish `Index` from `Offset`
+
+Packed layouts do **not** store “empty blocks”; stored blocks are placed contiguously in a packed buffer.
+
+### Use the correct terms
+- **`packedKernelIndex`**: “the N-th kernel

src/Nitload.Mathematics/.discard/MatrixCompute.Factorize.cs.old

@@ -0,0 +1,160 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixCompute.Factorize"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using Nitload.Common.Internal;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static partial class MatrixCompute
+    {
+
+        public static unsafe void FactorizePLU(
+                                  Matrix<float> matrix,
+                                  Span<int> pivot,
+                                  Matrix<float> lower,
+                                  Matrix<float> upper
+        )
+        {
+            if (matrix.ColumnLength != matrix.RowLength) {
+                throw new InvalidOperationException("Matrix must be square.");
+            }
+            if (matrix.Layout != MatrixLayout.DenseRectangle) {
+                throw new NotSupportedException("Only dense square matrices are supported.");
+            }
+            int n = matrix.RowLength;
+            if (pivot.Length < n) {
+                throw new ArgumentException("Pivot span length is insufficient.", nameof(pivot));
+            }
+            if (lower.RowLength != n ||
+                lower.ColumnLength != n ||
+                lower.Layout != MatrixLayout.UpperTriangle) {
+                throw new ArgumentException("Lower matrix must be lower-triangular with same dimensions as input.", nameof(lower));
+            }
+            if (upper.RowLength != n ||
+                upper.ColumnLength != n ||
+                upper.Layout != MatrixLayout.UpperTriangle) {
+                throw new ArgumentException("Upper matrix must be upper-triangular with same dimensions as input.", nameof(upper));
+            }
+
+            int size = matrix.KernelColumnLength * matrix.KernelRowLength;
+            if (size < 256)
+            {
+                for (int i = 0; i < n; i++)
+                {
+                    pivot[i] = i; // identity permutation
+                }
+
+                // Aggressive: always call native pivot to fetch four-column candidates
+                ulong* pivots = stackalloc ulong[4];
+                for (int kBlock = 0; kBlock < n; kBlock += 4)
+                {
+                    NumericalModule.GlobalContext
+                                   .mat_sk[(int)MatrixElementType.F32].pivot(matrix.Buffer,
+                                                                            (ulong)(kBlock >> 2),
+                                                                            (ulong)(matrix.KernelRowLength - (kBlock >> 2)),
+                                                                            (ulong)pivots);
+                    int colsInBlock = Math.Min(4, n - kBlock);
+                    for (int offset = 0; offset < colsInBlock; offset++)
+                    {
+                        int k = kBlock + offset;
+                        int pivotRow = k;
+                        float maxAbs;
+
+                        int candidateRow = (int)pivots[offset];
+                        int idx = -1;
+                        for (int r = k; r < n; r++)
+                        {
+                            if (pivot[r] == candidateRow)
+                            {
+                                idx = r;
+                                break;
+                            }
+                        }
+
+                        if (idx != -1)
+                        {
+                            pivotRow = idx;
+                            maxAbs = MathF.Abs(matrix.GetElementUnsafe(pivot[idx], k));
+                        } else
+                        {
+                            maxAbs = MathF.Abs(matrix.GetElementUnsafe(pivot[k], k));
+                            for (int r = k + 1; r < n; r++)
+                            {
+                                float val = MathF.Abs(matrix.GetElementUnsafe(pivot[r], k));
+                                if (val > maxAbs)
+                                {
+                                    maxAbs = val;
+                                    pivotRow = r;
+                                }
+                            }
+                        }
+
+                        if (maxAbs == 0f) {
+                            throw new InvalidOperationException("Zero pivot encountered in PLU.");
+                        }
+                        if (pivotRow != k) {
+                            (pivot[k], pivot[pivotRow]) = (pivot[pivotRow], pivot[k]);
+                        }
+                        lower[k, k] = 1.0f;    // unit diagonal for L
+
+                        // Compute U[k, j] for j >= k
+                        for (int j = k; j < n; j++)
+                        {
+                            float sum = 0f;
+                            for (int s = 0; s < k; s++) {
+                                sum += lower[k, s] * upper[s, j];
+                            }
+                            upper[k, j] = matrix.GetElementUnsafe(pivot[k], j) - sum;
+                        }
+
+                        // Compute L[i, k] for i > k
+                        for (int i = k + 1; i < n; i++)
+                        {
+                            float sum = 0f;
+                            for (int s = 0; s < k; s++) {
+                                sum += lower[i, s] * upper[s, k];
+                            }
+                            float u = upper[k, k];
+                            if (u == 0f) {
+                                throw new InvalidOperationException("Zero pivot encountered in PLU.");
+                            }
+                            lower[i, k] = (matrix.GetElementUnsafe(pivot[i], k) - sum) / u;
+                        }
+                    }
+                }
+
+                // TODO: For larger kernels, consider calling native pivot helper to extract 4 pivots at once.
+            } else if (size < 1024) { } else { }
+        }
+
+    }
+}

src/Nitload.Mathematics/.discard/MatrixCompute.Single.cs.old

@@ -0,0 +1,314 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixCompute.Float"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using Nitload.Common.Collections;
+using Nitload.Common.Internal;
+using Nitload.Mathematics;
+using Nitload.Mathematics.Internal;
+using System;
+using System.Collections.Generic;
+using System.ComponentModel.DataAnnotations.Schema;
+using System.Linq;
+using System.Numerics;
+using System.Security.AccessControl;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static partial class MatrixCompute
+    {
+
+        public static unsafe void Add(
+                                  Matrix<float> left,
+                                  Matrix<float> right,
+                                  Matrix<float> ans
+        )
+        {
+            MatrixElementType f32 = MatrixElementType.F32;
+            if (left.ColumnLength != right.ColumnLength ||
+                left.RowLength != right.RowLength ||
+                left.ColumnLength != ans.ColumnLength ||
+                left.RowLength != ans.RowLength) {
+                throw new IndexOutOfRangeException("Matrix dimensions must agree.");
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            if (count <= 256)
+            {
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                MatrixKernel<float>* rightPtr = right.Buffer;
+                MatrixKernel<float>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_mk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
+                                                                     count);
+
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                MatrixKernel<float>* rightPtr = right.Buffer;
+                MatrixKernel<float>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_mk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
+                                                                     count);
+
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Mat1 = (MatrixKernel*)left.Buffer;
+            package->Mat2 = (MatrixKernel*)right.Buffer;
+            package->Mat3 = (MatrixKernel*)ans.Buffer;
+            package->ElementSize = sizeof(float);
+            package->PrimCount = ans.KernelColumnLength;
+            package->SecCount = ans.KernelRowLength;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_add;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            return;
+        }
+
+        public static unsafe void Fill(
+                                  Matrix<float> ans,
+                                  float value
+        )
+        {
+            MatrixElementType f32 = MatrixElementType.F32;
+            long count = (long)ans.KernelColumnLength * ans.KernelRowLength;
+            ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
+            scalar->Data<float>(0) = value;
+            if (count <= 256)
+            {
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Answer = (nint)ans.Buffer;
+            package->ElementSize = sizeof(float);
+            package->PrimCount = ans.KernelColumnLength;
+            package->SecCount = ans.KernelRowLength;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.ANS_SCALAR_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_fill;
+            package->Scalar = scalar;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            MatrixParallelFactory.Release(scalar);
+            taskGroup.Dispose();
+        }
+
+        public static unsafe void Mul(
+                                  Matrix<float> left,
+                                  Matrix<float> right,
+                                  Matrix<float> ans
+        )
+        {
+            MatrixElementType f32 = MatrixElementType.F32;
+
+            // For matrix multiplication (left: MxK) * (right: KxN) => ans: MxN
+            // ColumnLength = number of columns (N), RowLength = number of rows (M)
+            if (left.ColumnLength != right.RowLength ||
+                ans.ColumnLength != right.ColumnLength ||
+                ans.RowLength != left.RowLength)
+            {
+                throw new IndexOutOfRangeException("Matrix dimensions must agree.");
+            }
+            if (left.Layout != MatrixLayout.DenseRectangle ||
+                right.Layout != MatrixLayout.DenseRectangle ||
+                ans.Layout != MatrixLayout.DenseRectangle) {
+                throw new NotSupportedException("Only dense rectangle matrices are supported.");
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+
+            if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = null;
+                if (count > 256) {
+                    pool = MatrixParallel.LeadParallel(1);
+                }
+                M128* enumeration = stackalloc M128[1];
+                enumeration->Member<int>(2) = left.KernelColumnLength;
+                enumeration->Member<int>(3) = ans.KernelColumnLength;
+                for (int i = 0; i < ans.KernelRowLength; i++)
+                {
+                    for (int j = 0; j < ans.KernelColumnLength; j++)
+                    {
+                        enumeration->Member<int>(0) = i;
+                        enumeration->Member<int>(1) = j;
+
+                        NumericalModule.GlobalContext
+                                           .mat_sk[f32].kernel_tile_fma((MatrixKernel*)left.Buffer,
+                                                                        (MatrixKernel*)right.Buffer,
+                                                                        (MatrixKernel*)ans.Buffer,
+                                                                        enumeration);
+                    }
+                }
+
+                if (count > 256) {
+                    MatrixParallel.ReturnParallelResources(ref pool!);
+                }
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Left = (nint)left.Buffer;
+            package->Right = (nint)right.Buffer;
+            package->Answer = (nint)ans.Buffer;
+            package->AnsPrimOffset = ans.KernelRowLength;
+            package->AnsSecOffset = ans.KernelColumnLength;
+            package->LeftPrimOffset = left.KernelColumnLength;
+            package->PrimCount = ans.KernelColumnLength;
+            package->SecCount = ans.KernelRowLength;
+            package->ElementSize = (int)MatrixElementType.F32;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.SPECIAL, MatrixOperationSpecial.SEQ_FMA);
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            return;
+        }
+
+        public static unsafe void Mul(
+                                  Matrix<float> left,
+                                  float right,
+                                  Matrix<float> ans
+        )
+        {
+            MatrixElementType f32 = MatrixElementType.F32;
+            if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
+                throw new IndexOutOfRangeException("Matrix dimensions must agree.");
+            }
+            if (left.Layout != MatrixLayout.DenseRectangle ||
+                ans.Layout != MatrixLayout.DenseRectangle) {
+                throw new NotSupportedException("Only dense rectangle matrices are supported.");
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
+            scalar->Data<float>(0) = right;
+            if (count <= 256)
+            {
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                MatrixKernel<float>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                                                                            count);
+
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                MatrixKernel<float>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                                                                            count);
+
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Right = (nint)left.Buffer;
+            package->Answer = (nint)ans.Buffer;
+            package->ElementSize = sizeof(float);
+            package->PrimCount = left.KernelColumnLength;
+            package->SecCount = left.KernelRowLength;
+            package->Scalar = scalar;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            MatrixParallelFactory.Release(scalar);
+            return;
+        }
+
+        public static unsafe void Sub(
+                                  Matrix<float> left,
+                                  Matrix<float> right,
+                                  Matrix<float> ans
+        )
+        {
+            MatrixElementType f32 = MatrixElementType.F32;
+            if (left.ColumnLength != right.ColumnLength ||
+                left.RowLength != right.RowLength ||
+                left.ColumnLength != ans.ColumnLength ||
+                left.RowLength != ans.RowLength) {
+                throw new IndexOutOfRangeException("Matrix dimensions must agree.");
+            }
+            if (left.Layout != MatrixLayout.DenseRectangle ||
+                right.Layout != MatrixLayout.DenseRectangle ||
+                ans.Layout != MatrixLayout.DenseRectangle) {
+                throw new NotSupportedException("Only dense rectangle matrices are supported.");
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            if (count <= 256)
+            {
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                MatrixKernel<float>* rightPtr = right.Buffer;
+                MatrixKernel<float>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
+                                                                     count);
+
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+
+                MatrixKernel<float>* leftPtr = left.Buffer;
+                MatrixKernel<float>* rightPtr = right.Buffer;
+                MatrixKernel<float>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
+                                                                     count);
+
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Left = (nint)left.Buffer;
+            package->Right = (nint)right.Buffer;
+            package->Answer = (nint)ans.Buffer;
+            package->ElementSize = sizeof(float);
+            package->PrimCount = left.KernelColumnLength;
+            package->SecCount = left.KernelRowLength;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_sub;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            return;
+        }
+
+    }
+}

src/Nitload.Mathematics/.discard/MatrixCompute.double.cs.old

@@ -0,0 +1,320 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixCompute.double"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using Nitload.Common.Internal;
+using Nitload.Mathematics;
+using Nitload.Mathematics.Internal;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+using unsafe MkProc = delegate* unmanaged[Cdecl]<Nitload.Mathematics.Internal.MatrixParallelTask, void>;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public partial class MatrixCompute
+    {
+
+        public static unsafe void Add(
+                                  Matrix<double> left,
+                                  Matrix<double> right,
+                                  Matrix<double> ans
+        )
+        {
+            MatrixElementType f64 = MatrixElementType.F64;
+            if (left.ColumnLength != right.ColumnLength ||
+                left.RowLength != right.RowLength ||
+                left.ColumnLength != ans.ColumnLength ||
+                left.RowLength != ans.RowLength) {
+                throw new MatrixDimensionMismatchException<double>(
+                    left, right, ans
+                    );
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            if (count <= 256)
+            {
+                MatrixKernel<double>* leftPtr = left.Buffer;
+                MatrixKernel<double>* rightPtr = right.Buffer;
+                MatrixKernel<double>* ansPtr = ans.Buffer;
+                MkProc add = (MkProc)NumericalModule.GlobalContext.
+                    mat_mk[f64.IntrinsicNode].
+                    Get(MatrixIntrinsicHandle.Add);
+
+                for (int i = 0; i < count; i++) {
+                    add(new MatrixParallelTask
+                    {
+                        Mat1 = (MatrixKernel*)(leftPtr + i),
+                        Mat2 = (MatrixKernel*)(rightPtr + i),
+                        Mat3 = (MatrixKernel*)(ansPtr + i),
+                        ElementType = MatrixElementType.F64,
+                    });
+                }
+
+                // .buffer_add(leftPtr, rightPtr, ansPtr, count);
+
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+
+                MatrixKernel<double>* leftPtr = left.Buffer;
+                MatrixKernel<double>* rightPtr = right.Buffer;
+                MatrixKernel<double>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_add(leftPtr, rightPtr, ansPtr,
+                                                                     count);
+
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Left = (nint)left.Buffer;
+            package->Right = (nint)right.Buffer;
+            package->Answer = (nint)ans.Buffer;
+            package->ElementSize = sizeof(double);
+            package->PrimCount = ans.KernelColumnLength;
+            package->SecCount = ans.KernelRowLength;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_add;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            return;
+        }
+
+        public static unsafe void Fill(
+                                  Matrix<double> ans,
+                                  double value
+        )
+        {
+            int f64 = (int)MatrixElementType.F64;
+            long count = (long)ans.KernelColumnLength * ans.KernelRowLength;
+            ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
+            scalar->Data<double>(0) = value;
+            if (count <= 256)
+            {
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Answer = (nint)ans.Buffer;
+            package->ElementSize = sizeof(double);
+            package->PrimCount = ans.KernelColumnLength;
+            package->SecCount = ans.KernelRowLength;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.ANS_SCALAR_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_fill;
+            package->Scalar = scalar;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            MatrixParallelFactory.Release(scalar);
+            taskGroup.Dispose();
+        }
+
+        public static unsafe void Mul(
+                                  Matrix<double> left,
+                                  Matrix<double> right,
+                                  Matrix<double> ans
+        )
+        {
+            int f64 = (int)MatrixElementType.F64;
+
+            // For matrix multiplication (left: MxK) * (right: KxN) => ans: MxN
+            // ColumnLength = number of columns (N), RowLength = number of rows (M)
+            if (left.ColumnLength != right.RowLength ||
+                ans.ColumnLength != right.ColumnLength ||
+                ans.RowLength != left.RowLength)
+            {
+                throw new IndexOutOfRangeException("Matrix dimensions must agree.");
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = null;
+                if (count > 256) {
+                    pool = MatrixParallel.LeadParallel(1);
+                }
+                M128* enumeration = stackalloc M128[1];
+                enumeration->Member<int>(2) = left.KernelColumnLength;
+                enumeration->Member<int>(3) = ans.KernelColumnLength;
+                for (int i = 0; i < ans.KernelRowLength; i++)
+                {
+                    for (int j = 0; j < ans.KernelColumnLength; j++)
+                    {
+                        enumeration->Member<int>(0) = i;
+                        enumeration->Member<int>(1) = j;
+
+                        NumericalModule.GlobalContext
+                                           .mat_sk[f64].kernel_tile_fma((MatrixKernel*)left.Buffer,
+                                                                        (MatrixKernel*)right.Buffer,
+                                                                        (MatrixKernel*)ans.Buffer,
+                                                                        enumeration);
+                    }
+                }
+                if (count > 256) {
+                    MatrixParallel.ReturnParallelResources(ref pool!);
+                }
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Left = (nint)left.Buffer;
+            package->Right = (nint)right.Buffer;
+            package->Answer = (nint)ans.Buffer;
+            package->AnsPrimOffset = ans.KernelRowLength;
+            package->AnsSecOffset = ans.KernelColumnLength;
+            package->LeftPrimOffset = left.KernelColumnLength;
+            package->PrimCount = ans.KernelColumnLength;
+            package->SecCount = ans.KernelRowLength;
+            package->ElementSize = (int)MatrixElementType.F64;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.SPECIAL, MatrixOperationSpecial.SEQ_FMA);
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            return;
+        }
+
+        public static unsafe void Mul(
+                                  Matrix<double> left,
+                                  double right,
+                                  Matrix<double> ans
+        )
+        {
+            int f64 = (int)MatrixElementType.F64;
+            if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
+                throw new IndexOutOfRangeException("Matrix dimensions must agree.");
+            }
+            if (left.Layout != MatrixLayout.DenseRectangle ||
+                ans.Layout != MatrixLayout.DenseRectangle) {
+                throw new NotSupportedException("Only dense rectangle matrices are supported.");
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
+            scalar->Data<double>(0) = right;
+            if (count <= 256)
+            {
+                MatrixKernel<double>* leftPtr = left.Buffer;
+                MatrixKernel<double>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                                                                            count);
+
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+
+                MatrixKernel<double>* leftPtr = left.Buffer;
+                MatrixKernel<double>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
+                                                                            count);
+
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Right = (nint)left.Buffer;
+            package->Answer = (nint)ans.Buffer;
+            package->ElementSize = sizeof(double);
+            package->PrimCount = left.KernelColumnLength;
+            package->SecCount = left.KernelRowLength;
+            package->Scalar = scalar;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            MatrixParallelFactory.Release(scalar);
+            return;
+        }
+
+        public static unsafe void Sub(
+                                  Matrix<double> left,
+                                  Matrix<double> right,
+                                  Matrix<double> ans
+        )
+        {
+            int f64 = (int)MatrixElementType.F64;
+            if (left.ColumnLength != right.ColumnLength ||
+                left.RowLength != right.RowLength ||
+                left.ColumnLength != ans.ColumnLength ||
+                left.RowLength != ans.RowLength) {
+                throw new IndexOutOfRangeException("Matrix dimensions must agree.");
+            }
+            if (left.Layout != MatrixLayout.DenseRectangle ||
+                right.Layout != MatrixLayout.DenseRectangle ||
+                ans.Layout != MatrixLayout.DenseRectangle) {
+                throw new NotSupportedException("Only dense rectangle matrices are supported.");
+            }
+            long count = (long)left.KernelColumnLength * left.KernelRowLength;
+            if (count <= 256)
+            {
+                MatrixKernel<double>* leftPtr = left.Buffer;
+                MatrixKernel<double>* rightPtr = right.Buffer;
+                MatrixKernel<double>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
+                                                                     count);
+
+                return;
+            } else if (count <= 1024)
+            {
+                MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
+
+                MatrixKernel<double>* leftPtr = left.Buffer;
+                MatrixKernel<double>* rightPtr = right.Buffer;
+                MatrixKernel<double>* ansPtr = ans.Buffer;
+                NumericalModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
+                                                                     count);
+
+                MatrixParallel.ReturnParallelResources(ref pool);
+                return;
+            }
+            MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
+            package->Left = (nint)left.Buffer;
+            package->Right = (nint)right.Buffer;
+            package->Answer = (nint)ans.Buffer;
+            package->ElementSize = sizeof(double);
+            package->PrimCount = left.KernelColumnLength;
+            package->SecCount = left.KernelRowLength;
+            package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
+            package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_sub;
+            MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
+
+            taskGroup.LaunchParallel(SleepMode.DeepSleep);
+            taskGroup.Dispose();
+            return;
+        }
+
+    }
+}

src/Nitload.Mathematics/.discard/MatrixParallelHandle.cs.old

@@ -0,0 +1,197 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixParallelHandle"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using Nitload.Common.Collections;
+using System;
+using System.Collections.Generic;
+using System.ComponentModel.DataAnnotations;
+using System.Diagnostics.CodeAnalysis;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public class MatrixParallelHandle : IDisposable
+    {
+
+        private static readonly (int a, int b)[] LookupTable = [
+            (1, 1), // 7
+            (0, 2), // 8  (no non-zero solution)
+            (3, 0), // 9  (no non-zero solution)
+            (2, 1), // 10
+            (1, 2), // 11
+            (0, 3), // 12 (no non-zero solution)
+            (3, 1), // 13
+            (2, 2), // 14
+            (1, 3), // 15
+            (4, 1), // 16
+            (3, 2), // 17
+            (2, 3), // 18
+            (1, 4), // 19
+            (4, 2), // 20
+        ];
+        private MatrixParallelThread[] _threads;
+        private MatrixParallelWrap[] _wraps;
+        private bool _disposed = false;
+        private int _count;
+        private ManualResetEventSlim _leaderEvent = new ManualResetEventSlim(false);
+
+        internal MatrixParallelHandle(MatrixParallelThread[] threadArray)
+        {
+            _threads = threadArray;
+            IsMultithread = threadArray.Length > 1;
+            SliceToWraps();
+        }
+
+        public ReadOnlySpan<MatrixParallelThread> Threads
+        {
+            get
+            {
+                if (_disposed) {
+                    throw new ObjectDisposedException(typeof(MatrixParallelHandle).Name);
+                }
+                return _threads;
+            }
+        }
+
+        internal ReadOnlySpan<MatrixParallelWrap> Wraps
+        {
+            get
+            {
+                if (_disposed) {
+                    throw new ObjectDisposedException(typeof(MatrixParallelHandle).Name);
+                }
+                return _wraps;
+            }
+        }
+
+        private bool IsMultithread { get; init; }
+
+        public void Dispose()
+        {
+            if (_disposed) {
+                return;
+            }
+            MatrixParallel.ReturnParallelResources(ref _threads);
+            _disposed = true;
+        }
+
+        public MatrixParallelThread GetLeader()
+        {
+            MatrixParallelThread leader = Threads[1];
+            leader.Role = MatrixParallelThread.RoleType.Leader;
+            leader.LeaderEvent = _leaderEvent;
+            return leader;
+        }
+
+        public void WaitAllWraps()
+        {
+            if (_disposed) {
+                throw new ObjectDisposedException(typeof(MatrixParallelHandle).Name);
+            }
+            SpinWait.SpinUntil(() => Volatile.Read(ref _count) == 0);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal void ReportWrapAccomplished()
+        {
+            if (_disposed) {
+                throw new ObjectDisposedException(typeof(MatrixParallelHandle).Name);
+            }
+            Interlocked.Decrement(ref _count);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static (int count_3, int count_4) FastDecompose(int c)
+        {
+            /*
+            * Thread wrapping to accelerate matrix operations by dividing threads into groups (pool).
+            * For M pool, each containing an average of N threads, we satisfy the condition M * N = T - 4,
+            * where T is the total thread _count of the CPU. Typically, for modern consumer CPUs, T is around 16.
+            * For simplicity, let M * N = 16.
+            *
+            * The total time complexity for publishing tasks is O(M + N):
+            * - O(M) accounts for distributing tasks across M pool.
+            * - O(N) accounts for threads within each wrap fetching tasks from the wrap queue.
+            *
+            * Using the inequality Sqrt(M * N) <= (M + N) / 2 (derived from the arithmetic-geometric mean inequality),
+            * the optimal configuration occurs when M = N. For M * N = 16, this gives M = N = 4.
+            * 
+            * To ensure compatibility across a wide range of CPU architectures, 
+            * we allow each wrap to contain either 3 or 4 threads, providing flexibility while maintaining efficiency.
+            */
+            if (c < 7)
+            {
+                return (0, 0);
+            }
+
+            int t = c - 7;
+            int k = t / 14;
+            int d = t % 14;
+
+            (int a0, int b0) = LookupTable[d];
+
+            return (a0 + (2 * k), b0 + (2 * k));
+        }
+
+        private void SliceToWraps()
+        {
+            if (_disposed) {
+                throw new ObjectDisposedException(typeof(MatrixParallelHandle).Name);
+            }
+            if (_wraps is not null) {
+                return;
+            }
+            int count = _threads.Length;
+            if (count < 6)
+            {
+                MatrixParallelWrap wrap = new MatrixParallelWrap(this, 0, count);
+                _wraps = [ wrap ];
+            } else
+            {
+                (int c3, int c4) = FastDecompose(_threads.Length);
+                _wraps = new MatrixParallelWrap[c3 + c4];
+                int begin = 0;
+                int i = 0;
+                for (; i < c3; i++)
+                {
+                    _wraps[i] = new MatrixParallelWrap(this, begin, 3);
+                    begin += 3;
+                }
+                for (; i < c3 + c4; i++)
+                {
+                    _wraps[i] = new MatrixParallelWrap(this, begin, 4);
+                    begin += 4;
+                }
+            }
+            _count = _wraps.Length;
+        }
+
+    }
+}

src/Nitload.Mathematics/.discard/MatrixParallelWrap.cs.old

@@ -0,0 +1,189 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixParallelWrap"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using Nitload.Common.HighPerformance;
+using System;
+using System.Collections.Generic;
+using System.Diagnostics.CodeAnalysis;
+using System.Linq;
+using System.Net.WebSockets;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public unsafe class MatrixParallelWrap
+    {
+
+        private int _begin;
+
+        private int _remainThreadCount;
+        private readonly MatrixParallelHandle _handle;
+
+        internal MatrixParallelWrap(
+                 MatrixParallelHandle handle,
+                 int begin,
+                 int count
+        )
+        {
+            WrapTaskCapacity = count;
+            _handle = handle;
+            _begin = begin;
+            for (int i = 0; i < WrapTaskCapacity; i++)
+            {
+                Threads[i].Role = MatrixParallelThread.RoleType.Worker;
+                Threads[i].BindedWrap = this;
+            }
+            Threads[0].Role = MatrixParallelThread.RoleType.Scheduler;
+        }
+
+        public static SleepMode AfterWork
+        {
+            get;
+            set => Interlocked.Exchange(ref Unsafe.As<SleepMode, int>(ref field), (int)value);
+        }
+
+        internal int Magic { get; set; }
+
+        private ReadOnlySpan<MatrixParallelThread> Threads
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get => _handle.Threads[_begin..(_begin + WrapTaskCapacity)];
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void ReportThreadAccomplished()
+        {
+            Interlocked.Decrement(ref _remainThreadCount);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void ReportWrapAccomplished()
+        {
+            _handle.ReportWrapAccomplished();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void SetAllWorker()
+        {
+            for (int i = 1; i < WrapTaskCapacity; i++) {
+                Threads[i].Restart();
+            }
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void WaitAllWorker()
+        {
+            SpinWait.SpinUntil(() => Volatile.Read(ref _remainThreadCount) == 0);
+        }
+
+        internal void PublishTaskToWorkers()
+        {
+            Queue<nint> tasks = WrapTask;
+            while (tasks.TryDequeue(out nint handle))
+            {
+                MatrixParallelTaskPackage* task = (MatrixParallelTaskPackage*)handle;
+                if ((int)task->ComputeMode.OperationStyle > 0)
+                {
+                    PublishBufferTaskToWorkers(task);
+                } else
+                {
+                    throw new NotImplementedException();
+                }
+            }
+        }
+
+        internal void SetScheduler()
+        {
+            MatrixParallelThread thread = Threads[0];
+            thread.Restart();
+            _ = Interlocked.Exchange(ref _remainThreadCount, WrapTaskCapacity);
+        }
+
+        private void PublishBufferTaskToWorkers(
+                     MatrixParallelTaskPackage* package
+        )
+        {
+            ReadOnlySpan<MatrixParallelThread> pool = Threads;
+            if (pool.Length == 1)
+            {
+                pool[0].Tasks.Enqueue((nint)package);
+                return;
+            }
+            long count = Unsafe.As<int, long>(ref package->PrimCount);
+            long actual = count / pool.Length;
+            long extra = count % pool.Length;
+
+            long offset = 0;
+            MatrixParallelTaskPackage* current = package;
+            /*
+            string msg = $"""
+                Total: count={count}, length={Threads.Length}{Environment.NewLine}
+                """;
+            /**/
+            foreach (MatrixParallelThread thread in pool)
+            {
+                long wrapTaskLength = actual;
+                long ex = extra > 0 ? 1 : 0;
+                long length = wrapTaskLength + ex;
+                extra -= ex;
+                ref long offsetRef = ref Unsafe.As<int, long>(ref current->PrimTileOffset);
+                offsetRef += offset;
+                ref long lengthRef = ref Unsafe.As<int, long>(ref current->PrimCount);
+                lengthRef = length;
+                thread.Tasks.Enqueue((nint)current);
+                /*
+                msg += $"""
+                    Slice to thread: total offset={offsetRef} local offset={offset}, length={lengthRef}{Environment.NewLine}
+                    """;
+                /**/
+                current = MatrixParallelFactory.FromExistPackage(current);
+                offset = length;
+            }
+
+            // Console.WriteLine(msg);
+        }
+
+        #region wrap task
+
+        internal unsafe Queue<nint> WrapTask { get; set; } = [];
+
+        internal int WrapTaskCapacity { get ; init; }
+
+        #endregion
+    }
+
+    public enum SleepMode : int
+    {
+
+        QuickWakeup = ThreadPriority.AboveNormal,
+        NormalSleep = ThreadPriority.BelowNormal,
+        DeepSleep = ThreadPriority.Lowest,
+
+    }
+}

src/Nitload.Mathematics/AssemblyInfo.cs

@@ -0,0 +1,17 @@
+// -----------------------------------------------------------------------------
+// file="AssemblyInfo"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+global using Nitload.Common.Internal;
+global using Nitload.Mathematics.Internal;
+global using System;
+global using System.Numerics;
+
+global using System.Threading;
+using System.Runtime.CompilerServices;
+
+[assembly: InternalsVisibleTo("Nitload.Graphics")]
+[assembly: InternalsVisibleTo("StgSharpDebug")]

src/Nitload.Mathematics/Geometries/Coordination/CoordinationBase.cs

@@ -0,0 +1,56 @@
+// -----------------------------------------------------------------------------
+// file="CoordinationBase"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System.Runtime.CompilerServices;
+
+namespace Nitload.Geometries
+{
+    public abstract class CoordinationBase
+    {
+
+        private bool globalCoordAvailable;
+        private CoordinationBase localCoordination;
+        private GraphicsMatrix coordMat;
+
+        public CoordinationBase(
+               CoordinationBase localCoordination
+        )
+        {
+            this.localCoordination = localCoordination;
+            coordMat = new GraphicsMatrix();
+        }
+
+        public virtual Point LocalOrigin
+        {
+            get => new Point(coordMat.column[3].XYZ);
+            internal set => coordMat.column[3].XYZ = value.Coord;
+        }
+
+        public virtual Vec3 LocalX
+        {
+            get => coordMat.column[3].XYZ;
+            internal set => coordMat.column[3].XYZ = value;
+        }
+
+        public virtual Vec3 LocalY
+        {
+            get => coordMat.column[3].XYZ;
+            internal set => coordMat.column[3].XYZ = value;
+        }
+
+        public virtual Vec3 LocalZ
+        {
+            get => coordMat.column[3].XYZ;
+            set => coordMat.column[3].XYZ = value;
+        }
+
+        protected ref GraphicsMatrix CoordMat => ref coordMat;
+
+    }
+}

src/Nitload.Mathematics/Geometries/Coordination/PlainCoordinate.cs

@@ -0,0 +1,44 @@
+// -----------------------------------------------------------------------------
+// file="PlainCoordinate"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Collections.Generic;
+using System.Net.Http.Headers;
+using System.Text;
+
+namespace Nitload.Geometries
+{
+    public class PlainCoordinate : CoordinationBase
+    {
+
+        public PlainCoordinate(
+               CoordinationBase local,
+               Point origin,
+               Vec3 xAxis,
+               Vec3 yAxis
+        )
+            : base(local)
+        {
+            CoordMat = GraphicsMatrix.Unit;
+
+            CoordMat.column[0] = new Vec4(xAxis, 0);
+            CoordMat.column[1] = new Vec4(yAxis, 0);
+            CoordMat.column[2] = new Vec4(Linear.Normalize(Vec3.Cross(xAxis, yAxis)), 0);
+
+            LocalOrigin = origin;
+        }
+
+        public static PlainCoordinate StandardPlainCoordination { get; } = new(null,
+                                                                               new Point(0, 0, 0),
+                                                                               new Vec3(1, 0, 0),
+                                                                               new Vec3(0, 1, 0));
+
+    }
+}

src/Nitload.Mathematics/Geometries/Coordination/SpaceCoordinate.cs

@@ -0,0 +1,30 @@
+// -----------------------------------------------------------------------------
+// file="SpaceCoordinate"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace Nitload.Geometries
+{
+    public class SpaceCoordinate : CoordinationBase
+    {
+
+        public SpaceCoordinate(
+               CoordinationBase localCoordination,
+               Point origin,
+               Vec3 xAxis,
+               Vec3 yAxis,
+               Vec3 zAxis
+        )
+            : base(localCoordination) { }
+
+    }
+}

src/Nitload.Mathematics/Geometries/Coordination/StreamLocalCoordinate.cs

@@ -0,0 +1,23 @@
+// -----------------------------------------------------------------------------
+// file="StreamLocalCoordinate"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    /*
+    public class StreamLocalCoordinate : CoordinationBase
+    {
+    }
+    /**/
+}

src/Nitload.Mathematics/Geometries/GemoryException.cs

@@ -0,0 +1,63 @@
+// -----------------------------------------------------------------------------
+// file="GemoryException"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Runtime.Serialization;
+
+namespace Nitload.Geometries
+{
+    public unsafe class ToUpperDimensionException : Exception
+    {
+
+        protected ToUpperDimensionException(
+                  SerializationInfo info,
+                  StreamingContext context
+        )
+            : base(info, context) { }
+
+        public ToUpperDimensionException() { }
+
+        public ToUpperDimensionException(
+               string message
+        )
+            : base(message) { }
+
+        public ToUpperDimensionException(
+               string message,
+               Exception innerException
+        )
+            : base(message, innerException) { }
+
+    }
+
+    /// <summary>
+    ///   当几何体调用了过多的CalcVec方法时，产生UnusedVertexException异常
+    /// </summary>
+    public class UnusedVertexException : Exception
+    {
+
+        protected UnusedVertexException(
+                  SerializationInfo info,
+                  StreamingContext context
+        )
+            : base(info, context) { }
+
+        public UnusedVertexException() { }
+
+        public UnusedVertexException(
+               string message
+        )
+            : base(message) { }
+
+        public UnusedVertexException(
+               string message,
+               Exception innerException
+        )
+            : base(message, innerException) { }
+
+    }
+}

src/Nitload.Mathematics/Geometries/GeometryOperation/GeometryOperation.cs

@@ -0,0 +1,97 @@
+// -----------------------------------------------------------------------------
+// file="GeometryOperation"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.Timing;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Geometries
+{
+    [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
+    public delegate Vec3<float> GeometryMotionDelegate(
+                                TimeSpanProvider timeSource,
+                                out float rotation
+    );
+
+    public class GeometryMotion
+    {
+
+        private bool isIncrementMode;
+        private bool sourceIsStatic;
+        private GeometryMotionDelegate movement;
+        private IntPtr handle;
+        private TimeSpanProvider time;
+
+        public unsafe GeometryMotion(
+                      bool isIncreament,
+                      TimeSpanProvider timeSource,
+                      GeometryMotionDelegate motionDelegate
+        )
+        {
+            isIncrementMode = isIncreament;
+            time = timeSource;
+            if (motionDelegate == null)
+            {
+                handle = (IntPtr)(delegate*<TimeSpanProvider, out float, Vec3<float>>)&defaultMotion;
+                movement = defaultMotion;
+            } else if (motionDelegate.Method.IsStatic)
+            {
+                sourceIsStatic = true;
+                handle = Marshal.GetFunctionPointerForDelegate(motionDelegate);
+                movement = motionDelegate;
+            } else
+            {
+                movement = motionDelegate;
+            }
+
+            // sourceIsStatic = false;
+        }
+
+        public bool IsIncrement => isIncrementMode;
+
+        public static GeometryMotion DefaultMotion => new GeometryMotion(true, null, null);
+
+        public unsafe void RunMotion(
+                           ref Vec3<float> coord,
+                           ref float rotation
+        )
+        {
+            float rot;
+            Vec3<float> vec;
+            if (sourceIsStatic)
+            {
+                vec = ((delegate* unmanaged[Cdecl]<TimeSpanProvider, out float, Vec3<float>>)handle)(
+                    time, out rot);
+            } else
+            {
+                vec = movement(time, out rot);
+            }
+
+            if (isIncrementMode)
+            {
+                coord += vec;
+                rotation += rot;
+            } else
+            {
+                coord = vec;
+                rotation = rot;
+            }
+        }
+
+        private static Vec3<float> defaultMotion(
+                                   TimeSpanProvider provider,
+                                   out float rotation
+        )
+        {
+            rotation = 0;
+            return Vec3<float>.Zero;
+        }
+
+    }//------------------------------- End of Class -------------------------------
+}

src/Nitload.Mathematics/Geometries/GeometryScaler/Angle.cs

@@ -0,0 +1,197 @@
+// -----------------------------------------------------------------------------
+// file="Angle"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics
+{
+    public static partial class GeometryScaler
+    {
+
+        internal const float DegToRad = 180 / MathF.PI;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float Cos(
+                            Radius r
+        )
+        {
+            return MathF.Cos(r._radius);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float Cos(
+                            Degree r
+        )
+        {
+            return MathF.Cos(r._degree / DegToRad);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float Sin(
+                            Radius r
+        )
+        {
+            return MathF.Sin(r._radius);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float Sin(
+                            Degree r
+        )
+        {
+            return MathF.Sin(r._degree / DegToRad);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float Tan(
+                            Radius r
+        )
+        {
+            return MathF.Tan(r._radius);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float Tan(
+                            Degree r
+        )
+        {
+            return MathF.Tan(r._degree / DegToRad);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Degree ToDegree(
+                             float value
+        )
+        {
+            return new Degree(value);
+        }
+
+    }
+
+    public struct Degree
+    {
+
+        internal float _degree;
+
+        public Degree(
+               float dgree
+        )
+        {
+            _degree = dgree;
+        }
+
+        public static implicit operator Radius(
+                                        Degree d
+        )
+        {
+            return new Radius(d._degree / GeometryScaler.DegToRad);
+        }
+
+    }
+
+    public struct Radius
+    {
+
+        internal float _radius;
+
+        [Obsolete("Unsafe, you may forget this is radius.", true)]
+        public Radius() { }
+
+        public Radius(
+               float radius
+        )
+        {
+            _radius = radius;
+        }
+
+        public static Radius Zero => new Radius(0);
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Radius operator -(
+                                      Radius r
+        )
+        {
+            return new Radius(-r._radius);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Radius operator -(
+                                      Radius left,
+                                      Radius Right
+        )
+        {
+            return new Radius(left._radius - Right._radius);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator !=(
+                                    Radius left,
+                                    Radius right
+        )
+        {
+            return left._radius != right._radius;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Radius operator *(
+                                      Radius left,
+                                      float right
+        )
+        {
+            return new Radius(left._radius * right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Radius operator /(
+                                      Radius left,
+                                      float right
+        )
+        {
+            return new Radius(left._radius / right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float operator /(
+                                     Radius left,
+                                     Radius right
+        )
+        {
+            return left._radius / right._radius;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Radius operator +(
+                                      Radius left,
+                                      Radius Right
+        )
+        {
+            return new Radius(left._radius + Right._radius);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator ==(
+                                    Radius left,
+                                    Radius right
+        )
+        {
+            return left._radius == right._radius;
+        }
+
+        public static implicit operator Degree(
+                                        Radius r
+        )
+        {
+            return new Degree(r._radius * 180 / MathF.PI);
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/IGeometry.cs

@@ -0,0 +1,38 @@
+// -----------------------------------------------------------------------------
+// file="IGeometry"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Collections.Generic;
+using System.ComponentModel;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    public interface IGeometry
+    {
+
+        public Point this[
+                     int index
+        ] { get; set; }
+
+        #pragma warning disable CA1819 
+        public Vec4[] VertexStream { get; }
+ #pragma warning restore CA1819 
+
+        public CoordinationBase Coordination { get; }
+
+        public int VertexCount { get; }
+
+        public ReadOnlySpan<int> VertexIndices { get; }
+
+    }
+}

src/Nitload.Mathematics/Geometries/IGeomtryBuffer.cs

@@ -0,0 +1,17 @@
+// -----------------------------------------------------------------------------
+// file="IGeomtryBuffer"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    public interface IGeomtryBuffer { }
+}

src/Nitload.Mathematics/Geometries/IInstancing.cs

@@ -0,0 +1,55 @@
+// -----------------------------------------------------------------------------
+// file="IInstancing"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Numerics;
+using System.Security.Cryptography.X509Certificates;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    public interface IInstancing
+    {
+
+        public float Rotation //
+        { get => GlobalBuffer.CoordAndRotationList[BufferId].W; //
+        set
+        {
+            Vec4<float> temp = GlobalBuffer.CoordAndRotationList[BufferId];
+            temp.W = value;
+            GlobalBuffer.CoordAndRotationList[BufferId] = temp;
+        } }
+
+        public float Scale { get => GlobalBuffer.ScalingList[BufferId]; set => GlobalBuffer.ScalingList[BufferId] =
+                                                                               value; }
+
+        public GeometryMotion Motion { get; set; }
+
+        public IInstancingBuffer GlobalBuffer { get; internal set; }
+
+        public int BufferId { get; internal set; }
+
+        public Vec3<float> CenterPositionGlobal { get; internal set; }
+
+        public Vec3<float> Coord { get => GlobalBuffer.CoordAndRotationList[BufferId].XYZ; set => GlobalBuffer.CoordAndRotationList[BufferId] =
+                                                                                                  new Vec4<float>(value.reg); }
+
+        public void Move()
+        {
+            (Vec3<float> coord, float rotation) = (Coord, Rotation);
+            Motion.RunMotion(ref coord, ref rotation);
+            (Coord, Rotation) = (coord, rotation);
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/IInstancingBuffer.cs

@@ -0,0 +1,33 @@
+// -----------------------------------------------------------------------------
+// file="IInstancingBuffer"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace Nitload.Geometries
+{
+    public interface IInstancingBuffer
+    {
+
+        public Span<Vec4<float>> CoordAndRotationSpan { get; }
+
+        public Span<float> ScalingSpan { get; }
+
+        internal IGeometry TypicalShape { get; }
+
+        internal List<IInstancing> InstanceList { get; }
+
+        internal List<Vec4<float>> CoordAndRotationList { get; }
+
+        internal List<float> ScalingList { get; }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Cordinate.cs

@@ -0,0 +1,58 @@
+// -----------------------------------------------------------------------------
+// file="Cordinate"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    [StructLayout(LayoutKind.Sequential)]
+    public struct PartialPlainCoordinate
+    {
+        /*
+        internal Matrix32 axis;
+
+        /// <summary>
+        ///   The _origin of this coordinate.
+        /// </summary>
+        public Vec3 origin;
+
+        public PartialPlainCoordinate( Vec3 origin, Vec3 horiziontial, Vec3 vertical )
+        {
+            this.origin = origin;
+            axis = new Matrix32();
+            Unsafe.SkipInit( out axis.mat );
+            axis.colum0 = horiziontial;
+            axis.colum1 = vertical;
+        }
+
+        /// <summary>
+        ///   The horizontial axis of the coordinate. Point to the right direction.
+        /// </summary>
+        public Vec3 axis_H => axis.colum0;
+
+        /// <summary>
+        ///   The vertical axis of the coordinate. Point to the top direction.
+        /// </summary>
+        public Vec3 axis_V => axis.colum1;
+
+        [MethodImpl( MethodImplOptions.AggressiveInlining )]
+        public Vec3 GetRealLocate( Vec2 raletiveLocation )
+        {
+            return origin + axis.VerticalVecMultiply( raletiveLocation );
+        }
+        */
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Define/GeometryOperation.cs

@@ -0,0 +1,27 @@
+// -----------------------------------------------------------------------------
+// file="GeometryOperation"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System.Runtime.CompilerServices;
+
+namespace Nitload.Geometries
+{
+    public static unsafe class GeometryOperation
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3 DefaultMotion(
+                           int tick
+        )
+        {
+            return default(Vec3);
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Define/IPlainEntity.cs

@@ -0,0 +1,55 @@
+// -----------------------------------------------------------------------------
+// file="IPlainEntity"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    public interface IPlainEntity
+    {
+
+        public Vec2 CenterPosition { get; set; }
+
+        public bool CollideWith(
+                    IPlainEntity entity
+        );
+
+    }
+
+    public class PlainEntity : IPlainEntity
+    {
+
+        private PlainGeometry _shape;
+
+        public PlainEntity(
+               PlainGeometry shape
+        )
+        {
+            _shape = shape;
+        }
+
+        bool IPlainEntity.CollideWith(
+                          IPlainEntity entity
+        )
+        {
+            throw new NotImplementedException();
+        }
+
+        Vec2 IPlainEntity.CenterPosition
+        {
+            get => _shape.Coordination.LocalOrigin.Coord.XY;
+            set => throw new InvalidOperationException();
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Define/ParticleOnPlain.cs

@@ -0,0 +1,58 @@
+// -----------------------------------------------------------------------------
+// file="ParticleOnPlain"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+using Nitload.Mathematics.Numeric.Graphics;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Numerics;
+using System.Security.Cryptography.X509Certificates;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    public class ParticleOnPlain : IInstancing, IPlainEntity
+    {
+
+        public ParticleOnPlain(
+               PlainInstancingBuffer<ParticleOnPlain> buffer
+        )
+        {
+            if (buffer == null) {
+                throw new ArgumentNullException(paramName:nameof(buffer));
+            }
+            GlobalBuffer = buffer;
+            BufferId = buffer.CreateInstanceID();
+            ((IInstancingBuffer)buffer).InstanceList.Add(this);
+        }
+
+        public GeometryMotion Motion { get; set; }
+
+        public IInstancingBuffer GlobalBuffer { get; set; }
+
+        public int BufferId { get; set; }
+
+        Vec3 IInstancing.CenterPositionGlobal { get; set; }
+
+        bool IPlainEntity.CollideWith(
+                          IPlainEntity entity
+        )
+        {
+            throw new NotImplementedException();
+        }
+
+        Vec2 IPlainEntity.CenterPosition
+        {
+            get => GlobalBuffer.CoordAndRotationList[BufferId].XY;
+            set => throw new NotImplementedException();
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Define/PlainGeometry.cs

@@ -0,0 +1,51 @@
+// -----------------------------------------------------------------------------
+// file="PlainGeometry"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload;
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.ComponentModel;
+
+namespace Nitload.Geometries
+{
+    /// <summary>
+    ///   Plain geometry only contains straight sides.
+    /// </summary>
+    public abstract class PlainGeometry : IGeometry
+    {
+
+        protected Vec4[] vertexMat;
+        internal CoordinationBase coordinate;
+
+        internal PlainGeometry(
+                 CoordinationBase coordination
+        )
+        {
+            coordinate = coordination;
+        }
+
+        public Point this[
+                     int index
+        ]
+        {
+            get => new Point(vertexMat[index]);
+            set => vertexMat[index].vec = value.coordVec;
+        }
+
+        public int VertexCount => this.vertexMat.Length;
+
+        public abstract ReadOnlySpan<int> VertexIndices { get; }
+
+        #pragma warning disable CA1819 
+        public Vec4[] VertexStream => vertexMat;
+
+        public CoordinationBase Coordination => this.coordinate;
+#pragma warning restore CA1819
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Define/PlainInstancingBuffer.cs

@@ -0,0 +1,151 @@
+// -----------------------------------------------------------------------------
+// file="PlainInstancingBuffer"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+using Nitload.Mathematics.Numeric.Graphics;
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using System.ComponentModel;
+using System.Diagnostics.Contracts;
+using System.Linq;
+using System.Numerics;
+using System.Reflection;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Text;
+using System.Threading;
+using System.Threading.Tasks;
+
+namespace Nitload.Geometries
+{
+    public class PlainInstancingBuffer<T> : IDisposable, IEnumerable<T>, IInstancingBuffer
+        where T : IInstancing
+    {
+
+        private bool disposedValue;
+
+        private CoordinatePlain globalCoordinate;
+        private List<Vec4> coordList;
+        private List<IInstancing> instanceList;
+        private List<float> scalingList;
+        private Mutex mutex;
+
+        public PlainInstancingBuffer()
+        {
+            coordList = new List<Vec4>(capacity:100);
+            instanceList = new List<IInstancing>(capacity:100);
+            scalingList = new List<float>(capacity:100);
+            coordArray = Array.Empty<Vec4>();
+            scaleArray = Array.Empty<float>();
+            mutex = new Mutex();
+        }
+
+        public List<IInstancing> InstanceList => instanceList;
+
+        public List<Vec4> CoordAndRotationList => coordList;
+
+        public List<float> ScalingList => scalingList;
+
+        public Span<Vec4> CoordAndRotationSpan => CollectionsMarshal.AsSpan(coordList);
+
+        public Span<float> ScalingSpan => CollectionsMarshal.AsSpan(scalingList);
+
+        public int CreateInstanceID()
+        {
+            mutex.WaitOne();
+            int ret = instanceList.Count;
+            coordList.Add(new Vec4());
+            scalingList.Add(0);
+            mutex.ReleaseMutex();
+            return ret;
+        }
+
+        // ~PlainInstancingBuffer()
+        // {
+        // Dispose(disposing: false);
+        // }
+
+        public void Dispose()
+        {
+            Dispose(disposing:true);
+            GC.SuppressFinalize(this);
+        }
+
+        public IEnumerator<T> GetEnumerator()
+        {
+            foreach (IInstancing item in instanceList) {
+                yield return (T)item;
+            }
+        }
+
+        public void Remove(
+                    T instance
+        )
+        {
+            mutex.WaitOne();
+            if (instance.GlobalBuffer != this) {
+                throw new ArgumentException(
+                    message:"instance does not belong to current Particle Buffer",
+                    paramName:nameof(instance));
+            }
+            int index = instance.BufferId;
+            IInstancing temp = instanceList[^1];
+            instanceList[instance.BufferId] = temp;
+            temp.BufferId = instance.BufferId;
+
+            coordList[instance.BufferId] = coordList[^1];
+            scalingList[instance.BufferId] = scalingList[^1];
+
+            instanceList.RemoveAt(instanceList.Count - 1);
+            coordList.RemoveAt(coordList.Count - 1);
+            scalingList.RemoveAt(scalingList.Count - 1);
+            mutex.ReleaseMutex();
+        }
+
+        public void UpdateAllParticles()
+        {
+            foreach (IInstancing particle in instanceList) {
+                particle.Move();
+            }
+        }
+
+        protected virtual void Dispose(
+                               bool disposing
+        )
+        {
+            if (!disposedValue)
+            {
+                if (disposing) {
+                    mutex.Dispose();
+                }
+
+                coordList.Clear();
+                instanceList.Clear();
+                scalingList.Clear();
+                disposedValue = true;
+            }
+        }
+
+        IEnumerator IEnumerable.GetEnumerator()
+        {
+            foreach (IInstancing item in instanceList) {
+                yield return (T)item;
+            }
+        }
+
+        IGeometry IInstancingBuffer.TypicalShape => throw new NotImplementedException();
+
+        #region field of mem opt
+
+        private Vec4[] coordArray;
+        private float[] scaleArray;
+
+    #endregion
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Parallelogram.cs

@@ -0,0 +1,68 @@
+// -----------------------------------------------------------------------------
+// file="Parallelogram"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+
+namespace Nitload.Geometries
+{
+    public class Parallelogram : Polygon4
+    {
+
+        internal readonly Func<int, Vec3> movCenterOperation = GeometryOperation.DefaultMotion;
+
+        internal Point center;
+
+        public Parallelogram(
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y
+        )
+            : base(v0x, v0y, v1x, v1y, v2x, v2y, v3x, v3y)
+        {
+            #if DEBUG
+            if (this[0].Coord + this[2].Coord != this[1].Coord + this[3].Coord)
+            {
+                DefaultLog.InternalWriteLog(
+                    "Init of geometry item failed, because four vertices cannot form a rectangle.",
+                    LogType.Warning);
+            }
+            #endif
+        }
+
+        public Parallelogram(
+               PlainCoordinate coordination,
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y
+        )
+            : base(coordination, v0x, v0y, v1x, v1y, v2x, v2y, v3x, v3y)
+        {
+            #if DEBUG
+            if (this[0].Coord + this[2].Coord != this[1].Coord + this[3].Coord)
+            {
+                DefaultLog.InternalWriteLog(
+                    "Init of geometry item failed, because four vertices cannot form a rectangle.",
+                    LogType.Warning);
+            }
+            #endif
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Plain.cs

@@ -0,0 +1,74 @@
+// -----------------------------------------------------------------------------
+// file="Plain"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System.Numerics;
+
+namespace Nitload.Geometries
+{
+    public unsafe class Plain
+    {
+
+        internal readonly float d = -1;
+        internal Vec3 _origin;
+
+        internal Vec3 plainParameter;
+        internal Vec3 vec1;
+        internal Vec3 vec2;
+
+        public Plain(
+               Vec3 o,
+               Vec3 v1,
+               Vec3 v2
+        )
+        {
+            this._origin = o;
+            this.vec1 = v1 - (v2 * (v1.Y / v2.Y));   //xz
+            this.vec2 = v2 - (vec1 * (v2.X / v1.X)); //yz
+
+            // |_origin|-k1*|vec1|-k2*|vec2| --> plainParameter.Z；
+            plainParameter.Z = _origin.Z - (vec1.Z * (_origin.X / vec1.X)) - (vec2.Z * (_origin.Y / vec2.Y));
+
+            // plainParameter.Z -> plainParameter.X, plainParameter.Y
+            plainParameter.X = vec1.X * (vec1.Z / plainParameter.Z);
+            plainParameter.Y = vec2.Y * (vec2.Z / plainParameter.Z);
+
+            plainParameter.X = 1 / plainParameter.X;
+            plainParameter.Y = 1 / plainParameter.Y;
+            plainParameter.Z = 1 / plainParameter.Z;
+        }
+
+        public Plain(
+               Point p1,
+               Point p2,
+               Point p3
+        )
+        {
+            Vec3 origin = p1.Coord;
+            Vec3 v1 = p2.Coord - origin;
+            Vec3 v2 = p3.Coord - origin;
+
+
+            vec1 = v1 - (v2 * (v1.Y / v2.Y));   //xz
+            vec2 = v2 - (vec1 * (v2.X / v1.X)); //yz
+
+            // |_origin|-k1*|vec1|-k2*|vec2| --> plainParameter.Z；
+            plainParameter.Z = origin.Z - (vec1.Z * (origin.X / vec1.X)) - (vec2.Z * (origin.Y / vec2.Y));
+
+            // plainParameter.Z --> plainParameter.X, plainParameter.Y
+            plainParameter.X = vec1.X * (vec1.Z / plainParameter.Z);
+            plainParameter.Y = vec2.Y * (vec2.Z / plainParameter.Z);
+
+            plainParameter.X = 1 / plainParameter.X;
+            plainParameter.Y = 1 / plainParameter.Y;
+            plainParameter.Z = 1 / plainParameter.Z;
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Polygon4.cs

@@ -0,0 +1,119 @@
+// -----------------------------------------------------------------------------
+// file="Polygon4"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Numerics;
+
+namespace Nitload.Geometries
+{
+    public class Polygon4 : PlainGeometry
+    {
+
+        protected static int[] Indices = [
+            0,
+            1,
+            2,
+            0,
+            2,
+            3
+        ];
+
+        internal Polygon4()
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            vertexMat = new Vec4[4];
+        }
+
+        internal Polygon4(
+                 Vec4[] rawCoord
+        )
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            if (rawCoord.Length != 4) {
+                throw new ArgumentOutOfRangeException(nameof(rawCoord));
+            }
+            vertexMat = rawCoord;
+        }
+
+        public Polygon4(
+               Point topLeft,
+               Point topRight,
+               Point bottomRight,
+               Point bottomLeft
+        )
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            this[0] = topLeft;
+            this[1] = topRight;
+            this[2] = bottomRight;
+            this[3] = bottomLeft;
+        }
+
+        public Polygon4(
+               CoordinationBase coordinate,
+               PlainCoordinate coordination,
+               Point topLeft,
+               Point topRight,
+               Point bottomRight,
+               Point bottomLeft
+        )
+            : base(coordinate)
+        {
+            this[0] = topLeft;
+            this[1] = topRight;
+            this[2] = bottomRight;
+            this[3] = bottomLeft;
+        }
+
+        public Polygon4(
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y
+        )
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            vertexMat = new Vec4[4];
+
+            this[0] = new Point(v0x, v0y, 0);
+            this[1] = new Point(v1x, v1y, 0);
+            this[2] = new Point(v2x, v2y, 0);
+            this[3] = new Point(v3x, v3y, 0);
+        }
+
+        public Polygon4(
+               PlainCoordinate coordination,
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y
+        )
+            : base(coordination)
+        {
+            vertexMat = new Vec4[4];
+
+            this[0] = new Point(v0x, v0y, 0);
+            this[1] = new Point(v1x, v1y, 0);
+            this[2] = new Point(v2x, v2y, 0);
+            this[3] = new Point(v3x, v3y, 0);
+        }
+
+        public override ReadOnlySpan<int> VertexIndices => Indices;
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Polygon5.cs

@@ -0,0 +1,85 @@
+// -----------------------------------------------------------------------------
+// file="Polygon5"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Numerics;
+
+namespace Nitload.Geometries
+{
+    public unsafe class Polygon5 : PlainGeometry
+    {
+
+        internal static int[] Indices = [
+            0,
+            1,
+            2,
+            0,
+            2,
+            3,
+            0,
+            3,
+            4
+        ];
+
+        public unsafe Polygon5()
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            vertexMat = new Vec4[5];
+        }
+
+        public Polygon5(
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y,
+               float v4x,
+               float v4y
+        )
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            vertexMat = new Vec4[5];
+            this[0] = new Point(v0x, v0y, 0);
+            this[1] = new Point(v1x, v1y, 0);
+            this[2] = new Point(v2x, v2y, 0);
+            this[3] = new Point(v3x, v3y, 0);
+            this[4] = new Point(v4x, v4y, 0);
+        }
+
+        public Polygon5(
+               PlainCoordinate coordination,
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y,
+               float v4x,
+               float v4y
+        )
+            : base(coordination)
+        {
+            vertexMat = new Vec4[5];
+            this[0] = new Point(v0x, v0y, 0);
+            this[1] = new Point(v1x, v1y, 0);
+            this[2] = new Point(v2x, v2y, 0);
+            this[3] = new Point(v3x, v3y, 0);
+            this[4] = new Point(v4x, v4y, 0);
+        }
+
+        public override ReadOnlySpan<int> VertexIndices => Indices;
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Polygon6.cs

@@ -0,0 +1,65 @@
+// -----------------------------------------------------------------------------
+// file="Polygon6"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+using Nitload.Mathematics.Numeric.Graphics;
+using System;
+using System.Numerics;
+
+namespace Nitload.Geometries
+{
+    public class Polygon6 : PlainGeometry
+    {
+
+        internal static int[] Indices = [
+            0,
+            1,
+            2,
+            0,
+            2,
+            3,
+            0,
+            3,
+            4,
+            0,
+            4,
+            5
+        ];
+
+        public Polygon6()
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            vertexMat = new Vec4[6];
+        }
+
+        public Polygon6(
+               PlainCoordinate coordination,
+               Vec2 vert0,
+               Vec2 vert1,
+               Vec2 vert2,
+               Vec2 vert3,
+               Vec2 vert4,
+               Vec2 vert5,
+               Vec2 vert6
+        )
+            : base(coordination)
+        {
+            vertexMat = new Vec4[6];
+            vertexMat[0].reg = vert0.reg;
+            vertexMat[1].reg = vert1.reg;
+            vertexMat[2].reg = vert2.reg;
+            vertexMat[3].reg = vert3.reg;
+            vertexMat[4].reg = vert4.reg;
+            vertexMat[5].reg = vert5.reg;
+            vertexMat[6].reg = vert6.reg;
+        }
+
+        public override ReadOnlySpan<int> VertexIndices => Indices;
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Rectangle.cs

@@ -0,0 +1,79 @@
+// -----------------------------------------------------------------------------
+// file="Rectangle"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+using Nitload.Mathematics.Numeric.Graphics;
+using System;
+using System.Numerics;
+
+namespace Nitload.Geometries
+{
+    public class Rectangle : Parallelogram
+    {
+
+        public Rectangle(
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y
+        )
+            : base(v0x, v0y, v1x, v1y, v2x, v2y, v3x, v3y)
+        {
+            #if DEBUG
+            Vec2 side0 = (this[0].Coord - this[1].Coord).XY;
+            Vec2 side1 = (this[2].Coord - this[1].Coord).XY;
+            if (side0.Cross(side1) == 0) {
+                throw new ArgumentException();
+            }
+            #endif
+        }
+
+        public Rectangle(
+               PlainCoordinate coordination,
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y
+        )
+            : base(coordination, v0x, v0y, v1x, v1y, v2x, v2y, v3x, v3y)
+        {
+            #if DEBUG
+            Vec2 side0 = (this[0].Coord - this[1].Coord).XY;
+            Vec2 side1 = (this[2].Coord - this[1].Coord).XY;
+            if (side0.Cross(side1) == 0) {
+                throw new ArgumentException();
+            }
+            #endif
+        }
+
+        /// <summary>
+        ///   Calculate movement to center
+        /// </summary>
+        /// <param _label="tick">
+        ///   CurrentPipeline time tick
+        /// </param>
+        /// <returns>
+        ///
+        /// </returns>
+        public virtual Vec3 MovCenter(
+                            uint tick
+        )
+        {
+            return default(Vec3);
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Square.cs

@@ -0,0 +1,36 @@
+// -----------------------------------------------------------------------------
+// file="Square"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+
+using System.Numerics;
+
+namespace Nitload.Geometries
+{
+    /// <summary>
+    ///   A square
+    /// </summary>
+    public class Square : Rectangle
+    {
+
+        public Square(
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y,
+               float v3x,
+               float v3y
+        )
+            : base(v0x, v0y, v1x, v1y, v2x, v2y, v3x, v3y)
+        {
+            if (true) { }
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/PlainGeometry/Triangle.cs

@@ -0,0 +1,90 @@
+// -----------------------------------------------------------------------------
+// file="Triangle"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Numerics;
+
+namespace Nitload.Geometries
+{
+    public unsafe class Triangle : PlainGeometry
+    {
+
+        internal static int[] Indices = [
+            0,
+            1,
+            2
+        ];
+
+        public Triangle()
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            vertexMat = new Vec4[3];
+        }
+
+        public Triangle(
+               Point vertex01,
+               Point vertex02,
+               Point vertex03
+        )
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            this[0] = vertex01;
+            this[1] = vertex02;
+            this[2] = vertex03;
+        }
+
+        public Triangle(
+               PlainCoordinate coordination,
+               Point vertex01,
+               Point vertex02,
+               Point vertex03
+        )
+            : base(coordination)
+        {
+            this[0] = vertex01;
+            this[1] = vertex02;
+            this[2] = vertex03;
+        }
+
+        public Triangle(
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y
+        )
+            : base(PlainCoordinate.StandardPlainCoordination)
+        {
+            this[0] = new Point(v0x, v0y, 0);
+            this[1] = new Point(v1x, v1y, 0);
+            this[2] = new Point(v2x, v2y, 0);
+        }
+
+        public Triangle(
+               PlainCoordinate coordination,
+               float v0x,
+               float v0y,
+               float v1x,
+               float v1y,
+               float v2x,
+               float v2y
+        )
+            : base(coordination)
+        {
+            this[0] = new Point(v0x, v0y, 0);
+            this[1] = new Point(v1x, v1y, 0);
+            this[2] = new Point(v2x, v2y, 0);
+        }
+
+        public override ReadOnlySpan<int> VertexIndices => Indices;
+
+    }
+}

src/Nitload.Mathematics/Geometries/Point.cs

@@ -0,0 +1,57 @@
+// -----------------------------------------------------------------------------
+// file="Point"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Collections.Generic;
+using System.Numerics;
+using System.Runtime.InteropServices;
+using System.Text;
+
+namespace Nitload.Geometries
+{
+    [StructLayout(LayoutKind.Explicit, Size = 16)]
+    public struct Point
+    {
+
+        [FieldOffset(0)] private Vec3 coord;
+        [FieldOffset(0)] internal Vector4 coordVec;
+
+        internal Point(
+                 Vec4 vector
+        )
+        {
+            coordVec = vector.vec;
+            coordVec.W = 1;
+        }
+
+        public Point(
+               Vec3 coord
+        )
+        {
+            Coord = coord;
+            coordVec.W = 1;
+        }
+
+        public Point(
+               float x,
+               float y,
+               float z
+        )
+        {
+            coordVec = new Vector4(x, y, z, 1);
+        }
+
+        public Vec3 Coord
+        {
+            get => coord;
+            set => coord = value;
+        }
+
+    }
+}

src/Nitload.Mathematics/Geometries/SpaceGeometry/Class1.cs

@@ -0,0 +1,15 @@
+// -----------------------------------------------------------------------------
+// file="Class1"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace Nitload.Geometries.SpaceGeometry
+{
+    internal class Class1 { }
+}

src/Nitload.Mathematics/Geometries/SpaceGeometry/TriangulatedSurface.cs

@@ -0,0 +1,120 @@
+// -----------------------------------------------------------------------------
+// file="TriangulatedSurface"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Graphics;
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
+using System.Text;
+
+namespace Nitload.Geometries
+{
+    public class TriangulatedSurface
+    {
+
+        private HashSet<(int begin,int end)> _sides;
+
+        private List<List<int>> _sideIndices;
+        private List<int> _triangleIndices;
+        private List<Vec3> _vertices;
+
+        public TriangulatedSurface()
+        {
+            _triangleIndices = new List<int>();
+            _vertices = new List<Vec3>();
+            _sides = new HashSet<(int,int)>();
+            _sideIndices = new List<List<int>>();
+        }
+
+        /**/
+        public IEnumerable<List<int>> SideIndicesEnumeratble => _sideIndices;
+
+        /*
+        public Triangle GetTriangle( int index )
+        {
+            Vec3HP v1,v2,v3;
+            int t = index * 3;
+            v1 = _vertices[ t + 0 ];
+            v2 = _vertices[ t + 1 ];
+            v3 = _vertices[ t + 2 ];
+            return new Triangle( v1, v2, v3 );
+        }
+        /**/
+
+        public Span<Vec3> VertexSpan => CollectionsMarshal.AsSpan(_vertices);
+
+        public Span<int> TriangleIndicesSpan => CollectionsMarshal.AsSpan(_triangleIndices);
+
+        /**/
+        public void AddIndicesByOrder(
+                    int v1,
+                    int v2,
+                    int v3
+        )
+        {
+            _triangleIndices.Add(v1);
+            _triangleIndices.Add(v2);
+            _triangleIndices.Add(v3);
+            if (!(_sides.Remove((v1,v2)) || _sides.Remove((v2, v1)))) {
+                _ = _sides.Add((v1, v2));
+            }
+            if (!(_sides.Remove((v2, v3)) || _sides.Remove((v3, v2)))) {
+                _ = _sides.Add((v2, v3));
+            }
+            if (!(_sides.Remove((v3, v1)) || _sides.Remove((v1, v3)))) {
+                _ = _sides.Add((v3, v1));
+            }
+        }
+
+        public void AddVertexByOrder(
+                    Vec3 vertex
+        )
+        {
+            _vertices.Add(vertex);
+        }
+
+        public int TriangleCount()
+        {
+            return _triangleIndices.Count / 3;
+        }
+
+        public void TrimExcess()
+        {
+            _triangleIndices.TrimExcess();
+            _vertices.TrimExcess();
+            Dictionary<int, int> side_dictionary = _sides.ToDictionary(x => x.begin, x => x.end);
+            int begin = side_dictionary.First().Key;
+            List<int> sideLoop = new(begin);
+            _sideIndices.Add(sideLoop);
+            while (side_dictionary.Count > 0)
+            {
+                if (side_dictionary.Remove(begin, out begin))
+                {
+                    sideLoop.Add(begin);
+                } else
+                {
+                    begin = side_dictionary.First().Key;
+                    sideLoop = new List<int>(begin);
+                    _sideIndices.Add(sideLoop);
+                }
+            }
+            _sideIndices.TrimExcess();
+        }
+
+        public int VertexCount()
+        {
+            return _vertices.Count;
+        }
+
+    }
+}

src/Nitload.Mathematics/Internal/ColumnSet2.cs

@@ -0,0 +1,77 @@
+// -----------------------------------------------------------------------------
+// file="ColumnSet2"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Text;
+
+namespace Nitload.Mathematics.Internal
+{
+    [StructLayout(LayoutKind.Explicit, Size = 12 * sizeof(float), Pack = 16)]
+    internal struct ColumnSet2 : IEquatable<ColumnSet2>
+    {
+
+        [FieldOffset(0 * sizeof(float))] internal float m00;
+        [FieldOffset(4 * sizeof(float))] internal float m01;
+        [FieldOffset(1 * sizeof(float))] internal float m10;
+        [FieldOffset(5 * sizeof(float))] internal float m11;
+        [FieldOffset(2 * sizeof(float))] internal float m20;
+        [FieldOffset(6 * sizeof(float))] internal float m21;
+        [FieldOffset(3 * sizeof(float))] internal float m30;
+        [FieldOffset(7 * sizeof(float))] internal float m31;
+
+        [FieldOffset(0)] internal Vector4 colum0;
+        [FieldOffset(4 * sizeof(float))] internal Vector4 colum1;
+
+        internal ColumnSet2(
+                 Vector4 c0,
+                 Vector4 c1
+        )
+        {
+            Unsafe.SkipInit(out this);
+            colum0 = c0;
+            colum1 = c1;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public override bool Equals(
+                             object? obj
+        )
+        {
+            return (obj is ColumnSet2 mat) && Equals(mat);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public bool Equals(
+                    ColumnSet2 other
+        )
+        {
+            return colum0.Equals(other.colum0) && colum1.Equals(other.colum1);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator !=(
+                                    ColumnSet2 left,
+                                    ColumnSet2 right
+        )
+        {
+            return !(left == right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator ==(
+                                    ColumnSet2 left,
+                                    ColumnSet2 right
+        )
+        {
+            return left.Equals(right);
+        }
+
+    }
+}

src/Nitload.Mathematics/Internal/ColumnSet3.cs

@@ -0,0 +1,75 @@
+// -----------------------------------------------------------------------------
+// file="ColumnSet3"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Internal;
+
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Internal
+{
+    [StructLayout(LayoutKind.Explicit, Size = 12 * sizeof(float), Pack = 16)]
+    internal struct ColumnSet3 : IEquatable<ColumnSet3>
+    {
+
+        [FieldOffset(0 * sizeof(float))] internal float m00;
+        [FieldOffset(4 * sizeof(float))] internal float m01;
+        [FieldOffset(8 * sizeof(float))] internal float m02;
+        [FieldOffset(1 * sizeof(float))] internal float m10;
+        [FieldOffset(5 * sizeof(float))] internal float m11;
+        [FieldOffset(9 * sizeof(float))] internal float m12;
+        [FieldOffset(2 * sizeof(float))] internal float m20;
+        [FieldOffset(6 * sizeof(float))] internal float m21;
+        [FieldOffset(10 * sizeof(float))] internal float m22;
+        [FieldOffset(3 * sizeof(float))] internal float m30;
+        [FieldOffset(7 * sizeof(float))] internal float m31;
+        [FieldOffset(11 * sizeof(float))] internal float m32;
+
+        [FieldOffset(0)] internal Vector4 colum0;
+        [FieldOffset(4 * sizeof(float))] internal Vector4 colum1;
+        [FieldOffset(8 * sizeof(float))] internal Vector4 colum2;
+
+        public override bool Equals(
+                             object? obj
+        )
+        {
+            return (obj is ColumnSet3 mat) && Equals(mat);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public bool Equals(
+                    ColumnSet3 other
+        )
+        {
+            return colum0.Equals(other.colum0) &&
+                   colum1.Equals(other.colum1) &&
+                   colum2.Equals(other.colum2);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator !=(
+                                    ColumnSet3 left,
+                                    ColumnSet3 right
+        )
+        {
+            return !(left == right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator ==(
+                                    ColumnSet3 left,
+                                    ColumnSet3 right
+        )
+        {
+            return left.Equals(right);
+        }
+
+    }
+}

src/Nitload.Mathematics/Internal/ColumnSet4.cs

@@ -0,0 +1,96 @@
+// -----------------------------------------------------------------------------
+// file="ColumnSet4"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload;
+using Nitload.Mathematics.Internal;
+
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Internal
+{
+    [StructLayout(LayoutKind.Explicit, Size = 16 * sizeof(float))]
+    internal struct ColumnSet4 : IEquatable<ColumnSet4>
+    {
+
+        [FieldOffset(0)]  internal Vector4 colum0;
+        [FieldOffset(16)] internal Vector4 colum1;
+        [FieldOffset(32)] internal Vector4 colum2;
+        [FieldOffset(48)] internal Vector4 colum3;
+
+        [FieldOffset(0)] public float m00;
+        [FieldOffset(16)] public float m01;
+        [FieldOffset(32)] public float m02;
+        [FieldOffset(48)] public float m03;
+        [FieldOffset(4)] public float m10;
+        [FieldOffset(20)] public float m11;
+        [FieldOffset(36)] public float m12;
+        [FieldOffset(52)] public float m13;
+        [FieldOffset(8)] public float m20;
+        [FieldOffset(24)] public float m21;
+        [FieldOffset(40)] public float m22;
+        [FieldOffset(56)] public float m23;
+        [FieldOffset(12)] public float m30;
+        [FieldOffset(28)] public float m31;
+        [FieldOffset(44)] public float m32;
+        [FieldOffset(60)] public float m33;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal ColumnSet4(
+                 Vector4 colum0,
+                 Vector4 colum1,
+                 Vector4 colum2,
+                 Vector4 colum3
+        )
+        {
+            Unsafe.SkipInit(out this);
+            this.colum0 = colum0;
+            this.colum1 = colum1;
+            this.colum2 = colum2;
+            this.colum3 = colum3;
+        }
+
+        public override bool Equals(
+                             object? obj
+        )
+        {
+            return (obj is ColumnSet4 mat) && Equals(mat);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public bool Equals(
+                    ColumnSet4 other
+        )
+        {
+            return colum0.Equals(other.colum0) &&
+                   colum1.Equals(other.colum1) &&
+                   colum2.Equals(other.colum2) &&
+                   colum3.Equals(other.colum3);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator !=(
+                                    ColumnSet4 left,
+                                    ColumnSet4 right
+        )
+        {
+            return !(left == right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator ==(
+                                    ColumnSet4 left,
+                                    ColumnSet4 right
+        )
+        {
+            return left.Equals(right);
+        }
+
+    }
+}

src/Nitload.Mathematics/Internal/IntrinsicContext.cs

@@ -0,0 +1,183 @@
+// -----------------------------------------------------------------------------
+// file="IntrinsicContext"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric.Runtime;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Internal
+{
+    internal enum MatrixIntrinsicHandle : int
+    {
+
+        Add = 0,
+        Fill = 1,
+        ScalarMul = 2,
+        Sub = 3,
+        Transpose = 4,
+        ClearPanel = 5,
+        Fma = 6,
+        PackFmaLeft = 7,
+        PackFmaRight = 8,
+        Pivot = 9,
+
+    }
+
+    [StructLayout(LayoutKind.Explicit, Size = 4)]
+    public struct MatrixElementType
+    {
+
+        public const int IntrinsicCount = 0x0a;
+
+        public static readonly MatrixElementType F32 = new(0x00000004);
+        public static readonly MatrixElementType F64 = new(0x00010008);
+        public static readonly MatrixElementType I08 = new(0x00020001);
+        public static readonly MatrixElementType I16 = new(0x00030002);
+        public static readonly MatrixElementType I32 = new(0x00040004);
+        public static readonly MatrixElementType I64 = new(0x00050008);
+        public static readonly MatrixElementType U08 = new(0x00060001);
+        public static readonly MatrixElementType U16 = new(0x00070002);
+        public static readonly MatrixElementType U32 = new(0x00080004);
+        public static readonly MatrixElementType U64 = new(0x00090008);
+
+        [FieldOffset(0)] internal readonly int Mask;
+        [FieldOffset(2)] public readonly ushort IntrinsicNode;
+
+        [FieldOffset(0)] public readonly ushort Size;
+
+        private MatrixElementType(
+                int mask
+        )
+        {
+            Unsafe.SkipInit(out this);
+            Unsafe.SkipInit(out Size);
+            Unsafe.SkipInit(out IntrinsicNode);
+            this.Mask = mask;
+        }
+
+        [Obsolete("Do not use this constructor. Use pre-defined ones from static fields.", true)]
+        public MatrixElementType()
+        {
+            throw new NotSupportedException("Do not use this constructor. Use pre-defined ones from static fields.");
+        }
+
+        public override readonly bool Equals(
+                                      object? obj
+        )
+        {
+            return (obj is MatrixElementType type) && (Mask == type.Mask);
+        }
+
+        public override readonly int GetHashCode()
+        {
+            return Mask;
+        }
+
+    }
+
+    [StructLayout(LayoutKind.Sequential)]
+    internal unsafe struct IntrinsicContext
+    {
+
+        public delegate* unmanaged[Cdecl]<char*, int, int> city_hash_simplify;              // city_hash_simplify
+        public delegate* unmanaged[Cdecl]<int, ulong> factorial_simd;                       // factorial_simd
+        public delegate* unmanaged[Cdecl]<char*, int, int, int> index_pair;                 // index_pair
+        public delegate* unmanaged[Cdecl]<Vector4*, Vector4*, void> vec_f32_normalize_3;   // f32_normalize_3
+        public MatrixMultiKernelIntrinsicArray mat_mk;                                              // mat[4] in native order
+        public MatrixSingleKernelIntrinsicArray mat_sk;                                              // mat[4] in native order
+
+        public IntrinsicContext()
+        {
+            Unsafe.SkipInit(out this);
+        }
+
+    }
+
+    #region mk
+
+    [InlineArray(10)]
+    internal unsafe struct MatrixMultiKernelIntrinsic
+    {
+
+        private nint handle;              //f32_panel_mul
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly void Call(
+                             MatrixIntrinsicHandle index,
+                             MatrixParallelTask* task
+        )
+        {
+            nint ptr = this[(int)index];
+            ((delegate* unmanaged[Cdecl]<MatrixParallelTask*, void>)ptr)(task);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly nint Get(
+                             MatrixIntrinsicHandle index
+        )
+        {
+            return this[(int)index];
+        }
+
+    }
+
+    [InlineArray(MatrixElementType.IntrinsicCount)]
+    internal unsafe struct MatrixMultiKernelIntrinsicArray
+    {
+
+        private MatrixMultiKernelIntrinsic e0;
+
+    }
+
+        #endregion
+
+    #region sk
+
+    [InlineArray(10)]
+    internal unsafe struct MatrixSingleKernelIntrinsic
+    {
+
+        private nint handle;              //f32_panel_mul
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly void Call(
+                             MatrixIntrinsicHandle index,
+                             MatrixKernel* left,
+                             MatrixKernel* right,
+                             MatrixKernel* ans,
+                             ulong scalar
+        )
+        {
+            nint ptr = this[(int)index];
+            ((delegate* unmanaged[Cdecl]<MatrixKernel*, MatrixKernel*, MatrixKernel*, ulong, void>)ptr)(left,
+                                                                                                        right,
+                                                                                                        ans,
+                                                                                                        scalar);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly nint Get(
+                             MatrixIntrinsicHandle index
+        )
+        {
+            return this[(int)index];
+        }
+
+    }
+
+    [InlineArray(MatrixElementType.IntrinsicCount)]
+    internal unsafe struct MatrixSingleKernelIntrinsicArray
+    {
+
+        private MatrixSingleKernelIntrinsic e0;
+
+    }
+
+    #endregion
+
+}

src/Nitload.Mathematics/Internal/MatrixKernel.cs

@@ -0,0 +1,55 @@
+// -----------------------------------------------------------------------------
+// file="MatrixKernel"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics
+{
+    public unsafe struct MatrixKernel
+    {
+
+        private byte Head;
+
+    }
+
+    [StructLayout(LayoutKind.Sequential)]
+    [InlineArray(16)]
+    public unsafe struct MatrixKernel<T> where T : unmanaged, INumber<T>
+    {
+
+        public static readonly int Size = 16 * sizeof(T);
+
+        private T c0r0;
+
+        public MatrixKernel()
+        {
+            Unsafe.SkipInit(out this);
+        }
+
+        public MatrixKernel(
+               bool clear
+        )
+        {
+            if (!clear) {
+                Unsafe.SkipInit(out this);
+            }
+        }
+
+        public ref T this[
+                     int x,
+                     int y
+        ]
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get => ref Unsafe.Add(ref c0r0, (x * 4) + y);
+        }
+
+    }
+}

src/Nitload.Mathematics/Internal/MatrixPanel.cs

@@ -0,0 +1,177 @@
+// -----------------------------------------------------------------------------
+// file="MatrixPanel"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Internal;
+using Nitload.Mathematics.Numeric;
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public unsafe struct MatrixPanel
+    {
+
+        public const int AlignmentBytes = 64;
+
+        private byte Head;
+
+        // Allocate aligned memory for one panel; caller must free via Destroy.
+        public static unsafe MatrixPanel<T>* Create<T>() where T : unmanaged, INumber<T>
+        {
+            return (MatrixPanel<T>*)NativeMemory.AlignedAlloc(Size<T>(), AlignmentBytes);
+        }
+
+        public static unsafe void Destroy<T>(
+                                  MatrixPanel<T>* panel
+        ) where T : unmanaged, INumber<T>
+        {
+            NativeMemory.AlignedFree(panel);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe Span<T> Elements<T>(
+                                     ref T head
+        ) where T : unmanaged, INumber<T>
+        {
+            void* p = Unsafe.AsPointer(ref head);
+            return new Span<T>(p, MaxElementCount<T>());
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static int ElementSideLength<T>() where T : unmanaged, INumber<T>
+        {
+            return PanelSideLength<T>();
+        }
+
+        public static int MaxElementCount<T>() where T : unmanaged, INumber<T>
+        {
+            int side = ElementSideLength<T>();
+            return side * side;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static nuint Size<T>() where T : unmanaged, INumber<T>
+        {
+            return (nuint)(MaxElementCount<T>() * sizeof(T));
+        }
+
+        /// <summary>
+        ///   Number of 4x4 kernels that fit in a panel for the given ISA and element T.
+        /// </summary>
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static int KernelCapacity<T>() where T : unmanaged, INumber<T>
+        {
+            int perSide = KernelSideCount<T>();
+            return perSide * perSide;
+        }
+
+        /// <summary>
+        ///   How many 4x4 kernels fit along a single panel dimension for the given ISA.
+        /// </summary>
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static int KernelSideCount<T>() where T : unmanaged, INumber<T>
+        {
+            SIMDID id = NumericalModule.GlobalIntrinsicMask;
+            switch (id.MaskByte[0] & 0b_00001111)
+            {
+                case 1:
+                    return id.MaskByte[1] switch
+                    {
+                        1 => KernelSideCountSse<T>(),
+                        2 => KernelSideCountAvx2<T>(),
+                        3 => KernelSideCountAvx512<T>(),
+                        4 => KernelSideCountAvx512<T>(),
+                        _ => throw new NotSupportedException("Unsupported SIMD level."),
+                    };
+                default:
+                    throw new NotSupportedException("Current hardware platform is not supported");
+            }
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static unsafe nuint Align64(
+                                    void* p
+        )
+        {
+            return ((nuint)p + (AlignmentBytes - 1)) & ~(nuint)(AlignmentBytes - 1);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int KernelSideCountAvx2<T>() where T : unmanaged, INumber<T>
+        {
+            return PanelSideFromBits<T>(256) / 4;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int KernelSideCountAvx512<T>() where T : unmanaged, INumber<T>
+        {
+            return PanelSideFromBits<T>(512) / 4;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int KernelSideCountSse<T>() where T : unmanaged, INumber<T>
+        {
+            return PanelSideFromBits<T>(128) / 4;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int PanelSideFromBits<T>(
+                           int vectorBits
+        ) where T : unmanaged, INumber<T>
+        {
+            int elemBits = sizeof(T) * 8;
+
+            // Q = max(4, vectorBits/64, vectorBits/elemBits)
+            int by64 = vectorBits / 64;
+            int byElem = vectorBits / elemBits;
+            int side = Math.Max(4, Math.Max(by64, byElem));
+            return side;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int PanelSideLength<T>() where T : unmanaged, INumber<T>
+        {
+            SIMDID id = NumericalModule.GlobalIntrinsicMask;
+            switch (id.MaskByte[0] & 0b_00001111)
+            {
+                case 1:
+                    return id.MaskByte[1] switch
+                    {
+                        1 => PanelSideFromBits<T>(128),
+                        2 => PanelSideFromBits<T>(256),
+                        3 => PanelSideFromBits<T>(512),
+                        4 => PanelSideFromBits<T>(512),
+                        _ => throw new NotSupportedException("Unsupported SIMD level."),
+                    };
+                default:
+                    throw new NotSupportedException("Current hardware platform is not supported");
+            }
+        }
+
+    }
+
+    // Single POD-style panel: ZMM-compatible size. Alignment handled by allocator, not by the struct.
+    public unsafe struct MatrixPanel<T> where T : unmanaged, INumber<T>
+    {
+
+        internal T Head; // first element of the panel buffer 
+
+        public ref T this[
+                     int column,
+                     int row
+        ]
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            get => ref Unsafe.Add(ref Head, row * MatrixPanel.ElementSideLength<T>() + column);
+        }
+
+        public Span<T> Elements => MatrixPanel.Elements(ref Head);
+
+    }
+}

src/Nitload.Mathematics/Internal/SIMDID.cs

@@ -0,0 +1,54 @@
+// -----------------------------------------------------------------------------
+// file="SIMDID"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Internal
+{
+    [StructLayout(LayoutKind.Explicit)]
+    internal unsafe struct SIMDID
+    {
+
+        [FieldOffset(0)] public fixed byte MaskByte[8];
+        [FieldOffset(0)] public ulong Mask;
+
+        internal SIMDID(
+                 ulong data
+        )
+        {
+            Mask = data;
+        }
+
+        public SIMDID()
+        {
+            Unsafe.SkipInit(out this);
+        }
+
+        public static bool operator !=(
+                                    SIMDID left,
+                                    SIMDID right
+        )
+        {
+            return left.Mask != right.Mask;
+        }
+
+        public static bool operator ==(
+                                    SIMDID left,
+                                    SIMDID right
+        )
+        {
+            return left.Mask == right.Mask;
+        }
+
+    }
+}

src/Nitload.Mathematics/Internal/ScalarPacket.cs

@@ -0,0 +1,49 @@
+// -----------------------------------------------------------------------------
+// file="ScalarPacket"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.HighPerformance.ProcessorAbstraction;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics
+{
+    public unsafe struct ScalarPacket
+    {
+
+        private fixed ulong data[8];
+
+        public void BroadCast<TScalar, TRegister>(
+                    int slotIndex,
+                    TScalar value
+        ) where TScalar : unmanaged, INumber<TScalar>
+            where TRegister : unmanaged, IRegisterPresentation
+        {
+            if (slotIndex is < 0 or >= 8) {
+                throw new IndexOutOfRangeException();
+            }
+            ref TRegister register = ref Unsafe.As<ulong, TRegister>(ref data[slotIndex]);
+            register.BroadCastFrom(value);
+        }
+
+        public ref T Data<T>(
+                     int index
+        ) where T : unmanaged, INumber<T>
+        {
+            if (index is < 0 or >= 8) {
+                throw new IndexOutOfRangeException();
+            }
+
+            return ref Unsafe.As<ulong, T>(ref data[index]);
+        }
+
+    }
+}

src/Nitload.Mathematics/Internal/UnsafeCompute.cs

@@ -0,0 +1,46 @@
+// -----------------------------------------------------------------------------
+// file="UnsafeCompute"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric;
+using System.Runtime.CompilerServices;
+
+namespace Nitload.Mathematics.Internal
+{
+    public static class UnsafeCompute
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe int CityHashSimplify(
+                                 char* span,
+                                 int length
+        )
+        {
+            return NumericalModule.GlobalContext.city_hash_simplify(span, length);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe int CityHashSimplify(
+                                 char* span,
+                                 int begin,
+                                 int length
+        )
+        {
+            return NumericalModule.GlobalContext.city_hash_simplify(span + begin, length);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe int IndexOfCharPair(
+                                 char* source,
+                                 int pairValue,
+                                 int length
+        )
+        {
+            return NumericalModule.GlobalContext.index_pair(source, pairValue, length);
+        }
+
+    }
+}

src/Nitload.Mathematics/Nitload.Mathematics.csproj

@@ -0,0 +1,11 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <ItemGroup>
+    <None Include="Numeric\Runtime\MatrixParallel.Publish.cs" />
+  </ItemGroup>
+
+  <ItemGroup>
+    <ProjectReference Include="..\Nitload.Common\Nitload.Common.csproj" />
+  </ItemGroup>
+
+</Project>

src/Nitload.Mathematics/Numeric/Graphics/Cordination.cs

@@ -0,0 +1,18 @@
+// -----------------------------------------------------------------------------
+// file="Cordination"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    public enum CoordinatePlain
+    {
+
+        XY,
+        YZ,
+        XZ,
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/GraphicVector.cs

@@ -0,0 +1,42 @@
+// -----------------------------------------------------------------------------
+// file="GraphicVector"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    internal static class GraphicVector
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool IsNormalVectorElement<T>() where T : unmanaged, INumber<T>
+        {
+            return typeof(T) == typeof(float) ||
+                   typeof(T) == typeof(int) ||
+                   typeof(T) == typeof(uint);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool IsWideVectorElement<T>() where T : unmanaged, INumber<T>
+        {
+            return typeof(T) == typeof(double) ||
+                   typeof(T) == typeof(long) ||
+                   typeof(T) == typeof(ulong);
+        }
+
+        public static NotSupportedException ThrowUnsupportedTypeException<T>()
+        {
+            return new NotSupportedException($"Operation is not supported for type {typeof(T)}");
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/GraphicsMatrix.cs

@@ -0,0 +1,261 @@
+// -----------------------------------------------------------------------------
+// file="GraphicsMatrix"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Internal;
+using Nitload.Mathematics.Numeric;
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics.X86;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    [StructLayout(LayoutKind.Explicit, Size = 16 * sizeof(float), Pack = 16)]
+    public unsafe struct GraphicsMatrix : IEquatable<GraphicsMatrix>
+    {
+
+        [FieldOffset(0 * sizeof(float))] internal Column column;
+
+        [FieldOffset(0)] internal ColumnSet4 mat;
+        [FieldOffset(0)] internal MatrixKernel<float> kernel;
+
+        internal GraphicsMatrix(
+                 ColumnSet4 mat
+        )
+        {
+            Unsafe.SkipInit(out this);
+            this.mat = mat;
+        }
+
+        internal GraphicsMatrix(
+                 Vector4 vec0,
+                 Vector4 vec1,
+                 Vector4 vec2,
+                 Vector4 vec3
+        )
+        {
+            Unsafe.SkipInit(out this);
+            mat.colum0 = vec0;
+            mat.colum1 = vec1;
+            mat.colum2 = vec2;
+            mat.colum3 = vec3;
+        }
+
+        public GraphicsMatrix()
+        {
+            Unsafe.SkipInit(out this);
+        }
+
+        public unsafe GraphicsMatrix(
+                      float a00,
+                      float a01,
+                      float a02,
+                      float a03,
+                      float a10,
+                      float a11,
+                      float a12,
+                      float a13,
+                      float a20,
+                      float a21,
+                      float a22,
+                      float a23,
+                      float a30,
+                      float a31,
+                      float a32,
+                      float a33
+        )
+        {
+            Unsafe.SkipInit(out this);
+            mat.colum0 = new Vector4(a00, a10, a20, a30);
+            mat.colum1 = new Vector4(a01, a11, a21, a31);
+            mat.colum2 = new Vector4(a02, a12, a22, a32);
+            mat.colum3 = new Vector4(a03, a13, a23, a33);
+        }
+
+        public unsafe float this[
+                            int rowNum,
+                            int columNum
+        ]
+        {
+            get
+            {
+#if DEBUG
+                if (rowNum is > 3 or < 0)
+                {
+                    throw new ArgumentOutOfRangeException(nameof(columNum));
+                }
+                if (columNum is > 3 or < 0) {
+                    throw new ArgumentOutOfRangeException(nameof(columNum));
+                }
+#endif
+                return kernel[rowNum + (4 * columNum)];
+            }
+            set
+            {
+#if DEBUG
+                if (rowNum is > 3 or < 0)
+                {
+                    throw new ArgumentOutOfRangeException(nameof(columNum));
+                }
+                if (columNum is > 3 or < 0) {
+                    throw new ArgumentOutOfRangeException(nameof(columNum));
+                }
+#endif
+                kernel[rowNum + (4 * columNum)] = value;
+            }
+        }
+
+        public GraphicsMatrix Transpose
+        {
+            get
+            {
+                GraphicsMatrix transpose = new();
+                InternalTranspose(ref transpose);
+                return transpose;
+            }
+        }
+
+        public static GraphicsMatrix Unit => new GraphicsMatrix(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0,
+                                                                0, 0, 0, 1);
+
+        public static GraphicsMatrix WNegativeUnit => new GraphicsMatrix(1, 0, 0, 0, 0, 1, 0, 0, 0,
+                                                                         0, 1, 0, 0, 0, 0, -1);
+
+        public static GraphicsMatrix XNegativeUnit => new GraphicsMatrix(-1, 0, 0, 0, 0, 1, 0, 0, 0,
+                                                                         0, 1, 0, 0, 0, 0, 1);
+
+        public static GraphicsMatrix YNegativeUnit => new GraphicsMatrix(1, 0, 0, 0, 0, -1, 0, 0, 0,
+                                                                         0, 1, 0, 0, 0, 0, 1);
+
+        public static GraphicsMatrix ZNegativeUnit => new GraphicsMatrix(1, 0, 0, 0, 0, 1, 0, 0, 0,
+                                                                         0, -1, 0, 0, 0, 0, 1);
+
+        public ReadOnlySpan<float> AsSpan()
+        {
+            return MemoryMarshal.CreateReadOnlySpan(ref Unsafe.As<ColumnSet4, float>(ref mat), 16);
+        }
+
+        public override bool Equals(
+                             object? obj
+        )
+        {
+            return (obj is GraphicsMatrix x) && Equals(x);
+        }
+
+        public bool Equals(
+                    GraphicsMatrix other
+        )
+        {
+            return mat.Equals(other.mat);
+        }
+
+        public override int GetHashCode()
+        {
+            return HashCode.Combine(mat.colum0.GetHashCode(), mat.colum1.GetHashCode(),
+                                    mat.colum2.GetHashCode(), mat.colum3.GetHashCode());
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal unsafe void InternalTranspose(
+                             ref GraphicsMatrix transpose
+        )
+        {
+            const int f32 = 0;
+            fixed (MatrixKernel<float>* source = &this.kernel, target = &transpose.kernel) {
+                NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Transpose,
+                                                               null, (MatrixKernel*)source,
+                                                               (MatrixKernel*)target, 0);
+            }
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe GraphicsMatrix operator -(
+                                                     GraphicsMatrix left,
+                                                     GraphicsMatrix right
+        )
+        {
+            GraphicsMatrix ret = new();
+            int f32 = MatrixElementType.F32.IntrinsicNode;
+            NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Sub,
+                                                           (MatrixKernel*)&left.kernel,
+                                                           (MatrixKernel*)&right.kernel,
+                                                           (MatrixKernel*)&ret.kernel, 1);
+            return ret;
+        }
+
+        public static bool operator !=(
+                                    GraphicsMatrix left,
+                                    GraphicsMatrix right
+        )
+        {
+            return !(left == right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe GraphicsMatrix operator *(
+                                                     GraphicsMatrix mat,
+                                                     float value
+        )
+        {
+            int f32 = MatrixElementType.F32.IntrinsicNode;
+            return new GraphicsMatrix(mat.mat.colum0 * value, mat.mat.colum1 * value,
+                                      mat.mat.colum2 * value, mat.mat.colum3 * value);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static GraphicsMatrix operator *(
+                                              GraphicsMatrix left,
+                                              GraphicsMatrix right
+        )
+        {
+            GraphicsMatrix ret = new();
+            int f32 = MatrixElementType.F32.IntrinsicNode;
+            NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Fma,
+                                                           (MatrixKernel*)&left,
+                                                           (MatrixKernel*)&right,
+                                                           (MatrixKernel*)&ret, 0);
+            return ret;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static unsafe GraphicsMatrix operator +(
+                                                     GraphicsMatrix left,
+                                                     GraphicsMatrix right
+        )
+        {
+            GraphicsMatrix ret = new();
+            int f32 = MatrixElementType.F32.IntrinsicNode;
+            NumericalModule.GlobalContext.mat_sk[f32].Call(MatrixIntrinsicHandle.Add,
+                                                           (MatrixKernel*)&left,
+                                                           (MatrixKernel*)&right,
+                                                           (MatrixKernel*)&ret, 1);
+            return ret;
+        }
+
+        public static bool operator ==(
+                                    GraphicsMatrix left,
+                                    GraphicsMatrix right
+        )
+        {
+            return left.Equals(right);
+        }
+
+        [InlineArray(4)]
+        internal struct Column
+        {
+
+            private Vec4<float> value;
+
+        }
+
+    }
+}
+

src/Nitload.Mathematics/Numeric/Graphics/Vec2.cs

@@ -0,0 +1,70 @@
+// -----------------------------------------------------------------------------
+// file="Vec2"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.Intrinsics;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    public static class Vec2
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector128<T> AsVector128Safe<T>(this Vec2<T> source) where T: unmanaged, INumber<T>
+        {
+            Vector128<T> result = Vector128<T>.Zero;
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector128<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector128<T> AsVector128Unsafe<T>(this Vec2<T> source) where T: unmanaged, INumber<T>
+        {
+            Unsafe.SkipInit(out Vector128<T> result);
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector128<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vec2<T> FromVector128Safe<T>(Vector128<T> source) where T: unmanaged, INumber<T>
+        {
+            return Unsafe.As<Vector128<T>, Vec2<T>>(ref source);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vec2<T> FromVector128Unsafe<T>(Vector128<T> source) where T: unmanaged, INumber<T>
+        {
+            return Unsafe.As<Vector128<T>, Vec2<T>>(ref source);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<float> ConvertToSingle(Vec2<int> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToSingle(source.AsVector128Safe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<float> ConvertToSingle(Vec2<uint> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToSingle(source.AsVector128Safe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<double> ConvertToDouble(Vec2<long> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToDouble(source.AsVector128Safe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<double> ConvertToDouble(Vec2<ulong> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToDouble(source.AsVector128Safe()));
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vec2_T.cs

@@ -0,0 +1,129 @@
+// -----------------------------------------------------------------------------
+// file="Vec2_T"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    /// <summary>
+    ///   A two dimension vector defined by two elements. Vec2 in World are default used as colum
+    ///   vector.
+    /// </summary>
+    /// <typeparam name="T">
+    ///
+    /// </typeparam>
+    [StructLayout(LayoutKind.Sequential, Pack = 16)]
+    public struct Vec2<T>(T x, T y) : IUnmanagedVector<Vec2<T>> where T: unmanaged, INumber<T>
+    {
+
+        public T X = x;
+        public T Y = y;
+
+        public static Vec2<T> Unit => new(T.One, T.One);
+
+        public static Vec2<T> Zero => new(T.Zero, T.Zero);
+
+        public static Vec2<T> One => new(T.One, T.One);
+
+        public Vec2<T> XY
+        {
+            readonly get => this;
+            set => this = value;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public T Cross(Vec2<T> right)
+        {
+            return (X * right.Y) - (Y * right.X);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly T Dot(Vec2<T> vec)
+        {
+            Vector128<T> leftVec = this.AsVector128Unsafe();
+            Vector128<T> rightVec = vec.AsVector128Unsafe();
+            Vector128<T> result_v128 = leftVec * rightVec;
+            Vec2<T> result = Unsafe.As<Vector128<T>, Vec2<T>>(ref result_v128);
+            return result.X + result.Y;
+        }
+
+        public override readonly bool Equals(object? obj)
+        {
+            return (obj is Vec2<T> other) && (other == this);
+        }
+
+        public override readonly int GetHashCode()
+        {
+            return HashCode.Combine(X, Y);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<T> operator -(Vec2<T> left, Vec2<T> right)
+        {
+            Vector128<T> leftVec = left.AsVector128Unsafe();
+            Vector128<T> rightVec = right.AsVector128Unsafe();
+            return Vec2.FromVector128Unsafe(leftVec - rightVec);
+        }
+
+        public static bool operator !=(Vec2<T> left, Vec2<T> right)
+        {
+            return !(left == right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<T> operator *(Vec2<T> vec, T value)
+        {
+            Vector128<T> leftVec = vec.AsVector128Unsafe();
+            return Vec2.FromVector128Unsafe(leftVec * value);
+        }
+
+        public static Vec2<T> operator *(T value, Vec2<T> vec)
+        {
+            Vector128<T> leftVec = vec.AsVector128Unsafe();
+            return Vec2.FromVector128Unsafe(leftVec * value);
+        }
+
+        public static Vec2<T> operator /(Vec2<T> vec, T value)
+        {
+            Vector128<T> leftVec = vec.AsVector128Unsafe();
+            return Vec2.FromVector128Unsafe(leftVec / value);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<T> operator +(Vec2<T> left, Vec2<T> right)
+        {
+            Vector128<T> leftVec = left.AsVector128Unsafe();
+            Vector128<T> rightVec = right.AsVector128Unsafe();
+            return Vec2.FromVector128Unsafe(leftVec + rightVec);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator ==(Vec2<T> left, Vec2<T> right)
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return Unsafe.As<Vec2<T>, ulong>(ref left) == Unsafe.As<Vec2<T>, ulong>(ref right);
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return Unsafe.As<Vec2<T>, Vector128<T>>(ref left) == Unsafe.As<Vec2<T>, Vector128<T>>(ref right);
+            }
+#pragma warning disable CA1065
+            throw new NotSupportedException($"Element type {typeof(T).Name} is not supported.");
+#pragma warning restore CA1065
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static implicit operator Vec2<T>((T, T) tuple)
+        {
+            return new Vec2<T>(tuple.Item1, tuple.Item2);
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vec3.cs

@@ -0,0 +1,214 @@
+// -----------------------------------------------------------------------------
+// file="Vec3"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.Intrinsics;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    public static class Vec3
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<double> ConvertToDouble(Vec3<long> source)
+        {
+            return FromVector256Unsafe(Vector256.ConvertToDouble(source.AsVector256Safe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<double> ConvertToDouble(Vec3<ulong> source)
+        {
+            return FromVector256Unsafe(Vector256.ConvertToDouble(source.AsVector256Safe()));
+        }
+
+        #region cast
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<float> ConvertToSingle(Vec3<int> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToSingle(source.AsVector128Safe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<float> ConvertToSingle(Vec3<uint> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToSingle(source.AsVector128Safe()));
+        }
+
+        public static bool IsParallel<T>(Vec3<T> left, Vec3<T> right) where T: unmanaged,INumber<T>
+        {
+            return Cross(left, right).Equals(T.Zero);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector128<T> AsVector128Safe<T>(this Vec3<T> source) where T: unmanaged, INumber<T>
+        {
+            Vector128<T> result = Vector128<T>.Zero;
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector128<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector128<T> AsVector128Unsafe<T>(this Vec3<T> source) where T: unmanaged, INumber<T>
+        {
+            Unsafe.SkipInit(out Vector128<T> result);
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector128<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector256<T> AsVector256Safe<T>(this Vec3<T> source) where T: unmanaged, INumber<T>
+        {
+            Vector256<T> result = Vector256<T>.Zero;
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector256<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector256<T> AsVector256Unsafe<T>(this Vec3<T> source) where T: unmanaged, INumber<T>
+        {
+            Unsafe.SkipInit(out Vector256<T> result);
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector256<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vec3<T> FromVector128Unsafe<T>(Vector128<T> source) where T: unmanaged, INumber<T>
+        {
+            return Unsafe.As<Vector128<T>, Vec3<T>>(ref source);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vec3<T> FromVector256Unsafe<T>(Vector256<T> source) where T: unmanaged, INumber<T>
+        {
+            return Unsafe.As<Vector256<T>, Vec3<T>>(ref source);
+        }
+
+#ednregion
+
+        #region cross
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> Cross<T>(Vec3<T> left, Vec3<T> right) where T: unmanaged, INumber<T>
+        {
+            if (typeof(T) == typeof(float))
+            {
+                Vec3<float> result = CrossSingle(
+                    Unsafe.As<Vec3<T>, Vec3<float>>(ref left),
+                    Unsafe.As<Vec3<T>, Vec3<float>>(ref right));
+                return Unsafe.As<Vec3<float>, Vec3<T>>(ref result);
+            }
+            if (typeof(T) == typeof(double))
+            {
+                Vec3<double> result = CrossDouble(
+                    Unsafe.As<Vec3<T>, Vec3<double>>(ref left),
+                    Unsafe.As<Vec3<T>, Vec3<double>>(ref right));
+                return Unsafe.As<Vec3<double>, Vec3<T>>(ref result);
+            }
+            if (typeof(T) == typeof(int))
+            {
+                Vec3<int> result = CrossInt32(
+                    Unsafe.As<Vec3<T>, Vec3<int>>(ref left),
+                    Unsafe.As<Vec3<T>, Vec3<int>>(ref right));
+                return Unsafe.As<Vec3<int>, Vec3<T>>(ref result);
+            }
+            if (typeof(T) == typeof(uint))
+            {
+                Vec3<uint> result = CrossUInt32(
+                    Unsafe.As<Vec3<T>, Vec3<uint>>(ref left),
+                    Unsafe.As<Vec3<T>, Vec3<uint>>(ref right));
+                return Unsafe.As<Vec3<uint>, Vec3<T>>(ref result);
+            }
+            if (typeof(T) == typeof(ulong))
+            {
+                Vec3<ulong> result = CrossUInt64(
+                    Unsafe.As<Vec3<T>, Vec3<ulong>>(ref left),
+                    Unsafe.As<Vec3<T>, Vec3<ulong>>(ref right));
+                return Unsafe.As<Vec3<ulong>, Vec3<T>>(ref result);
+            }
+            if (typeof(T) == typeof(long))
+            {
+                Vec3<long> result = CrossInt64(
+                    Unsafe.As<Vec3<T>, Vec3<long>>(ref left),
+                    Unsafe.As<Vec3<T>, Vec3<long>>(ref right));
+                return Unsafe.As<Vec3<long>, Vec3<T>>(ref result);
+            }
+
+            throw new NotSupportedException($"Cross product is not supported for {typeof(T)}.");
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static Vec3<double> CrossDouble(Vec3<double> left, Vec3<double> right)
+        {
+            Vector256<double> l = left.AsVector256Safe();
+            Vector256<double> r = right.AsVector256Safe();
+            Vector256<long> yzx = Vector256.Create(1L, 2L, 0L, 3L);
+            Vector256<long> zxy = Vector256.Create(2L, 0L, 1L, 3L);
+            Vector256<double> result = (Vector256.Shuffle(l, yzx) * Vector256.Shuffle(r, zxy)) - (Vector256.Shuffle(l, zxy) * Vector256.Shuffle(r, yzx));
+            return FromVector256Unsafe(result);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static Vec3<int> CrossInt32(Vec3<int> left, Vec3<int> right)
+        {
+            Vector128<int> l = left.AsVector128Safe();
+            Vector128<int> r = right.AsVector128Safe();
+            Vector128<int> yzx = Vector128.Create(1, 2, 0, 3);
+            Vector128<int> zxy = Vector128.Create(2, 0, 1, 3);
+            Vector128<int> result = (Vector128.Shuffle(l, yzx) * Vector128.Shuffle(r, zxy)) - (Vector128.Shuffle(l, zxy) * Vector128.Shuffle(r, yzx));
+            return FromVector128Unsafe(result);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static Vec3<float> CrossSingle(Vec3<float> left, Vec3<float> right)
+        {
+            Vector128<float> l = left.AsVector128Safe();
+            Vector128<float> r = right.AsVector128Safe();
+            Vector128<int> yzx = Vector128.Create(1, 2, 0, 3);
+            Vector128<int> zxy = Vector128.Create(2, 0, 1, 3);
+            Vector128<float> result = (Vector128.Shuffle(l, yzx) * Vector128.Shuffle(r, zxy)) - (Vector128.Shuffle(l, zxy) * Vector128.Shuffle(r, yzx));
+            return FromVector128Unsafe(result);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static Vec3<uint> CrossUInt32(Vec3<uint> left, Vec3<uint> right)
+        {
+            Vector128<uint> l = left.AsVector128Safe();
+            Vector128<uint> r = right.AsVector128Safe();
+            Vector128<uint> yzx = Vector128.Create(1U, 2U, 0U, 3U);
+            Vector128<uint> zxy = Vector128.Create(2U, 0U, 1U, 3U);
+            Vector128<uint> result = (Vector128.Shuffle(l, yzx) * Vector128.Shuffle(r, zxy)) - (Vector128.Shuffle(l, zxy) * Vector128.Shuffle(r, yzx));
+            return FromVector128Unsafe(result);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static Vec3<ulong> CrossUInt64(Vec3<ulong> left, Vec3<ulong> right)
+        {
+            Vector256<ulong> l = left.AsVector256Safe();
+            Vector256<ulong> r = right.AsVector256Safe();
+            Vector256<ulong> yzx = Vector256.Create(1UL, 2UL, 0UL, 3UL);
+            Vector256<ulong> zxy = Vector256.Create(2UL, 0UL, 1UL, 3UL);
+            Vector256<ulong> result = (Vector256.Shuffle(l, yzx) * Vector256.Shuffle(r, zxy)) - (Vector256.Shuffle(l, zxy) * Vector256.Shuffle(r, yzx));
+            return FromVector256Unsafe(result);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static Vec3<long> CrossInt64(Vec3<long> left, Vec3<long> right)
+        {
+            Vector256<long> l = left.AsVector256Safe();
+            Vector256<long> r = right.AsVector256Safe();
+            Vector256<long> yzx = Vector256.Create(1U, 2U, 0U, 3U);
+            Vector256<long> zxy = Vector256.Create(2U, 0U, 1U, 3U);
+            Vector256<long> result = (Vector256.Shuffle(l, yzx) * Vector256.Shuffle(r, zxy)) - (Vector256.Shuffle(l, zxy) * Vector256.Shuffle(r, yzx));
+            return FromVector256Unsafe(result);
+        }
+
+        #endregion
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vec3_T.cs

@@ -0,0 +1,260 @@
+// -----------------------------------------------------------------------------
+// file="Vec3_T"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric;
+using System;
+using System.Collections;
+using System.Collections.Generic;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    /// <summary>
+    ///   A tightly packed three-dimensional numeric vector.
+    /// </summary>
+    [StructLayout(LayoutKind.Sequential, Pack = 16)]
+    public struct Vec3<T> : IUnmanagedVector<Vec3<T>>, IEquatable<Vec3<T>>, IEnumerable<T>
+        where T: unmanaged, INumber<T>
+    {
+
+        public T X;
+        public T Y;
+        public T Z;
+
+        public Vec3(Vec2<T> xy, T z)
+        {
+            X = xy.X;
+            Y = xy.Y;
+            Z = z;
+        }
+
+        public Vec3(T x, T y, T z)
+        {
+            X = x;
+            Y = y;
+            Z = z;
+        }
+
+        public Vec2<T> XY
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            readonly get => new(X, Y);
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            set
+            {
+                X = value.X;
+                Y = value.Y;
+            }
+        }
+
+        public Vec3<T> XYZ
+        {
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            readonly get => this;
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            set => this = value;
+        }
+
+        public static Vec3<T> Zero => new(T.Zero, T.Zero, T.Zero);
+
+        public static Vec3<T> One => new(T.One, T.One, T.One);
+
+        public static Vec3<T> UnitX => new(T.One, T.Zero, T.Zero);
+
+        public static Vec3<T> UnitY => new(T.Zero, T.One, T.Zero);
+
+        public static Vec3<T> UnitZ => new(T.Zero, T.Zero, T.One);
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly T Dot(Vec3<T> right)
+        {
+            if (Unsafe.SizeOf<T>() == 8)
+            {
+                Vec4<T> v = Vec4.FromVector256Unsafe(this.AsVector256Unsafe() * right.AsVector256Unsafe());
+                return v.X + v.Y + v.Z;
+            }
+            if (Unsafe.SizeOf<T>() == 4)
+            {
+                Vec4<T> v = Vec4.FromVector128Unsafe(this.AsVector128Unsafe() * right.AsVector128Unsafe());
+                return v.X + v.Y + v.Z;
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly bool Equals(Vec3<T> other)
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return this.AsVector128Safe() == other.AsVector128Safe();
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return this.AsVector256Safe() == other.AsVector256Safe();
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        public override readonly bool Equals(object? obj)
+        {
+            return obj is Vec3<T> other && Equals(other);
+        }
+
+        public static Vec3<T> FromSpan(ReadOnlySpan<T> span)
+        {
+            if (span.Length < 3) {
+                throw new ArgumentException("Span length must be at least 3.", nameof(span));
+            }
+
+            return new(span[0], span[1], span[2]);
+        }
+
+        public readonly IEnumerator<T> GetEnumerator()
+        {
+            yield return X;
+            yield return Y;
+            yield return Z;
+        }
+
+        public override readonly int GetHashCode()
+        {
+            return HashCode.Combine(X, Y, Z);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly T GetMaxValue()
+        {
+            return T.Max(X, T.Max(Y, Z));
+        }
+
+        public static bool IsParallel(Vec3<T> left, Vec3<T> right)
+        {
+            return Vec3.Cross(left, right) == Zero;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> operator -(Vec3<T> value)
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return Vec3.FromVector128Unsafe(Vector128<T>.Zero - value.AsVector128Unsafe());
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return Vec3.FromVector256Unsafe(Vector256<T>.Zero - value.AsVector256Unsafe());
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> operator -(Vec3<T> left, Vec3<T> right)
+        {
+            if (Unsafe.SizeOf<T>() == 4)
+            {
+                Vector128<T> left_v128 = left.AsVector128Unsafe();
+                Vector128<T> right_v128 = right.AsVector128Unsafe();
+                return Vec3.FromVector128Unsafe(left_v128 - right_v128);
+            }
+            if (Unsafe.SizeOf<T>() == 8)
+            {
+                Vector256<T> left_v256 = left.AsVector256Unsafe();
+                Vector256<T> right_v256 = right.AsVector256Unsafe();
+                return Vec3.FromVector256Unsafe(left_v256 - right_v256);
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator !=(Vec3<T> left, Vec3<T> right)
+        {
+            return !left.Equals(right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> operator *(Vec3<T> vector, T scalar)
+        {
+            if (Unsafe.SizeOf<T>() == 4)
+            {
+                Vector128<T> vector_v128 = vector.AsVector128Unsafe();
+                Vector128<T> scalar_v128 = Vector128.Create(scalar);
+                return Vec3.FromVector128Unsafe(vector_v128 * scalar_v128);
+            }
+            if (Unsafe.SizeOf<T>() == 8)
+            {
+                Vector256<T> vector_v256 = vector.AsVector256Unsafe();
+                Vector256<T> scalar_v256 = Vector256.Create(scalar);
+                return Vec3.FromVector256Unsafe(vector_v256 * scalar_v256);
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> operator *(T scalar, Vec3<T> vector)
+        {
+            return vector * scalar;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static T operator *(Vec3<T> left, Vec3<T> right)
+        {
+            return left.Dot(right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> operator /(Vec3<T> vector, T scalar)
+        {
+            if (Unsafe.SizeOf<T>() == 4)
+            {
+                Vector128<T> vector_v128 = vector.AsVector128Unsafe();
+                Vector128<T> scalar_v128 = Vector128.Create(scalar);
+                return Vec3.FromVector128Unsafe(vector_v128 / scalar_v128);
+            }
+            if (Unsafe.SizeOf<T>() == 8)
+            {
+                Vector256<T> vector_v256 = vector.AsVector256Unsafe();
+                Vector256<T> scalar_v256 = Vector256.Create(scalar);
+                return Vec3.FromVector256Unsafe(vector_v256 / scalar_v256);
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> operator +(Vec3<T> left, Vec3<T> right)
+        {
+            if (Unsafe.SizeOf<T>() == 4)
+            {
+                Vector128<T> left_v128 = left.AsVector128Unsafe();
+                Vector128<T> right_v128 = right.AsVector128Unsafe();
+                return Vec3.FromVector128Unsafe(left_v128 + right_v128);
+            }
+            if (Unsafe.SizeOf<T>() == 8)
+            {
+                Vector256<T> left_v256 = left.AsVector256Unsafe();
+                Vector256<T> right_v256 = right.AsVector256Unsafe();
+                return Vec3.FromVector256Unsafe(left_v256 + right_v256);
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static bool operator ==(Vec3<T> left, Vec3<T> right)
+        {
+            return left.Equals(right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static implicit operator Vec3<T>((T X, T Y, T Z) tuple)
+        {
+            return new(tuple.X, tuple.Y, tuple.Z);
+        }
+
+        readonly IEnumerator IEnumerable.GetEnumerator()
+        {
+            return GetEnumerator();
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vec4.cs

@@ -0,0 +1,70 @@
+// -----------------------------------------------------------------------------
+// file="Vec4"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.Intrinsics;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    public static class Vec4
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector128<T> AsVector128Unsafe<T>(this Vec4<T> source) where T: unmanaged, INumber<T>
+        {
+            Unsafe.SkipInit(out Vector128<T> result);
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector128<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vector256<T> AsVector256Unsafe<T>(this Vec4<T> source) where T: unmanaged, INumber<T>
+        {
+            Unsafe.SkipInit(out Vector256<T> result);
+            Unsafe.WriteUnaligned(ref Unsafe.As<Vector256<T>, byte>(ref result), source);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vec4<T> FromVector128Unsafe<T>(Vector128<T> source) where T: unmanaged, INumber<T>
+        {
+            return Unsafe.As<Vector128<T>, Vec4<T>>(ref source);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static Vec4<T> FromVector256Unsafe<T>(Vector256<T> source) where T: unmanaged, INumber<T>
+        {
+            return Unsafe.As<Vector256<T>, Vec4<T>>(ref source);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<float> ConvertToSingle(Vec4<int> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToSingle(source.AsVector128Unsafe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<float> ConvertToSingle(Vec4<uint> source)
+        {
+            return FromVector128Unsafe(Vector128.ConvertToSingle(source.AsVector128Unsafe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<double> ConvertToDouble(Vec4<long> source)
+        {
+            return FromVector256Unsafe(Vector256.ConvertToDouble(source.AsVector256Unsafe()));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<double> ConvertToDouble(Vec4<ulong> source)
+        {
+            return FromVector256Unsafe(Vector256.ConvertToDouble(source.AsVector256Unsafe()));
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vec4_T.cs

@@ -0,0 +1,245 @@
+// -----------------------------------------------------------------------------
+// file="Vec4_T"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    /// <summary>
+    ///   A four-dimensional numeric vector.
+    /// </summary>
+    [StructLayout(LayoutKind.Sequential, Pack = 16)]
+    public struct Vec4<T> : IUnmanagedVector<Vec4<T>> where T : unmanaged, INumber<T>
+    {
+
+        public T W;
+
+        public T X;
+        public T Y;
+        public T Z;
+
+        public Vec4(
+               Vec3<T> v3,
+               T w
+        )
+        {
+            X = v3.X;
+            Y = v3.Y;
+            Z = v3.Z;
+            W = w;
+        }
+
+        public Vec4(
+               T x,
+               T y,
+               T z,
+               T w
+        )
+        {
+            X = x;
+            Y = y;
+            Z = z;
+            W = w;
+        }
+
+        public Vec2<T> XY
+        {
+            get => Unsafe.As<T, Vec2<T>>(ref X);
+            set
+            {
+                X = value.X;
+                Y = value.Y;
+            }
+        }
+
+        public Vec3<T> XYZ
+        {
+            get => Unsafe.As<T, Vec3<T>>(ref X);
+            set
+            {
+                T w = W;
+                Unsafe.As<T, Vec3<T>>(ref X) = value;
+                W = w;
+            }
+        }
+
+        public static Vec4<T> One => new(T.One, T.One, T.One, T.One);
+
+        public static Vec4<T> Zero => new(T.Zero, T.Zero, T.Zero, T.Zero);
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> Add(
+                              Vec4<T> left,
+                              Vec4<T> right
+        )
+        {
+            return left + right;
+        }
+
+        public ReadOnlySpan<T> AsSpan()
+        {
+            return MemoryMarshal.CreateReadOnlySpan(ref X, 4);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public readonly T Dot(
+                          Vec4<T> right
+        )
+        {
+            if (Unsafe.SizeOf<T>() == 4)
+            {
+                Vec4<T> product = Vec4.FromVector128Unsafe(this.AsVector128Unsafe() * right.AsVector128Unsafe());
+                return product.X + product.Y + product.Z + product.W;
+            }
+            if (Unsafe.SizeOf<T>() == 8)
+            {
+                Vec4<T> product = Vec4.FromVector256Unsafe(this.AsVector256Unsafe() * right.AsVector256Unsafe());
+                return product.X + product.Y + product.Z + product.W;
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        public override bool Equals(
+                             object? obj
+        )
+        {
+            return (obj is Vec4<T> v) && (v == this);
+        }
+
+        public static Vec4<T> FromSpan(
+                              ReadOnlySpan<T> span
+        )
+        {
+            if (span.Length < 4) {
+                throw new ArgumentException("Span length must be at least 4.", nameof(span));
+            }
+
+            Unsafe.SkipInit(out Vec4<T> result);
+            Span<T> target = MemoryMarshal.CreateSpan(ref result.X, 4);
+            span[..4].CopyTo(target);
+            return result;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public ReadOnlySpan<T>.Enumerator GetEnumerator()
+        {
+            return AsSpan().GetEnumerator();
+        }
+
+        public override int GetHashCode()
+        {
+            return HashCode.Combine(X, Y, Z, W);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> operator -(
+                                       Vec4<T> left,
+                                       Vec4<T> right
+        )
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return Vec4.FromVector128Unsafe(left.AsVector128Unsafe() - right.AsVector128Unsafe());
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return Vec4.FromVector256Unsafe(left.AsVector256Unsafe() - right.AsVector256Unsafe());
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        public static bool operator !=(
+                                    Vec4<T> left,
+                                    Vec4<T> right
+        )
+        {
+            return !(left == right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> operator *(
+                                       Vec4<T> vector,
+                                       T scalar
+        )
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return Vec4.FromVector128Unsafe(vector.AsVector128Unsafe() * Vector128.Create(scalar));
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return Vec4.FromVector256Unsafe(vector.AsVector256Unsafe() * Vector256.Create(scalar));
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> operator *(
+                                       T scalar,
+                                       Vec4<T> vector
+        )
+        {
+            return vector * scalar;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static T operator *(
+                                 Vec4<T> left,
+                                 Vec4<T> right
+        )
+        {
+            return left.Dot(right);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> operator /(
+                                       Vec4<T> vector,
+                                       T scalar
+        )
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return Vec4.FromVector128Unsafe(vector.AsVector128Unsafe() / Vector128.Create(scalar));
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return Vec4.FromVector256Unsafe(vector.AsVector256Unsafe() / Vector256.Create(scalar));
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> operator +(
+                                       Vec4<T> left,
+                                       Vec4<T> right
+        )
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return Vec4.FromVector128Unsafe(left.AsVector128Unsafe() + right.AsVector128Unsafe());
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return Vec4.FromVector256Unsafe(left.AsVector256Unsafe() + right.AsVector256Unsafe());
+            }
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+        }
+
+        public static bool operator ==(
+                                    Vec4<T> left,
+                                    Vec4<T> right
+        )
+        {
+            if (Unsafe.SizeOf<T>() == 4) {
+                return left.AsVector128Unsafe() == right.AsVector128Unsafe();
+            }
+            if (Unsafe.SizeOf<T>() == 8) {
+                return left.AsVector256Unsafe() == right.AsVector256Unsafe();
+            }
+#pragma warning disable CA1065
+            throw GraphicVector.ThrowUnsupportedTypeException<T>();
+#pragma warning restore CA1065
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vector.GetLength.cs

@@ -0,0 +1,136 @@
+// -----------------------------------------------------------------------------
+// file="Vector.GetLength"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System.Runtime.CompilerServices;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    public static unsafe partial class Linear
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static T GetLength<T>(
+                        this Vec3<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return T.Sqrt(vec.Dot(vec));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static T GetLength<T>(
+                        this Vec2<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return T.Sqrt(vec.Dot(vec));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static T GetLength<T>(
+                        this Vec4<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return T.Sqrt(vec.Dot(vec));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float GetLength(
+                            this Vec2<int> vec
+        )
+        {
+            return Vec2.ConvertToSingle(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float GetLength(
+                            this Vec2<uint> vec
+        )
+        {
+            return Vec2.ConvertToSingle(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static double GetLength(
+                             this Vec2<long> vec
+        )
+        {
+            return Vec2.ConvertToDouble(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static double GetLength(
+                             this Vec2<ulong> vec
+        )
+        {
+            return Vec2.ConvertToDouble(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float GetLength(
+                            this Vec3<int> vec
+        )
+        {
+            return Vec3.ConvertToSingle(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float GetLength(
+                            this Vec3<uint> vec
+        )
+        {
+            return Vec3.ConvertToSingle(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static double GetLength(
+                             this Vec3<long> vec
+        )
+        {
+            return Vec3.ConvertToDouble(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static double GetLength(
+                             this Vec3<ulong> vec
+        )
+        {
+            return Vec3.ConvertToDouble(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float GetLength(
+                            this Vec4<int> vec
+        )
+        {
+            return Vec4.ConvertToSingle(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static float GetLength(
+                            this Vec4<uint> vec
+        )
+        {
+            return Vec4.ConvertToSingle(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static double GetLength(
+                             this Vec4<long> vec
+        )
+        {
+            return Vec4.ConvertToDouble(vec).GetLength();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static double GetLength(
+                             this Vec4<ulong> vec
+        )
+        {
+            return Vec4.ConvertToDouble(vec).GetLength();
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vector.Normalize.cs

@@ -0,0 +1,125 @@
+// -----------------------------------------------------------------------------
+// file="Vector.Normalize"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    public static partial class Linear
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec2<T> Normalize<T>(
+                              this Vec2<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return vec / T.Sqrt(vec.Dot(vec));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> Normalize<T>(
+                              this Vec3<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return vec / T.Sqrt(vec.Dot(vec));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> Normalize<T>(
+                              this Vec4<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return vec / T.Sqrt(vec.Dot(vec));
+        }
+
+        public static Vec2<float> Normalize(
+                                  this Vec2<int> vec
+        )
+        {
+            return Vec2.ConvertToSingle(vec).Normalize();
+        }
+
+        public static Vec2<float> Normalize(
+                                  this Vec2<uint> vec
+        )
+        {
+            return Vec2.ConvertToSingle(vec).Normalize();
+        }
+
+        public static Vec2<double> Normalize(
+                                   this Vec2<long> vec
+        )
+        {
+            return Vec2.ConvertToDouble(vec).Normalize();
+        }
+
+        public static Vec2<double> Normalize(
+                                   this Vec2<ulong> vec
+        )
+        {
+            return Vec2.ConvertToDouble(vec).Normalize();
+        }
+
+        public static Vec3<float> Normalize(
+                                  this Vec3<int> vec
+        )
+        {
+            return Vec3.ConvertToSingle(vec).Normalize();
+        }
+
+        public static Vec3<float> Normalize(
+                                  this Vec3<uint> vec
+        )
+        {
+            return Vec3.ConvertToSingle(vec).Normalize();
+        }
+
+        public static Vec3<double> Normalize(
+                                   this Vec3<long> vec
+        )
+        {
+            return Vec3.ConvertToDouble(vec).Normalize();
+        }
+
+        public static Vec3<double> Normalize(
+                                   this Vec3<ulong> vec
+        )
+        {
+            return Vec3.ConvertToDouble(vec).Normalize();
+        }
+
+        public static Vec4<float> Normalize(
+                                  this Vec4<int> vec
+        )
+        {
+            return Vec4.ConvertToSingle(vec).Normalize();
+        }
+
+        public static Vec4<float> Normalize(
+                                  this Vec4<uint> vec
+        )
+        {
+            return Vec4.ConvertToSingle(vec).Normalize();
+        }
+
+        public static Vec4<double> Normalize(
+                                   this Vec4<long> vec
+        )
+        {
+            return Vec4.ConvertToDouble(vec).Normalize();
+        }
+
+        public static Vec4<double> Normalize(
+                                   this Vec4<ulong> vec
+        )
+        {
+            return Vec4.ConvertToDouble(vec).Normalize();
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Graphics/Vector.Orthogonalize.cs

@@ -0,0 +1,63 @@
+// -----------------------------------------------------------------------------
+// file="Vector.Orthogonalize"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+namespace Nitload.Mathematics.Numeric.Graphics
+{
+    public static partial class Linear
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec4<T> Orthogonalize<T>(
+                              this Vec4<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return vec.Normalize();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static Vec3<T> Orthogonalize<T>(
+                              this Vec3<T> vec
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            return vec.Normalize();
+        }
+
+        /// <summary>
+        ///   Normalizes two vectors and makes <paramref name="target" /> perpendicular to <paramref
+        ///   name="source" />.
+        /// </summary>
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static void Orthogonalize<T>(
+                           this ref Vec4<T> source,
+                           ref Vec4<T> target
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            source = source.Normalize();
+            target -= source * target.Dot(source);
+            target = target.Normalize();
+        }
+
+        /// <summary>
+        ///   Normalizes two vectors and makes <paramref name="target" /> perpendicular to <paramref
+        ///   name="source" />.
+        /// </summary>
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static void Orthogonalize<T>(
+                           this ref Vec3<T> source,
+                           ref Vec3<T> target
+        ) where T : unmanaged, INumber<T>, IRootFunctions<T>
+        {
+            source = source.Normalize();
+            target -= source * target.Dot(source);
+            target = target.Normalize();
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/Exception.cs

@@ -0,0 +1,45 @@
+// -----------------------------------------------------------------------------
+// file="Exception"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Numerics;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public class MatrixDimensionMismatchException<T> : Exception where T : unmanaged, INumber<T>
+    {
+
+        public MatrixDimensionMismatchException(
+               Matrix<T> right,
+               Matrix<T> ans
+        )
+            : base($"""
+                    Matrix dimensions must agree but current:
+                    right({right?.ColumnLength ?? 0}×{right?.RowLength ?? 0}),
+                    ans({ans?.ColumnLength ?? 0}×{ans?.RowLength ?? 0})
+                    """
+                ) { }
+
+        public MatrixDimensionMismatchException(
+               Matrix<T> left,
+               Matrix<T> right,
+               Matrix<T> ans
+        )
+            : base($"""
+                    Matrix dimensions must agree but current:
+                    left({left?.ColumnLength ?? 0}×{left?.RowLength ?? 0}),
+                    right({right?.ColumnLength ?? 0}×{right?.RowLength ?? 0}),
+                    ans({ans?.ColumnLength ?? 0}×{ans?.RowLength ?? 0})
+                    """
+                ) { }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/IUnmanagedVector.cs

@@ -0,0 +1,26 @@
+// -----------------------------------------------------------------------------
+// file="IUnmanagedVector"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public interface IUnmanagedVector<TSelf> where TSelf : unmanaged, IUnmanagedVector<TSelf>
+    {
+
+        public static abstract TSelf Zero { get; }
+
+        public static abstract TSelf One { get; }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/Matrix.Create.cs

@@ -0,0 +1,90 @@
+// -----------------------------------------------------------------------------
+// file="Matrix.Create"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Mathematics;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+using Tlsf = global::Nitload.Common.HighPerformance.Memory.TwoLayerSegregatedFitAllocator;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public unsafe partial class Matrix<T>
+    {
+
+        private const int KernelSideLength = 4;
+
+        public static Matrix<T> Create(
+                                int columnLength,
+                                int rowLength,
+                                MatrixLayout layout,
+                                TwoLayerSegregatedFitAllocator allocator
+        )
+        {
+            ArgumentNullException.ThrowIfNull(allocator);
+
+            ArgumentOutOfRangeException.ThrowIfNegativeOrZero(columnLength);
+
+            ArgumentOutOfRangeException.ThrowIfNegativeOrZero(rowLength);
+
+            int kernelColumns = DivideRoundUp(columnLength, KernelSideLength);
+            int kernelRows = DivideRoundUp(rowLength, KernelSideLength);
+
+            if (layout is MatrixLayout.UpperTriangle or MatrixLayout.LowerTriangle)
+            {
+                if (kernelColumns != kernelRows) {
+                    throw new ArgumentException("Triangular layouts require square matrices.");
+                }
+            }
+
+            long kernelCount = layout switch
+            {
+                MatrixLayout.DenseRectangle => (long)kernelColumns * kernelRows,
+                MatrixLayout.UpperTriangle => TriangularKernelCount(kernelColumns),
+                MatrixLayout.LowerTriangle => TriangularKernelCount(kernelRows),
+                _ => throw new ArgumentOutOfRangeException(nameof(layout))
+            };
+
+            long requiredBytes = kernelCount * MatrixKernel<T>.Size;
+            ArgumentOutOfRangeException.ThrowIfNegativeOrZero(requiredBytes, nameof(requiredBytes));
+
+            Tlsf.Handle handle = allocator.Alloc((uint)requiredBytes);
+            MatrixStorageHLSF<T> storage = new MatrixStorageHLSF<T>(allocator, handle);
+
+            Matrix<T> matrix = new Matrix<T>(MatrixKernelBoarder.Create(layout), new MatrixSize(columnLength,
+                                                                                                rowLength,
+                                                                                                kernelColumns,
+                                                                                                kernelRows), storage);
+
+            return matrix;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static int DivideRoundUp(
+                           int value,
+                           int divisor
+        )
+        {
+            return (value + (divisor - 1)) / divisor;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static long TriangularKernelCount(
+                            int order
+        )
+        {
+            long n = order;
+            return (n * (n + 1L)) / 2L;
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/Matrix.cs

@@ -0,0 +1,184 @@
+// -----------------------------------------------------------------------------
+// file="Matrix"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric;
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+
+// using hlsfHandle = global::Nitload.Common.HighPerformance.Memory.HybridLayerSegregatedFitAllocationHandle;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public unsafe partial class Matrix<T> where T : unmanaged, INumber<T>
+    {
+
+        private protected static readonly int t_count = 16 / sizeof(T);
+
+        private protected static readonly int t_size = sizeof(T);
+        private readonly MatrixKernelBoarder _boarder;
+        private readonly MatrixSize _size;
+        private readonly MatrixStorage<T> _buffer;
+
+        private protected Matrix(
+                          MatrixKernelBoarder boarder,
+                          MatrixSize size,
+                          MatrixStorage<T> buffer
+
+        )
+        {
+            _boarder = boarder;
+            _size = size;
+            _buffer = buffer;
+        }
+
+        public T this[
+                 int x,
+                 int y
+        ]
+        {
+            get
+            {
+                if ((y < 0) ||
+                    (y > _size.ElementColumnLength) ||
+                    (x < 0) ||
+                    (x > _size.ElementRowLength)) {
+                    throw new ArgumentOutOfRangeException($"Index [{x},{y}] is out of matrix index range.");
+                }
+                long offset = GetElementOffset(x, y);
+                if (offset < 0) {
+                    return T.Zero;
+                }
+                return ((T*)Buffer)[offset];
+            }
+            set
+            {
+                if ((y < 0) ||
+                    (y > _size.ElementColumnLength) ||
+                    (x < 0) ||
+                    (x > _size.ElementRowLength)) {
+                    throw new ArgumentOutOfRangeException($"Index [{x},{y}] is out of matrix index range.");
+                }
+                long offset = GetElementOffset(x, y);
+                if (offset < 0) {
+                    return;
+                }
+                ((T*)Buffer)[offset] = value;
+            }
+        }
+
+        /// <summary>
+        ///   Matlab style indexer to generate a new matrix
+        /// </summary>
+        /// <param name="column">
+        ///
+        /// </param>
+        /// <param name="row">
+        ///
+        /// </param>
+        /// <returns>
+        ///
+        /// </returns>
+        /// <exception cref="NotImplementedException">
+        ///
+        /// </exception>
+        public unsafe Matrix<T> this[
+                                ReadOnlySpan<MatrixIndex> column,
+                                ReadOnlySpan<MatrixIndex> row
+        ]
+        {
+            get => throw new NotImplementedException();
+            set => throw new NotImplementedException();
+        }
+
+        /// <summary>
+        ///   Length of the matrix in columns.
+        /// </summary>
+        public int ColumnLength => _size.ElementColumnLength;
+
+        /// <summary>
+        ///   Length of the matrix in rows.
+        /// </summary>
+        public int RowLength => _size.ElementRowLength;
+
+        public MatrixLayout Layout { get; init; }
+
+        internal MatrixKernel<T>* Buffer => (MatrixKernel<T>*)_buffer.BufferPointer;
+
+        protected internal int KernelColumnLength => _size.KernelColumnLength;
+
+        protected internal int KernelRowLength => _size.KernelRowLength;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public unsafe ref MatrixKernel<T> GetMatrixKernel(
+                                          int x,
+                                          int y
+        )
+        {
+            return ref (_boarder.IsInBuffer(x, y, in _size) ?
+                        ref _buffer[_boarder.GetIndexOffset(x, y, in _size)] :
+                        ref Unsafe.NullRef<MatrixKernel<T>>());
+        }
+
+        internal long GetElementOffset(
+                      int x,
+                      int y
+        )
+        {
+            int
+                k_c = y / 4,
+                c = y % 4,
+                k_r = x / 4,
+                r = x % 4;
+            if (!_boarder.IsInBuffer(k_r, k_c, in _size)) {
+                return -1;
+            }
+            long kernelOffset = _boarder.GetIndexOffset(k_r, k_c, in _size);
+
+            return (kernelOffset * 16) + (4 * r) + c;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal unsafe ref T GetElementUnsafe(
+                              int x,
+                              int y
+        )
+        {
+            // Caller guarantees: indices in range, layout is dense square.
+            int k_c = y / 4;
+            int k_r = x / 4;
+            int c = y % 4;
+            int r = x % 4;
+
+            ref MatrixKernel<T> kernel = ref _buffer[_boarder.GetIndexOffset(k_r, k_c, in _size)];
+            return ref kernel[c, r];
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal unsafe MatrixKernel<T>* GetKernelAddressUnsafe(
+                                         int x,
+                                         int y
+        )
+        {
+            return Buffer + _boarder.GetIndexOffset(x, y, in _size);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal unsafe ref MatrixKernel<T> GetKernelUnsafe(
+                                            int x,
+                                            int y
+        )
+        {
+            // Caller guarantees: indices in range, layout is dense square.
+            int k_c = y / 4;
+            int k_r = x / 4;
+            return ref _buffer[_boarder.GetIndexOffset(k_c, k_r, in _size)];
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/MatrixCompute.Double.cs

@@ -0,0 +1,22 @@
+// -----------------------------------------------------------------------------
+// file="MatrixCompute.Double"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Mathematics.Numeric.Runtime;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+using unsafe MkProc = delegate* unmanaged[Cdecl]<Nitload.Mathematics.Numeric.Runtime.MatrixParallelTask*, void>;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static unsafe partial class MatrixCompute { }
+}

src/Nitload.Mathematics/Numeric/Linear/MatrixCompute.Single.cs

@@ -0,0 +1,61 @@
+// -----------------------------------------------------------------------------
+// file="MatrixCompute.Single"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Numeric.Runtime;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+using unsafe MkProc = delegate* unmanaged[Cdecl]<Nitload.Mathematics.Numeric.Runtime.MatrixParallelTask*, void>;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static unsafe partial class MatrixCompute
+    {
+
+        public static void Add(
+                           Matrix<float> left,
+                           Matrix<float> right,
+                           Matrix<float> ans
+        )
+        {
+            ArgumentNullException.ThrowIfNull(left, nameof(left));
+            ArgumentNullException.ThrowIfNull(right, nameof(right));
+            ArgumentNullException.ThrowIfNull(ans, nameof(ans));
+            if ((left.ColumnLength != right.ColumnLength) ||
+                (left.RowLength != right.RowLength) ||
+                (left.ColumnLength != ans.ColumnLength) ||
+                (left.RowLength != ans.RowLength)) {
+                throw new MatrixDimensionMismatchException<float>(
+                    left, right, ans
+                    );
+            }
+            MatrixElementType f32 = MatrixElementType.F32;
+            long count = ((long)left.KernelColumnLength) * left.KernelRowLength;
+
+            MatrixKernel<float>* leftPtr = left.Buffer;
+            MatrixKernel<float>* rightPtr = right.Buffer;
+            MatrixKernel<float>* ansPtr = ans.Buffer;
+            MkProc add = (MkProc)NumericalModule.GlobalContext.
+                mat_mk[f32.IntrinsicNode].
+                Get(MatrixIntrinsicHandle.Add);
+            MatrixParallelTask* t = stackalloc MatrixParallelTask[1];
+            /*
+            t->Mat1 = (MatrixKernel*)leftPtr;
+            t->Mat2 = (MatrixKernel*)rightPtr;
+            t->Mat3 = (MatrixKernel*)ansPtr;
+            t->ElementType = f32;
+            t->ComputeHandle = MatrixIntrinsicHandle.Add;
+            t->Scalar.Data<long>(0) = count;
+            /**/
+            SequentialOperation(t, count);
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/MatrixCompute.cs

@@ -0,0 +1,60 @@
+// -----------------------------------------------------------------------------
+// file="MatrixCompute"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Common.Internal;
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Runtime;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+using unsafe MkProc = delegate* unmanaged[Cdecl]<Nitload.Mathematics.Numeric.Runtime.MatrixParallelTask*, void>;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static unsafe partial class MatrixCompute
+    {
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal static unsafe MatrixKernel* GetKernelAddressUnsafe<T>(
+                                             MatrixKernel<T>* buffer,
+                                             int colLength,
+                                             int x,
+                                             int y
+        ) where T : unmanaged, INumber<T>
+        {
+            return (MatrixKernel*)(buffer + ((colLength * x) + y));
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static void SequentialOperation(
+                            MatrixParallelTask* task,
+                            long count
+        )
+        {
+            MatrixElementType elementType = task->ElementType;
+            MatrixParallelThread[] pool;
+            MkProc op = (MkProc)NumericalModule.GlobalContext.
+                mat_mk[elementType.IntrinsicNode].
+                Get(task->ComputeHandle);
+            if ((count <= 256) || !MatrixParallel.IsParallelMeaningful)
+            {
+                op(task);
+                return;
+            }
+            MatrixParallelWrap taskWrap = /**/ null/*MatrixParallel.ScheduleTask<BufferQueueLeftRightAns_NoNuma>(task)/**/;
+            taskWrap.LaunchParallel(SleepMode.DeepSleep);
+            taskWrap.Dispose();
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/MatrixIndex.cs

@@ -0,0 +1,61 @@
+// -----------------------------------------------------------------------------
+// file="MatrixIndex"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Numeric
+{
+    [StructLayout(LayoutKind.Sequential)]
+    public readonly struct MatrixIndex
+    {
+
+        public readonly int Begin;
+        public readonly int End;
+
+        private MatrixIndex(
+                int begin,
+                int end,
+                bool safe = false
+        )
+        {
+            Begin = begin;
+            End = end;
+        }
+
+        internal MatrixIndex(
+                 int begin,
+                 int end
+        )
+        {
+            if (begin > end) {
+                throw new ArgumentException("Begin index cannot be greater than end index.");
+            }
+
+            Begin = begin;
+            End = end;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static implicit operator MatrixIndex(
+                                        (int, int) range
+        )
+        {
+            return new(range.Item1, range.Item2, false);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static implicit operator MatrixIndex(
+                                        int index
+        )
+        {
+            return new(index, -1, false);
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/MatrixMetadata.cs

@@ -0,0 +1,234 @@
+// -----------------------------------------------------------------------------
+// file="MatrixMetadata"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public enum MatrixLayout
+    {
+
+        DenseRectangle,
+        UpperTriangle,
+        LowerTriangle,
+
+    }
+
+    public abstract class MatrixKernelBoarder
+    {
+
+        private static readonly MatrixKernelBoarder DenseSquareInstance = new DenseSquareKernelBoarder();
+        private static readonly MatrixKernelBoarder LowerTriangleInstance = new LowerTriangleKernelBoarder();
+        private static readonly MatrixKernelBoarder UpperTriangleInstance = new UpperTriangleKernelBoarder();
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static MatrixKernelBoarder Create(
+                                          MatrixLayout layout
+        )
+        {
+            return layout switch
+            {
+                MatrixLayout.DenseRectangle => DenseSquareInstance,
+                MatrixLayout.UpperTriangle => UpperTriangleInstance,
+                MatrixLayout.LowerTriangle => LowerTriangleInstance,
+                _ => throw new ArgumentOutOfRangeException(nameof(layout)),
+            };
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public abstract long GetIndexOffset(
+                             int xKernel,
+                             int yKernel,
+                             in MatrixSize size
+        );
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public abstract bool IsInBuffer(
+                             int xKernel,
+                             int yKernel,
+                             in MatrixSize enumerator
+        );
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        protected static long GetTriangularOrder(
+                              in MatrixSize enumerator
+        )
+        {
+            return Math.Min(enumerator.KernelColumnLength, enumerator.KernelRowLength);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        protected static bool IsWithinKernelDomain(
+                              in MatrixSize enumerator,
+                              int xKernel,
+                              int yKernel
+        )
+        {
+            return (xKernel >= 0) &&
+                   (yKernel >= 0) &&
+                   (xKernel < enumerator.KernelColumnLength) &&
+                   (yKernel < enumerator.KernelRowLength);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        protected static bool IsWithinTriangularDomain(
+                              in MatrixSize enumerator,
+                              int xKernel,
+                              int yKernel
+        )
+        {
+            long order = GetTriangularOrder(enumerator);
+            return (xKernel >= 0) && (yKernel >= 0) && (xKernel < order) && (yKernel < order);
+        }
+
+        private sealed class DenseSquareKernelBoarder : MatrixKernelBoarder
+        {
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public override long GetIndexOffset(
+                                 int xKernel,
+                                 int yKernel,
+                                 in MatrixSize enumerator
+            )
+            {
+                long x = xKernel;
+                long y = yKernel;
+                long col = enumerator.KernelColumnLength;
+
+                // Console.Write($"Attempt to get offset at ({x},{y}), offset is {x * col + y}  ");
+                return (x * col) + y;
+            }
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public override bool IsInBuffer(
+                                 int xKernel,
+                                 int yKernel,
+                                 in MatrixSize enumerator
+            )
+            {
+                return IsWithinKernelDomain(enumerator, xKernel, yKernel);
+            }
+
+        }
+
+        private sealed class LowerTriangleKernelBoarder : MatrixKernelBoarder
+        {
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public override long GetIndexOffset(
+                                 int xKernel,
+                                 int yKernel,
+                                 in MatrixSize enumerator
+            )
+            {
+                long order = GetTriangularOrder(enumerator);
+                long x = xKernel;
+                long y = yKernel;
+
+                long prefixColumns = (x * order) - ((x * (x - 1L)) / 2L);
+                return prefixColumns + (y - x);
+            }
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public override bool IsInBuffer(
+                                 int xKernel,
+                                 int yKernel,
+                                 in MatrixSize enumerator
+            )
+            {
+                return IsWithinTriangularDomain(enumerator, xKernel, yKernel) &&
+                       (yKernel >= xKernel);
+            }
+
+        }
+
+        private sealed unsafe class UpperTriangleKernelBoarder : MatrixKernelBoarder
+        {
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public override long GetIndexOffset(
+                                 int xKernel,
+                                 int yKernel,
+                                 in MatrixSize enumerator
+            )
+            {
+                long x = xKernel;
+                long y = yKernel;
+
+                long prefixColumns = (x * (x + 1L)) / 2L;
+                return prefixColumns + y;
+            }
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public override bool IsInBuffer(
+                                 int xKernel,
+                                 int yKernel,
+                                 in MatrixSize enumerator
+            )
+            {
+                return IsWithinTriangularDomain(enumerator, xKernel, yKernel) &&
+                       (yKernel <= xKernel);
+            }
+
+        }
+
+    }
+
+    [StructLayout(LayoutKind.Explicit)]
+    public struct MatrixSize
+    {
+
+        [FieldOffset(0)] public readonly int ElementColumnLength;
+        [FieldOffset(4)] public readonly int ElementRowLength;
+
+        [FieldOffset(8)] public readonly int KernelColumnLength;
+        [FieldOffset(12)] public readonly int KernelRowLength;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public MatrixSize(
+               int elementColumnLength,
+               int elementRowLength,
+               int kernelColumnLength,
+               int kernelRowLength
+        )
+        {
+            ElementColumnLength = elementColumnLength;
+            ElementRowLength = elementRowLength;
+            KernelColumnLength = kernelColumnLength;
+            KernelRowLength = kernelRowLength;
+        }
+
+    }
+
+    [StructLayout(LayoutKind.Explicit)]
+    public struct MatrixKernelEnumerator
+    {
+
+        [FieldOffset(0)] public int CurrentColumn;
+        [FieldOffset(4)] public int CurrentRow;
+
+        [FieldOffset(8)] public readonly int KernelColumnLength;
+        [FieldOffset(12)] public readonly int KernelRowLength;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public MatrixKernelEnumerator(
+               int kernelColumnLength,
+               int kernelRowLength,
+               int currentColumn = 0,
+               int currentRow = 0
+        )
+        {
+            CurrentColumn = currentColumn;
+            CurrentRow = currentRow;
+            KernelColumnLength = kernelColumnLength;
+            KernelRowLength = kernelRowLength;
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Linear/MatrixStorage.cs

@@ -0,0 +1,53 @@
+// -----------------------------------------------------------------------------
+// file="MatrixStorage"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Mathematics;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+using Tlsf = global::Nitload.Common.HighPerformance.Memory.TwoLayerSegregatedFitAllocator;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public abstract class MatrixStorage<T> : IDisposable where T : unmanaged, INumber<T>
+    {
+
+        internal unsafe ref MatrixKernel<T> this[
+                                            long index
+        ] => ref Unsafe.AsRef<MatrixKernel<T>>(BufferPointer + (MatrixKernel<T>.Size * index));
+
+        protected internal abstract unsafe T* BufferPointer { get; }
+
+        public abstract void Dispose();
+
+    }
+
+    internal class MatrixStorageHLSF<T>(
+                   Tlsf allocator,
+                   Tlsf.Handle handle
+    ) : MatrixStorage<T> where T : unmanaged, INumber<T>
+    {
+
+        internal readonly Tlsf.Handle _handle = handle;
+
+        internal readonly Tlsf _allocator = allocator;
+
+        protected internal override unsafe T* BufferPointer => (T*)(ulong)_handle.Address;
+
+        public override void Dispose()
+        {
+            allocator.Free(_handle);
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/NumericalModule.cs

@@ -0,0 +1,75 @@
+// -----------------------------------------------------------------------------
+// file="NumericalModule"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Mathematics.Internal;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public class NumericalModule : IStaticModule
+    {
+
+        public static unsafe int DefaultSIMDAlignment
+        {
+            get
+            {
+                SIMDID id = GlobalIntrinsicMask;
+                return (id.MaskByte[0] & 0b_00001111) switch
+                {
+                    1 => id.MaskByte[1] switch
+                    {
+                        1 => 16,
+                        2 => 32,
+                        3 => 64,
+                        4 => 64,
+                        _ => 8,
+                    },
+                    _ => 8,
+                };
+            }
+        }
+
+        public string ModuleName => "Math.Numerical";
+
+        internal static IntrinsicContext GlobalContext { get; set; } = new();
+
+        internal static SIMDID GlobalIntrinsicMask { get; set; }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static void InitializeModule()
+        {
+            MatrixParallel.Init();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void InitializeModule(
+                    IModuleInitializeConfiguration profile
+        )
+        {
+            try
+            {
+                MatrixParallel.Init();
+            }
+            finally
+            {
+                GC.Collect();
+            }
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void UninitializeModule()
+        {
+            MatrixParallel.Deinit();
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixOperatorPolicy.cs

@@ -0,0 +1,45 @@
+// -----------------------------------------------------------------------------
+// file="MatrixOperatorPolicy"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric.Runtime
+{
+    internal enum MatrixCachePolicy : short
+    {
+
+        None = 0,
+
+        #region usage
+
+        Left = 1,
+        Right = 2,
+        Ans = 3,
+
+        #endregion
+
+        #region index direction
+
+        RowMajor = 4,
+        ColMajor = 5,
+        Sequential = 6,
+
+        #endregion
+
+        #region shape
+
+        Square = 8,
+        UpperTriangle = 9,
+        LowerTriangle = 10,
+
+        #endregion
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallel.Await.old

@@ -0,0 +1,157 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixParallel.Await"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using System;
+using System.Collections.Concurrent;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static partial class MatrixParallel
+    {
+
+        private static AwaitingQueue _awaitingLeaderHandle;
+
+        private static volatile int _createGeneration;
+        private static int _sortState;
+
+        public static MatrixParallelThread[] LeadParallel(
+                                             long recommendCount
+        )
+        {
+            int recommend = (recommendCount > int.MaxValue) ? int.MaxValue : ((int)recommendCount);
+            int min = (int)long.Min(recommend, _threads.Capacity);
+            min /= 2;
+            MatrixParallelThread[] ret;
+            if (_threads.Count < min)
+            {
+                WaitHandle handle = new(
+                    _createGeneration,
+                    min,
+                    recommend
+                );
+                do
+                {
+                    if (Interlocked.Exchange(ref _sortState, 1) == 0)
+                    {
+                        _awaitingLeaderHandle.Enqueue(handle);
+                        _ = Interlocked.Exchange(ref _sortState, 0);
+                        _ = Interlocked.Increment(ref _createGeneration);
+                        break;
+                    }
+                } while (true);
+                int cycleGeneration = _createGeneration;
+                do
+                {
+                    Thread.Sleep(0);
+
+                    if (_createGeneration != cycleGeneration)
+                    {
+                        if (Interlocked.Exchange(ref _sortState, 1) == 0)
+                        {
+                            if (_awaitingLeaderHandle.TryPopIfRefEqual(cycleGeneration, in handle))
+                            {
+                                // Sort is implemented in TryPopIfRefEqual
+                                ret = _threads.PopRange(handle.BestCount);
+                                _createGeneration++;
+                                _ = Interlocked.Exchange(ref _sortState, 0);
+                                return ret;
+                            }
+                            _ = Interlocked.Exchange(ref _sortState, 0);
+                        } else
+                        {
+                            cycleGeneration = _createGeneration;
+                        }
+                    }
+                } while (true);
+            } else
+            {
+                // directly get the leader
+                while (true)
+                {
+                    if (Interlocked.Exchange(ref _sortState, 1) == 0)
+                    {
+                        ret = _threads.PopRange(recommend);
+                        _createGeneration++;
+                        _ = Interlocked.Exchange(ref _sortState, 0);
+                        return ret;
+                    }
+                }
+            }
+
+            // return ret;
+        }
+
+        public static void ReturnParallelResources(
+                           ref MatrixParallelThread[] pool
+        )
+        {
+            while (true)
+            {
+                if (Interlocked.Exchange(ref _sortState, 1) == 0)
+                {
+                    _threads.PushRange(pool);
+                    _ = Interlocked.Exchange(ref _sortState, 0);
+                    break;
+                }
+            }
+
+            pool = null!;
+        }
+
+        private record class WaitHandle(
+                             int generation,
+                             int least,
+                             int best
+        )
+        {
+
+            public int CreateGeneration { get; init; } = generation;
+
+            public int LeastCount { get; init; } = least;
+
+            public int BestCount { get; init; } = best;
+
+            public int Priority { get; private set; }
+
+            [MethodImpl(MethodImplOptions.AggressiveInlining)]
+            public void UpdatePriority(
+                        int currentGeneration
+            )
+            {
+                Priority = (currentGeneration - CreateGeneration) * (MatrixParallel._threads.Capacity - LeastCount);
+            }
+
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallel.AwaitingQueue.old

@@ -0,0 +1,169 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixParallel.AwaitingQueue"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using System;
+using System.Collections.Concurrent;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static partial class MatrixParallel
+    {
+        /*
+         * Because caller of this class has its own sync object
+         * this queue will use no sync object internally.
+        */
+
+        private class AwaitingQueue
+        {
+
+            private Element _quickElements;
+            private volatile int _count;
+            private readonly Queue<WaitHandle> _fallback = [];
+
+            public AwaitingQueue() { }
+
+            public void Enqueue(
+                        WaitHandle handle
+            )
+            {
+                if (_count > 7)
+                {
+                    _fallback.Enqueue(handle);
+                    _ = Interlocked.Increment(ref _count);
+                    return;
+                }
+                _quickElements[_count] = handle;
+                _ = Interlocked.Increment(ref _count);
+            }
+
+            public bool TryDequeue(
+                        int generation,
+                        out WaitHandle handle
+            )
+            {
+                if (_count == 0)
+                {
+                    handle = null!;
+                    return false;
+                }
+                for (int i = 0; i < _count; i++) {
+                    _quickElements[i].UpdatePriority(generation);
+                }
+                Element.Sort(ref _quickElements);
+                handle = _quickElements[0];
+                _count--;
+                _quickElements[0] = default!;
+                if (_count > 7) {
+                    _ = _fallback.TryDequeue(out _quickElements[0]!);
+                }
+                return true;
+            }
+
+            public bool TryPeek(
+                        int generation,
+                        out WaitHandle handle
+            )
+            {
+                if (_count == 0)
+                {
+                    handle = null!;
+                    return false;
+                }
+                for (int i = 0; i < _count; i++) {
+                    _quickElements[i].UpdatePriority(generation);
+                }
+
+                Element.Sort(ref _quickElements);
+                handle = _quickElements[0];
+                return true;
+            }
+
+            public bool TryPopIfRefEqual(
+                        int generation,
+                        in WaitHandle handle
+            )
+            {
+                if (_count == 0) {
+                    return false;
+                }
+                for (int i = 0; i < _count; i++) {
+                    _quickElements[i].UpdatePriority(generation);
+                }
+                Element.Sort(ref _quickElements);
+                if (!ReferenceEquals(_quickElements[0], handle)) {
+                    return false;
+                }
+                _count--;
+                _quickElements[0] = null!;
+                if (_count > 7) {
+                    _ = _fallback.TryDequeue(out _quickElements[0]!);
+                }
+                return true;
+            }
+
+        }
+
+        [InlineArray(8)]
+        private struct Element
+        {
+
+            private WaitHandle value;
+
+            public static void Sort(
+                               ref Element e
+            )
+            {
+                int maxIdx = 0;
+                WaitHandle first = e[0];
+                int maxPr = first?.Priority ?? int.MinValue;
+
+                for (int i = 1; i < 8; i++)
+                {
+                    WaitHandle cur = e[i];
+                    int p = cur?.Priority ?? int.MinValue;
+                    if (p > maxPr)
+                    {
+                        maxPr = p;
+                        maxIdx = i;
+                    }
+                }
+
+                if (maxIdx != 0) {
+                    (e[0], e[maxIdx]) = (e[maxIdx], e[0]);
+                }
+            }
+
+        }
+
+    }
+}
+

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallel.Numa.cs

@@ -0,0 +1,44 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallel.Numa"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.Collections;
+using Nitload.Common.HighPerformance.Memory;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static unsafe partial class MatrixParallel
+    {
+
+        private static L4[] _numaBuffer;
+        private static NumaMeasurement[] _numaMeasurements;
+
+        private static int[] _numaPayload;
+        private static MatrixParallelThread[][] _numaSyncHeavy;
+        private static MatrixParallelThread[][] _numaWaveFront;
+        private static MatrixParallelThread[] _threads;
+
+        private static void UpdateThreadPayload()
+        {
+            for (int i = 0; i < _numaPayload.Length; i++) { }
+        }
+
+        private struct NumaMeasurement
+        {
+
+            public uint ActivePredictors;
+            public uint ActiveThreadCount;
+            public ulong ActiveTasks;
+
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallel.Publish.cs

@@ -0,0 +1,72 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallel.Publish"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.Collections;
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Mathematics;
+using Nitload.Mathematics.Numeric.Runtime;
+using System;
+using System.Collections.Generic;
+using System.ComponentModel;
+using System.Data;
+using System.Linq;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static partial class MatrixParallel
+    {
+
+        private const int MIN_TASK_PER_THREAD = 256;
+        private const int REC_TASK_PER_THREAD = 1024;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+
+        internal static unsafe MatrixParallelWrap ScheduleTask<TTaskQueue>(
+                                                  MatrixParallelTask* package
+        )
+
+        {
+            /*
+            int primCount = package->PrimCount;
+            int secCount = package->SecCount;
+            long count = (long)primCount * secCount;
+            /**/
+
+            // Console.WriteLine($"amount of kernels is {count}");
+            long payload = 0;
+            int payloadIndex = 0;
+            ulong additionalPayload = 0;
+            for (int i = 0; i < _numaMeasurements.Length; i++)
+            {
+                MatrixParallelThread[] numa = _numaWaveFront[i];
+                if ((numa is null) || (numa.Length == 0))
+                {
+                    continue;
+                }
+                ulong total = (ulong)_numaWaveFront[i].Length;
+                NumaMeasurement measure = _numaMeasurements[i];
+                ulong numaPayload = measure.ActiveTasks;
+
+                // uint numaActive = measure.ActiveThreadCount;
+                long curPayload = ((long)total) - (((long)(numaPayload + additionalPayload + 63)) / 64);
+                if (curPayload < payload)
+                {
+                    payload = curPayload;
+                    payloadIndex = i;
+                }
+            }
+            L4 cache = _numaBuffer[payloadIndex];
+            MatrixParallelThread[] selectedNuma = _numaWaveFront[payloadIndex];
+            throw new NotImplementedException();
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallel.cs

@@ -0,0 +1,199 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallel"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.Collections;
+using Nitload.Common.HighPerformance.ProcessorAbstraction;
+using System.Collections.Generic;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public static unsafe partial class MatrixParallel
+    {
+
+        private const int RESERVED_THREAD_COUNT = 4;
+
+        private static bool _meaningful;
+
+        private static CapacityFixedStack<MatrixParallelThread> _retiredThreads;
+
+        public static bool IsParallelSupported { get; } = Environment.ProcessorCount > 6;
+
+        public static bool IsParallelMeaningful => IsParallelSupported &&
+                                                   IsInitialized &&
+                                                   _meaningful;
+
+        public static int ContextCount { get; private set; }
+
+        private static bool IsInitialized { get; set; }
+
+        internal static void Deinit()
+        {
+            foreach (MatrixParallelThread t in _threads)
+            {
+                t.Pause();
+                t.Dispose();
+            }
+            foreach (MatrixParallelThread t in _retiredThreads) {
+                t.Dispose();
+            }
+            ContextCount = 0;
+            IsInitialized = false;
+        }
+
+        internal static void Init()
+        {
+            Dictionary<ulong, int> contextMap = [];
+
+            if (IsInitialized) {
+                return;
+            }
+            if (!IsParallelSupported)
+            {
+                IsInitialized = true;
+                return;
+            }
+            contextMap = [];
+            int count = Environment.ProcessorCount;
+            SIMDID globalSimd = new(0);
+            _retiredThreads = new(RESERVED_THREAD_COUNT);
+
+            // init thread on all logical cores
+            _numaWaveFront = Numa.CreateForeachThread<MatrixParallelThread>(
+               core => {
+                   MatrixParallelThread thread = new(core);
+                   SpinWait.SpinUntil(() => thread.State == MatrixThreadState.Ready);
+
+                   SIMDID simd = thread.Simd;
+                   ulong mask = simd.Mask;
+                   contextMap[mask] = contextMap.TryGetValue(mask, out int contextCount) ?
+                                      (contextCount + 1) :
+                                      1;
+                   globalSimd = (globalSimd.Mask != 0) ? UniteSIMDID(globalSimd, simd) : simd;
+                   return thread;
+               });
+
+            // make global simd  context
+            IntrinsicContext globalContext = new();
+            LoadIntrinsicFunction(&globalContext, globalSimd.Mask);
+            NumericalModule.GlobalContext = globalContext;
+            NumericalModule.GlobalIntrinsicMask = globalSimd;
+            ContextCount = contextMap.Count;
+
+            // retire threads
+            if (ContextCount == 1)
+            {
+                MatrixParallelThread[] threadArr = _numaWaveFront[0];
+                int retireCount = Math.Min(RESERVED_THREAD_COUNT, Math.Max(0, threadArr.Length));
+                for (int i = 1; i <= retireCount; i++) {
+                    RetireThread(threadArr[^i]);
+                }
+                _numaWaveFront[0] = threadArr[0..^retireCount];
+            } else
+            {
+                RetireWhenMultipleSIMDID(_numaWaveFront);
+            }
+            _meaningful = _threads.Length > 1;
+
+            // _awaitingLeaderHandle = new AwaitingQueue();
+            IsInitialized = true;
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static void RetireThread(
+                            MatrixParallelThread thread
+        )
+        {
+            thread.Pause();
+            thread.Dispose();
+            _retiredThreads.Push(thread);
+        }
+
+        private static void RetireWhenMultipleSIMDID(
+                            MatrixParallelThread[][] threads
+        )
+        {
+            MatrixParallelThread[] allThreads = new MatrixParallelThread[Environment.ProcessorCount];
+            int offset = 0;
+            foreach (MatrixParallelThread[] threadArr in threads)
+            {
+                Array.Copy(threadArr, 0, allThreads, offset, threadArr.Length);
+                offset += threadArr.Length;
+            }
+            Array.Sort(allThreads, static(
+                                   left,
+                                   right
+            ) => CompareSIMDID(left.Simd, right.Simd));
+            int retireCount = Math.Min(RESERVED_THREAD_COUNT, Math.Max(0, allThreads.Length));
+
+            MatrixParallelThread[] retireArr = new MatrixParallelThread[retireCount];
+            for (int i = 0; i < retireCount; i++)
+            {
+                MatrixParallelThread thread = allThreads[i];
+                retireArr[i] = thread;
+                RetireThread(thread);
+            }
+
+            // rebuild numa active arrays
+            for (int numa = 0; numa < _numaWaveFront.Length; numa++)
+            {
+                MatrixParallelThread[] source = threads[numa];
+                int keepCount = 0;
+
+                for (int i = 0; i < source.Length; i++)
+                {
+                    if (!retireArr.Contains(source[i])) {
+                        keepCount++;
+                    }
+                }
+
+                MatrixParallelThread[] active = new MatrixParallelThread[keepCount];
+                int index = 0;
+                for (int i = 0; i < source.Length; i++)
+                {
+                    MatrixParallelThread thread = source[i];
+                    if (!retireArr.Contains(thread)) {
+                        active[index++] = thread;
+                    }
+                }
+
+                _numaWaveFront[numa] = active;
+            }
+
+            _threads = allThreads[retireCount..];
+        }
+
+        #region C invoke
+
+        [LibraryImport(Native.LibName, EntryPoint = "sn_compare_simd")]
+        [UnmanagedCallConv(CallConvs = [typeof(CallConvCdecl)])]
+        [DefaultDllImportSearchPaths(DllImportSearchPath.AssemblyDirectory)]
+        private static partial int CompareSIMDID(
+                                   SIMDID left,
+                                   SIMDID right
+        );
+
+        [LibraryImport(Native.LibName, EntryPoint = "load_intrinsic_function")]
+        [UnmanagedCallConv(CallConvs = [typeof(CallConvCdecl)])]
+        [DefaultDllImportSearchPaths(DllImportSearchPath.AssemblyDirectory)]
+        private static partial void LoadIntrinsicFunction(
+                                    IntrinsicContext* context,
+                                    ulong id
+        );
+
+        [LibraryImport(Native.LibName, EntryPoint = "sn_get_unite_simd")]
+        [UnmanagedCallConv(CallConvs = [typeof(CallConvCdecl)])]
+        [DefaultDllImportSearchPaths(DllImportSearchPath.AssemblyDirectory)]
+        private static partial SIMDID UniteSIMDID(
+                                      SIMDID left,
+                                      SIMDID right
+        );
+
+    #endregion
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallelQueue.cs

@@ -0,0 +1,77 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallelQueue"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Microsoft.VisualBasic.FileIO;
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Mathematics.Numeric.Runtime;
+using System;
+using System.Collections.Concurrent;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices.JavaScript;
+using System.Text;
+using System.Threading;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    internal unsafe class MatrixParallelQueue : IDisposable
+    {
+
+        private MatrixParallelTask[] _currentPack;
+        private MatrixParallelTask[] _predictBuffer;
+        private bool disposedValue;
+
+        private int _currentCapacity;
+        private int _cursor;
+        private int _predictCapacity;
+
+        private ManualResetEventSlim _exchangeLock = new();
+
+        public MatrixParallelQueue()
+        {
+            _currentPack = new MatrixParallelTask[512];
+            _predictBuffer = new MatrixParallelTask[512];
+        }
+
+        public void Dispose()
+        {
+            Dispose(disposing:true);
+            GC.SuppressFinalize(this);
+        }
+
+        public bool TryGetTask(
+                    out MatrixParallelTask task
+        )
+        {
+            if (_cursor >= _currentCapacity)
+            {
+                task = default;
+                return false;
+            }
+            task = _currentPack[_cursor++];
+            return true;
+        }
+
+        protected virtual void Dispose(
+                               bool disposing
+        )
+        {
+            if (!disposedValue)
+            {
+                if (disposing) {
+                    _exchangeLock.Dispose();
+                }
+
+                disposedValue = true;
+            }
+        }
+
+    }
+}
+

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallelTask.cs

@@ -0,0 +1,70 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallelTask"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Common.HighPerformance.ProcessorAbstraction;
+using System;
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+
+namespace Nitload.Mathematics.Numeric.Runtime
+{
+    [StructLayout(LayoutKind.Sequential)]
+    internal struct MatrixParallelSource
+    {
+
+        public MatrixParallelTask Task;
+        public ulong Count;
+        public ulong Cursor;
+
+    }
+
+    [StructLayout(LayoutKind.Explicit, Size = 128)]
+    internal unsafe struct MatrixParallelTask
+    {
+
+        [FieldOffset(104)] private ulong _reserved0;
+
+        [FieldOffset(0)] public L4.Handle Buffer0;
+        [FieldOffset(16)] public L4.Handle Buffer1;
+        [FieldOffset(32)] public L4.Handle Buffer2;
+        [FieldOffset(48)] public L4.Handle Buffer3;
+        [FieldOffset(112)] public M128 ScalarPack;
+        [FieldOffset(64)] public M512 Profile;
+        [FieldOffset(64)] public MatrixElementType ElementType;         //记录了类型掩码的类型位宽的union，大小4byte
+        [FieldOffset(68)] public MatrixIntrinsicHandle ComputeHandle;    //记录了函数指针在类型组中的索引，实际就是int，但是为每一个索引起了名字，于是变成了enum，大小4byte
+
+        public MatrixParallelTask()
+        {
+            Unsafe.SkipInit(out this);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static void Copy(
+                           MatrixParallelTask* source,
+                           MatrixParallelTask* target
+        )
+        {
+            Buffer.MemoryCopy(source, target, sizeof(MatrixParallelTask),
+                              sizeof(MatrixParallelTask));
+        }
+
+        #region offset and length settings
+
+        [FieldOffset(72)] public int Offset0;
+        [FieldOffset(76)] public int Count0;
+        [FieldOffset(80)] public int Offset1;
+        [FieldOffset(84)] public int Count1;
+        [FieldOffset(88)] public int Offset2;
+        [FieldOffset(92)] public int Count2;
+        [FieldOffset(96)] public int Offset3;
+        [FieldOffset(100)] public int Count3;
+
+    #endregion
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallelThread.Fma.cs.old

@@ -0,0 +1,106 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixParallelThread.Fma"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using Nitload.Common.Internal;
+using Nitload.Mathematics;
+using Nitload.Mathematics.Internal;
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public partial class MatrixParallelThread
+    {
+
+        internal static unsafe void ComputeFmaPanelF32(
+                                    MatrixParallelTask* p
+        )
+        {
+            int f32 = MatrixElementType.F32.IntrinsicNode;
+            MatrixKernel* leftBuffer = (MatrixKernel*)p->Left;
+            MatrixKernel* rightBuffer = (MatrixKernel*)p->Right;
+            MatrixKernel* ansBuffer = (MatrixKernel*)p->Answer;
+            ref long offsetRef = ref Unsafe.As<int, long>(ref p->PrimTileOffset);
+            ref long lengthRef = ref Unsafe.As<int, long>(ref p->PrimCount);
+            int ansWidth = p->AnsPrimOffset;
+            int ansHeight = p->AnsSecOffset;
+            int commonK = p->LeftPrimOffset;
+
+            M128* enumeration = stackalloc M128[1];
+            enumeration->Member<int>(2) = commonK;
+            enumeration->Member<int>(3) = ansHeight;
+
+            // in this occasion size must be 1
+            for (long c = offsetRef; c < offsetRef + lengthRef; c++)
+            {
+                int i = (int)(c / ansHeight);
+                int j = (int)(c % ansHeight);
+                enumeration->Member<int>(0) = i;
+                enumeration->Member<int>(1) = j;
+                NumericalModule.GlobalContext.mat_sk[f32].kernel_tile_fma(leftBuffer, rightBuffer,
+                                                                         ansBuffer, enumeration);
+            }
+        }
+
+        internal static unsafe void ComputeFmaPanelF64(
+                                    MatrixParallelTask* p
+        )
+        {
+            int f64 = MatrixElementType.F64.IntrinsicNode;
+
+            // int size = 1;
+            MatrixKernel* leftBuffer = (MatrixKernel*)p->Left;
+            MatrixKernel* rightBuffer = (MatrixKernel*)p->Right;
+            MatrixKernel* ansBuffer = (MatrixKernel*)p->Answer;
+            ref long offsetRef = ref Unsafe.As<int, long>(ref p->PrimTileOffset);
+            ref long lengthRef = ref Unsafe.As<int, long>(ref p->PrimCount);
+            int ansWidth = p->AnsPrimOffset;
+            int ansHeight = p->AnsSecOffset;
+            int commonK = p->LeftPrimOffset;
+
+            M128* enumeration = stackalloc M128[1];
+            enumeration->Member<int>(2) = commonK;
+            enumeration->Member<int>(3) = ansHeight;
+
+            // in this occasion size must be 1
+            for (long c = offsetRef; c < offsetRef + lengthRef; c++)
+            {
+                int i = (int)(c / ansHeight);
+                int j = (int)(c % ansHeight);
+                enumeration->Member<int>(0) = i;
+                enumeration->Member<int>(1) = j;
+                NumericalModule.GlobalContext.mat_sk[f64].kernel_tile_fma(leftBuffer, rightBuffer,
+                                                                         ansBuffer, enumeration);
+            }
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallelThread.Model.cs.old

@@ -0,0 +1,76 @@
+//-----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+// file="MatrixParallelThread.Model"
+// Project: StgSharp
+// AuthorGroup: Nitload
+// Copyright (c) Nitload. All rights reserved.
+//     
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+//     
+// -----------------------------------------------------------------------
+// -----------------------------------------------------------------------
+using Nitload.Mathematics;
+using Nitload.Mathematics.Internal;
+using System;
+using System.Runtime.CompilerServices;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public partial class MatrixParallelThread
+    {
+
+        private const int KernelStealCount = 64;
+
+        internal unsafe void BufferCompute(
+                             MatrixParallelTask* p
+        )
+        {
+            // while () { }
+            int size = p->ElementSize * 16;
+            long count = Unsafe.As<int, long>(ref p->PrimCount);
+            long offset = Unsafe.As<int, long>(ref p->PrimTileOffset) * size;
+            ScalarPacket* scalar = p->Scalar;
+            IntPtr computeHandle = _context[p->ComputeHandle];
+            ((delegate* unmanaged[Cdecl]<void*, void*, void*, void*, long, void>)computeHandle)((void*)(p->Left + offset),
+                                                                                                (void*)(p->Right + offset),
+                                                                                                (void*)(p->Answer + offset),
+                                                                                                scalar,
+                                                                                                count);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal unsafe void SpecialBufferCompute(
+                             MatrixParallelTask* p
+        )
+        {
+            switch ((MatrixElementType)p->ElementSize)
+            {
+                case MatrixElementType.F64:
+                    ComputeFmaPanelF64(p);
+                    return;
+                case MatrixElementType.F32:
+                    ComputeFmaPanelF32(p);
+                    break;
+                default:
+                    break;
+            }
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallelThread.Sequential.cs

@@ -0,0 +1,17 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallelThread.Sequential"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric.Runtime
+{
+    internal class MatrixParallelThread { }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallelThread.cs

@@ -0,0 +1,179 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallelThread"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.HighPerformance.ProcessorAbstraction;
+using Nitload.Common.Internal;
+using Nitload.Mathematics.Numeric.Runtime;
+using Nitload.Threading;
+using System;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Threading;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    internal enum MatrixThreadState : int
+    {
+
+        Born = 0,
+        Ready = 1,
+        Running = 2,
+        Idle = 4,
+        Retired = 8,
+
+    }
+
+    public unsafe partial class MatrixParallelThread : IEquatable<MatrixParallelThread>, IDisposable
+    {
+
+        private IntrinsicContext* _context;
+        private readonly MatrixParallelTask* _currentTask;
+        private volatile bool _disposed;
+        private readonly ManualResetEventSlim _resetEvent;
+
+        private volatile MatrixThreadState _target;
+        private Thread _managedThread;
+
+        internal MatrixParallelThread(
+                 LogicalCore binded
+        )
+        {
+            _resetEvent = new ManualResetEventSlim(false);
+            _managedThread = new Thread(MatrixParallelWorkLoad);
+            _managedThread.Start();
+            BindedLogicalCore = binded;
+        }
+
+        public LogicalCore BindedLogicalCore { get; init; }
+
+        public ThreadPriority AfterWork { get; set; }
+
+        internal MatrixThreadState State { get; private set; }
+
+        internal SIMDID Simd { get; private set; }
+
+        public void Dispose()
+        {
+            Dispose(disposing:true);
+            GC.SuppressFinalize(this);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void Restart()
+        {
+            _resetEvent.Set();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal void Pause()
+        {
+            ObjectDisposedException.ThrowIf(_disposed, nameof(MatrixParallelThread));
+            if (State == MatrixThreadState.Idle) {
+                return;
+            }
+            _resetEvent.Set();
+            SpinWait.SpinUntil(() => State == MatrixThreadState.Idle);
+            _managedThread.Join();
+        }
+
+        protected virtual void Dispose(
+                               bool disposing
+        )
+        {
+            if (!_disposed)
+            {
+                if (disposing)
+                {
+                    _disposed = true;
+                    _resetEvent.Set();
+                    _managedThread.Join();
+                    _resetEvent.Dispose();
+                }
+            }
+        }
+
+        private unsafe void MatrixParallelWorkLoad()
+        {
+            State = MatrixThreadState.Born;
+
+            // Bind OS-thread affinity once when the worker starts running
+            ThreadHelper.BindToLogicalCore(BindedLogicalCore);
+            _ = Thread.Yield();
+            if (Simd.Mask != 0)
+            {
+                SIMDID feature = new(GetSimdId());
+                IntrinsicContext* contextPtr = stackalloc IntrinsicContext[1];
+                LoadIntrinsicFunction(contextPtr, feature.Mask);
+                Simd = feature;
+                _context = contextPtr;
+            }
+            MatrixParallelTask* taskPtr = stackalloc MatrixParallelTask[4];
+            State = MatrixThreadState.Ready;
+            while (true)
+            {
+                // Thread self = Thread.CurrentThread;
+                _resetEvent.Wait();
+                _resetEvent.Reset();
+                if (_disposed)
+                {
+                    break;
+                }
+                if (_target == MatrixThreadState.Idle)
+                {
+                    State = MatrixThreadState.Idle;
+                    continue;
+                }
+                State = MatrixThreadState.Idle;
+                State = MatrixThreadState.Running;
+                _managedThread.Priority = ThreadPriority.Highest;
+
+
+                // MatrixParallelQueue taskQueue = BindedWrap.TaskQueue;
+                /*
+                if (0 > (int)p ->ComputeHandle) {
+                    SpecialBufferCompute(p);
+                }
+                switch (p->ComputeMode)
+                {
+                    case MatrixIndexStyle.BUFFER_INDEX:
+                        BufferCompute(p);
+                        break;
+                    default:
+                        break;
+                }
+                /**/
+
+                _managedThread.Priority = AfterWork;
+                State = MatrixThreadState.Idle;
+            }
+            State = MatrixThreadState.Retired;
+        }
+
+#pragma warning disable CA5393
+        [LibraryImport(Native.LibName, EntryPoint = "sn_get_simd_level_local")]
+        [UnmanagedCallConv(CallConvs = [typeof(CallConvCdecl)])]
+        [DefaultDllImportSearchPaths(DllImportSearchPath.AssemblyDirectory)]
+        internal static partial ulong GetSimdId();
+
+        [LibraryImport(Native.LibName, EntryPoint = "load_intrinsic_function")]
+        [UnmanagedCallConv(CallConvs = [typeof(CallConvCdecl)])]
+        [DefaultDllImportSearchPaths(DllImportSearchPath.AssemblyDirectory)]
+        private static unsafe partial void LoadIntrinsicFunction(
+                                           IntrinsicContext* context,
+                                           ulong id
+        );
+
+        public bool Equals(
+                    MatrixParallelThread? other
+        )
+        {
+            return (other is not null) && ReferenceEquals(this, other);
+        }
+#pragma warning restore CA5393
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/MatrixParallelWrap.cs

@@ -0,0 +1,107 @@
+// -----------------------------------------------------------------------------
+// file="MatrixParallelWrap"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using Nitload.Common.Collections;
+using Nitload.Common.HighPerformance.Memory;
+using Nitload.Mathematics.Numeric.Runtime;
+using System;
+using System.Collections.Generic;
+using System.ComponentModel.DataAnnotations;
+using System.Diagnostics.CodeAnalysis;
+using System.Linq;
+using System.Runtime.CompilerServices;
+using System.Text;
+using System.Threading;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    public enum SleepMode : int
+    {
+
+        QuickWakeup = ThreadPriority.AboveNormal,
+        NormalSleep = ThreadPriority.BelowNormal,
+        DeepSleep = ThreadPriority.Lowest,
+
+    }
+
+    public class MatrixParallelWrap : IDisposable
+    {
+
+        private MatrixParallelThread[] _threads;
+        private bool _disposed;
+        private readonly CountdownEvent _countdownEvent;
+
+        internal MatrixParallelWrap(
+                 MatrixParallelThread[] threadArray
+        )
+        {
+            _threads = threadArray;
+            _countdownEvent = new CountdownEvent(threadArray.Length);
+        }
+
+        public ReadOnlySpan<MatrixParallelThread> Threads => _disposed ? [] : _threads;
+
+        public ThreadPriority AfterWork { get; private set; }
+
+        internal L4 DataCache { get; set; }
+
+        public void Dispose()
+        {
+            Dispose(true);
+
+            GC.SuppressFinalize(this);
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public void LaunchParallel(
+                    SleepMode afterWork
+        )
+        {
+            ObjectDisposedException.ThrowIf(_disposed, nameof(MatrixParallelWrap));
+            AfterWork = (ThreadPriority)afterWork;
+            foreach (MatrixParallelThread thread in _threads)
+            {
+                thread.AfterWork = AfterWork;
+
+                // thread.BindedWrap = this;
+            }
+            WaitAllWraps();
+        }
+
+        public void WaitAllWraps()
+        {
+            if (_disposed) {
+                return;
+            }
+            _countdownEvent.Wait();
+        }
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        internal void ReportThreadAccomplished()
+        {
+            ObjectDisposedException.ThrowIf(_disposed, nameof(MatrixParallelWrap));
+            _ = _countdownEvent.Signal();
+        }
+
+        private void Dispose(
+                     bool disposing
+        )
+        {
+            if (_disposed) {
+                return;
+            }
+            if (disposing)
+            {
+                // MatrixParallel.ReturnParallelResources(ref _threads);
+                _countdownEvent.Dispose();
+            }
+            _disposed = true;
+        }
+
+    }
+}

src/Nitload.Mathematics/Numeric/Runtime/QueueFactory.cs

@@ -0,0 +1,21 @@
+// -----------------------------------------------------------------------------
+// file="QueueFactory"
+// Project: StgSharp
+// Copyright (c) Nitload.
+// SPDX-License-Identifier: MIT
+// -----------------------------------------------------------------------------
+
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Nitload.Mathematics.Numeric
+{
+    internal interface IMatrixQueueFactory
+    {
+
+        // IMatrixParallelQueue<IMatrixParallelQueue> CreateSequential();
+    }
+}