feat(types): add i8/i16/u8 element types, storage only

CHANGELOG.md

@@ -15,6 +15,57 @@ against, and semantic-release manages only the **patch** component.
 
 ### Added
 
+- **Narrow integer element types: `i8`, `i16`, `u8`**
+  ([#88](https://github.com/stillwater-sc/mtl5-python/issues/88) phase 1).
+  These are the operand types of MTL5 v5.11.0's widening dot and integer GEMM —
+  the accumulator stays `int32`, which is what `vpdpbusd` / `vpmaddwd` / `SDOT`
+  take — so registering the containers is what makes those kernels reachable at
+  all. **Those kernels arrive with the v5.11.0 pin; this phase does not require
+  it** and is built and tested against the current v5.10.0.
+  `mtl5.vector(np.arange(4, dtype=np.int8))` now gives a
+  `DenseVector_i8`, zero-copy, alongside `DenseMatrix_i8` and the `u8`/`i16`
+  equivalents.
+
+  **This also fixes a silent wrong-dtype bug.** The native factories take
+  `nb::ndarray<T>` without `.noconvert()`, so nanobind's converting second pass
+  handed an int8 array to the *float* overload — registered first — and
+  `mtl5.vector(int8_array)` returned a `DenseVector_f32` that reported
+  `is_view=True` while being a view of the converted temporary, so writes
+  through the NumPy array were invisible. An exact match now resolves in the
+  first pass, before conversion is considered. The dtypes that remain
+  unregistered (`f16`, `uint16`, `uint32`, `uint64`) still behave the old way.
+
+  **`dot`, `norm` and `__matmul__` are deliberately not registered for these
+  types**, and that is a measurement rather than caution. All three default to
+  accumulating in the element type, so on 8- and 16-bit operands they overflow
+  almost immediately. One input covers all three, since each sums the same six
+  products — vectors of six 100s, and a 6×6 of 100s. Exact `dot` and exact
+  `A @ A` element are both 60000; exact `two_norm` is 244.9490:
+
+  | | `i8` | `i16` | `u8` |
+  |---|---|---|---|
+  | `dot` | 96 | −5536 | 96 |
+  | `two_norm` | 9.798 | **nan** | 9.798 |
+  | `A @ A` | 96 | −5536 | 96 |
+
+  `two_norm` is `sqrt(dot)`, which is what makes the `i16` nan legible: the sum
+  wrapped negative, and `sqrt(-5536)` has no answer to give.
+
+  These are MTL5's documented two's-complement wrapping and become correct the
+  moment an int32 accumulator is supplied. `i16` has real headroom — a 2×2 of
+  100s gives the exact 20000 — so its failure needs a longer `k` rather than
+  being immediate; that is a difference of degree, since `k=4` already wraps it
+  and nothing in the API tells a caller where the edge is. `i32` has the same
+  failure past ~46341, where an 8-bit sum of squares overflows at two elements
+  of 12 — an edge case there, essentially every input here.
+
+  `mtl5.dot(i8_vec, i8_vec)` and `i8_mat @ i8_mat` therefore raise `TypeError`
+  rather than returning a wrapped number, until phase 2 lands the accumulator.
+  `__matmul__` is registered on the matrix class itself, so this needed
+  `register_native_matrix` split from `register_native_matrix_matmul` — the
+  existing types are registered exactly as before.
+
+
 - **nanobind 3 is supported**, and the build requirement widens from
   `nanobind>=2.0,<3` to `>=2.0,<4`
   ([#85](https://github.com/stillwater-sc/mtl5-python/issues/85)). The cap added

mtl5/__init__.py

@@ -70,6 +70,8 @@
         # Universal cfloat types
         DenseMatrix_fp8,
         DenseMatrix_fp16,
+        DenseMatrix_i8,
+        DenseMatrix_i16,
         DenseMatrix_i32,
         DenseMatrix_i64,
         # Universal lns types
@@ -84,6 +86,7 @@
         DenseMatrix_qd_cascade,
         DenseMatrix_takum32,
         DenseMatrix_td_cascade,
+        DenseMatrix_u8,
         DenseVector_c64,
         DenseVector_c128,
         DenseVector_cfloat32,
@@ -94,6 +97,8 @@
         DenseVector_fixpnt16,
         DenseVector_fp8,
         DenseVector_fp16,
+        DenseVector_i8,
+        DenseVector_i16,
         DenseVector_i32,
         DenseVector_i64,
         DenseVector_lns16,
@@ -105,6 +110,7 @@
         DenseVector_qd_cascade,
         DenseVector_takum32,
         DenseVector_td_cascade,
+        DenseVector_u8,
         # QR / LQ factorization objects — complex only; the real ones are reached
         # through mtl5.qr()/mtl5.lq(), which dispatch on dtype.
         LQFactor_c64,
@@ -586,8 +592,11 @@ def wrapper(*args, dtype: str = "f64", **kwargs):
     "DenseVector_c128",
     "DenseVector_f32",
     "DenseVector_f64",
+    "DenseVector_i8",
+    "DenseVector_i16",
     "DenseVector_i32",
     "DenseVector_i64",
+    "DenseVector_u8",
     # Typed vector classes — Universal cfloat
     "DenseVector_fp8",
     "DenseVector_fp16",
@@ -613,8 +622,11 @@ def wrapper(*args, dtype: str = "f64", **kwargs):
     "DenseMatrix_c128",
     "DenseMatrix_f32",
     "DenseMatrix_f64",
+    "DenseMatrix_i8",
+    "DenseMatrix_i16",
     "DenseMatrix_i32",
     "DenseMatrix_i64",
+    "DenseMatrix_u8",
     # Typed matrix classes — Universal cfloat
     "DenseMatrix_fp8",
     "DenseMatrix_fp16",

python/src/mtl5_module.cpp

@@ -133,13 +133,17 @@ void register_native_vector(nb::module_& m) {
         });
 }
 
+/// The DenseMatrix_<T> class: storage, indexing, conversion. NO arithmetic --
+/// `__matmul__` is registered separately by register_native_matrix_matmul, so a
+/// type can carry data without also advertising a product it cannot compute
+/// correctly. Returns the class so the caller can add that back.
 template <typename T>
     requires std::is_arithmetic_v<T>
-void register_native_matrix(nb::module_& m) {
+nb::class_<MatrixView<T>> register_native_matrix(nb::module_& m) {
     using MV = MatrixView<T>;
     std::string name = std::string("DenseMatrix_") + type_suffix<T>();
 
-    nb::class_<MV>(m, name.c_str())
+    return nb::class_<MV>(m, name.c_str())
         .def_prop_ro("num_rows", [](const MV& mv) { return mv.mat.num_rows(); })
         .def_prop_ro("num_cols", [](const MV& mv) { return mv.mat.num_cols(); })
         .def_prop_ro("shape", [](const MV& mv) {
@@ -178,40 +182,6 @@ void register_native_matrix(nb::module_& m) {
             }
             return MV(std::move(AT));
         })
-        .def("__matmul__", [](const MV& A, const MV& B) {
-            if (A.mat.num_cols() != B.mat.num_rows())
-                throw std::invalid_argument("matmul: A.num_cols != B.num_rows");
-            mtl::mat::dense2D<T> C(A.mat.num_rows(), B.mat.num_cols());
-            {
-                nogil guard;
-                mtl::mat::dense2D<T> Ac(A.mat.num_rows(), A.mat.num_cols());
-                for (std::size_t i = 0; i < A.mat.num_rows(); ++i)
-                    for (std::size_t j = 0; j < A.mat.num_cols(); ++j)
-                        Ac(i, j) = A.mat(i, j);
-                mtl::mat::dense2D<T> Bc(B.mat.num_rows(), B.mat.num_cols());
-                for (std::size_t i = 0; i < B.mat.num_rows(); ++i)
-                    for (std::size_t j = 0; j < B.mat.num_cols(); ++j)
-                        Bc(i, j) = B.mat(i, j);
-                mtl::mult(Ac, Bc, C);
-            }
-            return MV(std::move(C));
-        })
-        .def("__matmul__", [](const MV& A, const VectorView<T>& x) {
-            if (A.mat.num_cols() != x.vec.size())
-                throw std::invalid_argument("matmul: A.num_cols != len(x)");
-            mtl::vec::dense_vector<T> y(A.mat.num_rows());
-            {
-                nogil guard;
-                mtl::mat::dense2D<T> Ac(A.mat.num_rows(), A.mat.num_cols());
-                for (std::size_t i = 0; i < A.mat.num_rows(); ++i)
-                    for (std::size_t j = 0; j < A.mat.num_cols(); ++j)
-                        Ac(i, j) = A.mat(i, j);
-                mtl::vec::dense_vector<T> xc(x.vec.size());
-                for (std::size_t i = 0; i < x.vec.size(); ++i) xc[i] = x.vec[i];
-                mtl::mult(Ac, xc, y);
-            }
-            return VectorView<T>(std::move(y));
-        })
         .def("copy", [](const MV& mv) {
             auto owned = mtl::mat::dense2D<T>(mv.mat.num_rows(), mv.mat.num_cols());
             {
@@ -274,6 +244,52 @@ void register_native_vector_factory(nb::module_& m) {
 // ---------------------------------------------------------------------------
 // Zero-copy matrix() factory — creates a view borrowing NumPy memory
 // ---------------------------------------------------------------------------
+/// The matrix products. Split out of register_native_matrix because they
+/// accumulate in the ELEMENT type: on an 8- or 16-bit operand that overflows
+/// almost immediately (a 2x2 of 100s gives 32 where the answer is 20000), so the
+/// narrow integer types get the container without it. mtl5.mixed is where an
+/// int32 accumulator will make these meaningful for them.
+template <typename T>
+    requires std::is_arithmetic_v<T>
+void register_native_matrix_matmul(nb::class_<MatrixView<T>>& cls) {
+    using MV = MatrixView<T>;
+    cls
+        .def("__matmul__", [](const MV& A, const MV& B) {
+            if (A.mat.num_cols() != B.mat.num_rows())
+                throw std::invalid_argument("matmul: A.num_cols != B.num_rows");
+            mtl::mat::dense2D<T> C(A.mat.num_rows(), B.mat.num_cols());
+            {
+                nogil guard;
+                mtl::mat::dense2D<T> Ac(A.mat.num_rows(), A.mat.num_cols());
+                for (std::size_t i = 0; i < A.mat.num_rows(); ++i)
+                    for (std::size_t j = 0; j < A.mat.num_cols(); ++j)
+                        Ac(i, j) = A.mat(i, j);
+                mtl::mat::dense2D<T> Bc(B.mat.num_rows(), B.mat.num_cols());
+                for (std::size_t i = 0; i < B.mat.num_rows(); ++i)
+                    for (std::size_t j = 0; j < B.mat.num_cols(); ++j)
+                        Bc(i, j) = B.mat(i, j);
+                mtl::mult(Ac, Bc, C);
+            }
+            return MV(std::move(C));
+        })
+        .def("__matmul__", [](const MV& A, const VectorView<T>& x) {
+            if (A.mat.num_cols() != x.vec.size())
+                throw std::invalid_argument("matmul: A.num_cols != len(x)");
+            mtl::vec::dense_vector<T> y(A.mat.num_rows());
+            {
+                nogil guard;
+                mtl::mat::dense2D<T> Ac(A.mat.num_rows(), A.mat.num_cols());
+                for (std::size_t i = 0; i < A.mat.num_rows(); ++i)
+                    for (std::size_t j = 0; j < A.mat.num_cols(); ++j)
+                        Ac(i, j) = A.mat(i, j);
+                mtl::vec::dense_vector<T> xc(x.vec.size());
+                for (std::size_t i = 0; i < x.vec.size(); ++i) xc[i] = x.vec[i];
+                mtl::mult(Ac, xc, y);
+            }
+            return VectorView<T>(std::move(y));
+        });
+}
+
 template <typename T>
     requires std::is_arithmetic_v<T>
 void register_native_matrix_factory(nb::module_& m) {
@@ -1004,13 +1020,30 @@ void register_universal_solve(nb::module_& m) {
 // ===========================================================================
 // Convenience: register all bindings for one type
 // ===========================================================================
+/// Containers and the zero-copy NumPy factories -- storage, no arithmetic.
+///
+/// Split out from register_native because the narrow integer types can carry
+/// data long before they can compute with it: their operations need an
+/// accumulator wider than the element type, and until that exists the generic
+/// element-typed forms are not merely inexact but wrong. See the registration
+/// site for the measurements.
 template <typename T>
     requires std::is_arithmetic_v<T>
-void register_native(nb::module_& m) {
+nb::class_<MatrixView<T>> register_native_storage(nb::module_& m) {
     register_native_vector<T>(m);
-    register_native_matrix<T>(m);
+    auto mat_cls = register_native_matrix<T>(m);
     register_native_vector_factory<T>(m);
     register_native_matrix_factory<T>(m);
+    // Returned rather than discarded so register_native can add `__matmul__`
+    // back. The vector class needs no such treatment: it carries no arithmetic.
+    return mat_cls;
+}
+
+template <typename T>
+    requires std::is_arithmetic_v<T>
+void register_native(nb::module_& m) {
+    auto mat_cls = register_native_storage<T>(m);
+    register_native_matrix_matmul<T>(mat_cls);
     register_native_norm_overload<T>(m);
     register_native_dot_overload<T>(m);
 }
@@ -1599,6 +1632,58 @@ NB_MODULE(_core, m) {
     register_native<int32_t>(m);              // i32
     register_native<int64_t>(m);              // i64
 
+    // Narrow integers: STORAGE ONLY, deliberately. These are the OPERAND types
+    // of MTL5's widening dot and integer GEMM -- the accumulator stays int32,
+    // which is what vpdpbusd / vpmaddwd / SDOT take. Registering the containers
+    // is what makes those kernels reachable at all; the accumulator policy that
+    // drives them is a separate step (mtl5.mixed).
+    //
+    // `norm`, `dot` and `__matmul__` are NOT registered for these, and that is a
+    // measurement rather than caution. All three default to accumulating in the
+    // ELEMENT type, so on 8- and 16-bit operands they overflow almost
+    // immediately. One input covers all three, because each sums the same six
+    // products: vectors of six 100s, and a 6x6 of 100s. Exact dot and exact
+    // A@A element are both 60000; exact two_norm is 244.9490.
+    //
+    //     op         i8        i16        u8
+    //     dot        96        -5536      96
+    //     two_norm   9.798     nan        9.798
+    //     A @ A      96        -5536      96
+    //
+    // two_norm is sqrt(dot), which is what makes the i16 nan legible: the sum
+    // wrapped NEGATIVE, and sqrt(-5536) has no answer to give.
+    //
+    // i16 has real headroom -- a 2x2 of 100s gives the exact 20000 -- so its
+    // failure needs a longer k rather than being immediate. That is a difference
+    // of degree, not of kind: k=4 already wraps it to -25536, and nothing in the
+    // API tells a caller where the edge is. The accumulator is what moves it.
+    //
+    // The dot results are MTL5's documented two's-complement wrapping and become
+    // useful the moment an int32 accumulator is supplied. `two_norm` is worse
+    // than wrapping: it takes sqrt of a sum that has wrapped, and of a NEGATIVE
+    // one for i16, so it yields nan. i32 has the same failure but only past
+    // ~46341, where an 8-bit sum of squares overflows at two elements of 12 --
+    // an edge case there, essentially every input here.
+    //
+    // Exposing an operation that is wrong for almost all inputs is worse than
+    // not exposing it, so these three carry data and nothing else until the
+    // accumulator lands. `mtl5.dot(i8_vec, i8_vec)` and `i8_mat @ i8_mat` raise
+    // TypeError rather than returning a wrapped number.
+    //
+    // Registering the containers also fixes a silent wrong-dtype bug. These
+    // factories take nb::ndarray<T> WITHOUT .noconvert(), so nanobind's second
+    // (converting) pass used to hand an int8 array to the float overload --
+    // registered first -- and `mtl5.vector(np.arange(4, dtype=np.int8))`
+    // returned a DenseVector_f32 that reported is_view=True while being a view
+    // of the converted temporary, so writes through the NumPy array were
+    // invisible. An exact match now resolves in the FIRST pass, before
+    // conversion is considered. Every dtype still unregistered (f16, uint16,
+    // uint32, uint64) keeps that behaviour -- see the note in mtl5_ndarray.cpp,
+    // where .noconvert() was added for exactly this.
+    register_native_storage<int8_t>(m);       // i8
+    register_native_storage<int16_t>(m);      // i16
+    register_native_storage<uint8_t>(m);      // u8
+
     // ----- Sparse matrices (CSR via mtl::compressed2D) -----------------------
     register_sparse_matrix<float>(m);
     register_sparse_matrix<double>(m);

python/src/mtl5_types.hpp

@@ -119,6 +119,14 @@ template <> constexpr const char* type_suffix<float>()   { return "f32"; }
 template <> constexpr const char* type_suffix<double>()  { return "f64"; }
 template <> constexpr const char* type_suffix<int32_t>() { return "i32"; }
 template <> constexpr const char* type_suffix<int64_t>() { return "i64"; }
+// Narrow integers, for MTL5's widening dot and integer GEMM: the accumulator is
+// int32 while the OPERANDS are 8- or 16-bit, which is what the hardware
+// multiply-accumulate instructions take. int8_t/uint8_t are signed/unsigned
+// char, so these are distinct specializations from the ones above rather than
+// aliases of them.
+template <> constexpr const char* type_suffix<int8_t>()  { return "i8"; }
+template <> constexpr const char* type_suffix<int16_t>() { return "i16"; }
+template <> constexpr const char* type_suffix<uint8_t>() { return "u8"; }
 template <> constexpr const char* type_suffix<fp8>()     { return "fp8"; }
 template <> constexpr const char* type_suffix<fp16>()    { return "fp16"; }
 template <> constexpr const char* type_suffix<posit8>()  { return "posit8"; }