Project 5: Const Generic Matrix - Type-Level Linear Algebra

Build a Matrix library where dimensions (rows/columns) are encoded in the type system, enabling compile-time dimension checking for all matrix operations. This eliminates an entire class of runtime errors by making invalid matrix operations impossible to express.

Quick Reference

Attribute Value
Difficulty Level 2: Intermediate
Time Estimate 1-2 weeks
Language Rust
Prerequisites Basic Rust (ownership, traits, generics), elementary linear algebra
Key Topics Const generics, monomorphization, operator overloading, zero-cost abstractions

1. Learning Objectives

By completing this project, you will:

  1. Understand const generics - How to use constant values as type parameters (const N: usize)
  2. Master monomorphization - How the compiler generates specialized code for each matrix size
  3. Implement operator traits - Overload +, -, * for custom types with proper trait bounds
  4. Apply type-level constraints - Enforce mathematical rules (M x N * N x P = M x P) at compile time
  5. Design zero-cost abstractions - Prove through benchmarks that type safety has no runtime overhead
  6. Use fixed-size arrays - Understand why [[T; C]; R] beats Vec<Vec<T>> for known dimensions
  7. Explore nightly features - Learn about generic_const_exprs for computed const values
  8. Build production-quality APIs - Create ergonomic interfaces with new, zeros, identity, transpose
  9. Write comprehensive tests - Including compile-fail tests that verify type errors
  10. Connect to real libraries - Understand how nalgebra, ndarray, and graphics engines use these patterns

2. Theoretical Foundation

2.1 Core Concepts

Type-Level Programming: Encoding Values in Types

In most programming, values exist at runtime:

// Runtime value - known when program runs
let rows: usize = 3;
let cols: usize = 4;

In type-level programming, values exist at compile time:

// Compile-time value - known when program is compiled
struct Matrix<const R: usize, const C: usize>;
let m: Matrix<3, 4>;  // The "3" and "4" ARE the type

This shift is profound. When dimensions are part of the type:

Regular Matrix:            Const Generic Matrix:
+------------------+       +------------------+
| rows: usize      |       | NO RUNTIME DATA  |
| cols: usize      |       | Type IS the data |
| data: Vec<f32>   |       | data: [[f32;C];R]|
+------------------+       +------------------+
     |                            |
     v                            v
  Runtime check:              Compile-time check:
  if rows != other.cols       The code won't compile
     panic!("Mismatch!")      if dimensions mismatch

Monomorphization: How the Compiler Specializes Your Code

When you write generic code, Rust doesn’t create one “generic” function. It creates a specialized copy for each concrete type used:

// You write this:
impl<const R: usize, const C: usize> Matrix<R, C> {
    pub fn zeros() -> Self {
        Self { data: [[0.0; C]; R] }
    }
}

// The compiler generates these (conceptually):
impl Matrix<2, 2> {
    pub fn zeros() -> Self {
        Self { data: [[0.0, 0.0], [0.0, 0.0]] }
    }
}

impl Matrix<3, 4> {
    pub fn zeros() -> Self {
        Self { data: [[0.0, 0.0, 0.0, 0.0],
                      [0.0, 0.0, 0.0, 0.0],
                      [0.0, 0.0, 0.0, 0.0]] }
    }
}

Source Code                      Compiled Binary
+------------------+             +------------------+
| Matrix<R, C>     |             | Matrix_2_2       |
| generic code     | ==========> | Matrix_3_4       |
|                  |   rustc     | Matrix_10_10     |
+------------------+             +------------------+
                                 (specialized copies)

Tradeoff: More sizes = larger binary. But each function is optimized perfectly for its specific size.

2.2 Why This Matters

Compile-time dimension checking eliminates entire categories of bugs:

RUNTIME CHECKING (Traditional)

fn multiply(a: DynMatrix, b: DynMatrix) -> Result<DynMatrix, Error> {
    if a.cols != b.rows {           // Check happens at RUNTIME
        return Err("Dimension mismatch!");
    }
    // ... perform multiplication
}

COMPILE-TIME CHECKING (Const Generics)

fn multiply<const R1: usize, const C1: usize, const C2: usize>
    (a: Matrix<R1, C1>, b: Matrix<C1, C2>) -> Matrix<R1, C2>
{                            ^^^      ^^^
    // No check needed!      These MUST match or code won't compile
}

Real-world impact:

  • NASA Mars Climate Orbiter: Lost due to unit conversion error (cost: $327 million)
  • Graphics programming: Shader matrix mismatches cause visual artifacts or crashes
  • Machine learning: Tensor shape errors are the #1 debugging pain point

2.3 Historical Context

Const generics were one of Rust’s most requested features:

2017: RFC 2000 proposed (Const Generics)
      "Allow generic parameters over constant values"

2019: Partial implementation lands (min_const_generics)
      Limited to integers, no const expressions

2021: Rust 1.51 stabilizes min_const_generics
      - Arrays can use const generics: [T; N]
      - Structs can have const parameters
      - Basic arithmetic NOT allowed yet

2023+: generic_const_exprs (nightly)
       - Allows R * C computations
       - Enables more complex type math

Before const generics, you had two bad options:

// Option 1: Macros (ugly, slow to compile)
macro_rules! matrix {
    ($rows:expr, $cols:expr) => {
        struct Matrix { data: [[f32; $cols]; $rows] }
    };
}

// Option 2: Runtime dimensions (slow, error-prone)
struct Matrix {
    rows: usize,
    cols: usize,
    data: Vec<f32>,
}

2.4 Common Misconceptions

Misconception 1: “Const generics add runtime overhead”

  • Reality: They add ZERO runtime overhead. The type parameters exist only at compile time.

Misconception 2: “I can do arithmetic with const generics on stable Rust”

  • Reality: Complex const expressions like R * C require the nightly generic_const_exprs feature.

Misconception 3: “Type State pattern is the same as const generics”

  • Reality: Type State uses phantom type parameters; const generics use actual values in types.

Misconception 4: “More matrix sizes = slower code”

  • Reality: Each specialized version is perfectly optimized. The code is FASTER, but binary is larger.

3. Project Specification

3.1 What You Will Build

A Matrix library where the dimensions (rows/cols) are part of the type. Multiplying a Matrix<3, 2> by a Matrix<2, 4> is allowed, but multiplying by a Matrix<5, 5> will fail to compile.

Example Code:

let a = Matrix::<3, 2>::new([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]]);
let b = Matrix::<2, 4>::new([[1.0, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0]]);
let c = a * b; // Compiles! Result is Matrix<3, 4>

let d = Matrix::<5, 5>::new();
let e = a * d; // Error: "expected Matrix<2, _>, found Matrix<5, 5>"

3.2 Functional Requirements

  • Define Matrix<const R: usize, const C: usize> struct
  • Implement new(), zeros(), ones() constructors
  • Implement Debug trait for pretty-printing
  • Implement Add trait for same-dimension matrices
  • Implement Sub trait for same-dimension matrices
  • Implement Mul trait with dimension checking (cols of A = rows of B)
  • Implement scalar multiplication
  • Implement transpose() returning Matrix<C, R>
  • Implement identity() for square matrices only
  • Implement Index and IndexMut for element access
  • Write compile-fail tests proving invalid operations don’t compile

3.3 Non-Functional Requirements

  • Zero runtime overhead: Prove via benchmarks that performance matches hand-written arrays
  • No heap allocation: Use fixed-size arrays, not Vec
  • Ergonomic API: Support natural syntax like a + b and a * b
  • Comprehensive documentation: Every public item documented

3.4 Example Usage / Output

$ cargo new --lib const-matrix
     Created library `const-matrix` package

$ cd const-matrix

$ cargo test
   Compiling const-matrix v0.1.0
    Finished test [unoptimized + debuginfo] target(s) in 1.23s
     Running unittests src/lib.rs

running 8 tests
test tests::test_zeros ... ok
test tests::test_ones ... ok
test tests::test_addition ... ok
test tests::test_subtraction ... ok
test tests::test_scalar_mul ... ok
test tests::test_matrix_mul_dimensions ... ok
test tests::test_transpose ... ok
test tests::test_identity ... ok

test result: ok. 8 passed; 0 failed; 0 ignored

$ cargo run --example demo
   Compiling const-matrix v0.1.0
    Finished dev [unoptimized + debuginfo] target(s) in 0.45s
     Running `target/debug/examples/demo`

=== Const Generic Matrix Demo ===

Creating a 3x2 matrix:
Matrix<3, 2> [
  [  1.0000,   2.0000]
  [  3.0000,   4.0000]
  [  5.0000,   6.0000]
]

Creating a 2x4 matrix:
Matrix<2, 4> [
  [  1.0000,   2.0000,   3.0000,   4.0000]
  [  5.0000,   6.0000,   7.0000,   8.0000]
]

Multiplying 3x2 * 2x4 = 3x4:
Matrix<3, 4> [
  [ 11.0000,  14.0000,  17.0000,  20.0000]
  [ 23.0000,  30.0000,  37.0000,  44.0000]
  [ 35.0000,  46.0000,  57.0000,  68.0000]
]

Dimension mismatch would be caught at compile time!
// let bad = matrix_3x2 * matrix_3x4;  // Won't compile!

4. Solution Architecture

4.1 High-Level Design

+-----------------------------------------------------------+
|                     Matrix<R, C>                          |
+-----------------------------------------------------------+
|  Compile-Time Constants:                                  |
|    - R: Number of rows (const generic)                    |
|    - C: Number of columns (const generic)                 |
+-----------------------------------------------------------+
|  Runtime Data:                                            |
|    - data: [[f64; C]; R] (fixed-size nested array)        |
+-----------------------------------------------------------+
|  Operations (Trait Implementations):                      |
|    - Add<Matrix<R,C>> -> Matrix<R,C>                      |
|    - Sub<Matrix<R,C>> -> Matrix<R,C>                      |
|    - Mul<Matrix<C,P>> -> Matrix<R,P>   <-- KEY!           |
|    - Mul<f64> -> Matrix<R,C> (scalar)                     |
+-----------------------------------------------------------+

4.2 Key Components

  1. Core Struct: Matrix<const R: usize, const C: usize>
  2. Storage: [[f64; C]; R] - row-major nested array
  3. Operator Traits: Add, Sub, Mul, Index, IndexMut
  4. Utility Methods: new(), zeros(), ones(), identity(), transpose()
  5. Display: Debug implementation for pretty-printing

4.3 Data Structures

/// A fixed-size matrix with dimensions known at compile time.
#[derive(Clone, Copy, PartialEq)]
pub struct Matrix<const R: usize, const C: usize> {
    data: [[f64; C]; R],
}

// Memory layout for Matrix<3, 4>:
// +------+------+------+------+
// | a00  | a01  | a02  | a03  |  <- Row 0
// +------+------+------+------+
// | a10  | a11  | a12  | a13  |  <- Row 1
// +------+------+------+------+
// | a20  | a21  | a22  | a23  |  <- Row 2
// +------+------+------+------+
// Contiguous in memory: [a00, a01, a02, a03, a10, a11, ...]

4.4 Algorithm Overview

Matrix Multiplication Algorithm:

For Matrix<R1, C1> * Matrix<C1, C2> = Matrix<R1, C2>:

result[i][j] = sum(left[i][k] * right[k][j]) for k in 0..C1

Example: Matrix<2,3> * Matrix<3,2> = Matrix<2,2>

    Left (2x3)          Right (3x2)         Result (2x2)
    [a b c]             [x u]               [ax+bz+cy  au+bw+cv]
    [d e f]      *      [z w]       =       [dx+ez+fy  du+ew+fv]
                        [y v]

5. Implementation Guide

5.1 Development Environment Setup

# Create project
cargo new --lib const-matrix
cd const-matrix

# Add dev dependencies for testing
cargo add --dev criterion proptest trybuild

Cargo.toml:

[package]
name = "const-matrix"
version = "0.1.0"
edition = "2021"

[dev-dependencies]
criterion = "0.5"
proptest = "1.4"
trybuild = "1.0"

[[bench]]
name = "matrix_bench"
harness = false

5.2 Project Structure

const-matrix/
├── Cargo.toml
├── src/
│   └── lib.rs              # Main library code
├── examples/
│   └── demo.rs             # Usage demonstration
├── benches/
│   └── matrix_bench.rs     # Performance benchmarks
└── tests/
    ├── integration.rs      # Integration tests
    └── ui/                 # Compile-fail tests
        ├── add_mismatch.rs
        └── mul_mismatch.rs

5.3 The Core Question You’re Answering

“Can I force the compiler to understand linear algebra?”

Yes. By moving values (like 3 and 2) into the type system, the compiler’s trait solver becomes a dimensional analysis engine. Every matrix operation is verified before your program runs.

5.4 Concepts You Must Understand First

Stop and research these before coding:

  1. Const Generics basics
    • How to define struct Matrix<const R: usize, const C: usize>
    • Book Reference: “Programming Rust” Ch. 11
  2. Monomorphization
    • Why does the compiler generate a new version of your function for every different size?
    • Book Reference: “Idiomatic Rust” Ch. 5
  3. Trait implementation for Generics
    • How to implement Mul only for matrices where columns of A match rows of B

5.5 Questions to Guide Your Design

  1. Storage
    • Should you use Vec<T> or [T; R * C]? (Hint: Since R and C are const, an array is faster!)
  2. Operations
    • How do you define the return type of a multiplication as Matrix<R1, C2>?
  3. Bounds
    • How do you handle cases where you need R * C to be calculated at compile time? (Hint: generic_const_exprs feature on Nightly)

5.6 Thinking Exercise

The Cost of Freedom

If you have 100 different matrix sizes in your program, how does that affect the binary size? Compare this to a library where dimensions are stored as runtime integers.

Answer: Each unique Matrix<R, C> type generates its own code. 100 sizes means 100 copies of each method. This is the “monomorphization cost” - larger binaries but perfect optimization for each size. Runtime-sized matrices have one code copy but add branching and bounds checks.

5.7 Hints in Layers

Hint 1: The definition

struct Matrix<const R: usize, const C: usize> { data: [[f32; C]; R] }

Hint 2: Implementation Use the impl block:

impl<const R: usize, const C: usize> Matrix<R, C> { ... }

Hint 3: Multiplication The Mul trait implementation:

impl<const R1: usize, const C1: usize, const C2: usize>
    Mul<Matrix<C1, C2>> for Matrix<R1, C1>
{
    type Output = Matrix<R1, C2>;
    ...
}

Hint 4: Verification Use compile-fail tests with trybuild crate to verify invalid operations don’t compile.

5.8 The Interview Questions They’ll Ask

  1. “What is the difference between a generic type parameter and a const generic parameter?”

  2. “How do const generics reduce runtime overhead?”

  3. “What are the limitations of const generics in Stable Rust currently?”

  4. “Can you explain ‘Type State’ pattern using const generics?”

  5. “How does monomorphization affect binary size vs runtime performance?”

5.9 Books That Will Help

Topic Book Chapter
Const Generics “Programming Rust” Ch. 11
Matrix Math “Math for Programmers” Ch. 4
Type-level logic “Idiomatic Rust” Ch. 5
Zero-Cost Abstractions “Rust for Rustaceans” Ch. 1

5.10 Implementation Phases

Phase 1: Foundation (Day 1)

Goal: Create the struct and basic constructors.

#[derive(Clone, Copy, PartialEq)]
pub struct Matrix<const R: usize, const C: usize> {
    data: [[f64; C]; R],
}

impl<const R: usize, const C: usize> Matrix<R, C> {
    pub fn new(data: [[f64; C]; R]) -> Self {
        Self { data }
    }

    pub fn zeros() -> Self {
        Self { data: [[0.0; C]; R] }
    }

    pub fn ones() -> Self {
        Self { data: [[1.0; C]; R] }
    }
}

Checkpoint: Can you create Matrix::<3, 4>::zeros()?

Phase 2: Debug and Access (Day 2-3)

Goal: Make matrices printable and add element access.

impl<const R: usize, const C: usize> std::fmt::Debug for Matrix<R, C> {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        writeln!(f, "Matrix<{}, {}> [", R, C)?;
        for row in &self.data {
            write!(f, "  [")?;
            for (i, val) in row.iter().enumerate() {
                if i > 0 { write!(f, ", ")?; }
                write!(f, "{:8.4}", val)?;
            }
            writeln!(f, "]")?;
        }
        write!(f, "]")
    }
}

Phase 3: Addition and Subtraction (Day 4-5)

Goal: Enable matrix_a + matrix_b syntax.

impl<const R: usize, const C: usize> std::ops::Add for Matrix<R, C> {
    type Output = Self;

    fn add(self, rhs: Self) -> Self::Output {
        let mut result = [[0.0; C]; R];
        for i in 0..R {
            for j in 0..C {
                result[i][j] = self.data[i][j] + rhs.data[i][j];
            }
        }
        Matrix { data: result }
    }
}

Phase 4: Multiplication with Dimension Checking (Day 6-8)

Goal: Matrix multiplication with compile-time dimension verification.

impl<const R1: usize, const C1: usize, const C2: usize>
    std::ops::Mul<Matrix<C1, C2>> for Matrix<R1, C1>
{
    type Output = Matrix<R1, C2>;

    fn mul(self, rhs: Matrix<C1, C2>) -> Self::Output {
        let mut result = [[0.0; C2]; R1];
        for i in 0..R1 {
            for j in 0..C2 {
                for k in 0..C1 {
                    result[i][j] += self.data[i][k] * rhs.data[k][j];
                }
            }
        }
        Matrix { data: result }
    }
}

5.11 Key Implementation Decisions

Decision Options Recommended Rationale
Element type Generic T vs fixed f64 Start with f64, generalize later Simpler bounds
Storage [[T;C];R] vs [T;R*C] [[T;C];R] Natural row access
Copy semantics Copy vs move Copy for small matrices Ergonomic
Row vs column major Row-major vs column-major Row-major C convention

6. Testing Strategy

6.1 Unit Tests

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_zeros() {
        let m: Matrix<3, 4> = Matrix::zeros();
        assert_eq!(m.data, [[0.0; 4]; 3]);
    }

    #[test]
    fn test_addition() {
        let a = Matrix::new([[1.0, 2.0], [3.0, 4.0]]);
        let b = Matrix::new([[5.0, 6.0], [7.0, 8.0]]);
        let c = a + b;
        assert_eq!(c.data, [[6.0, 8.0], [10.0, 12.0]]);
    }

    #[test]
    fn test_multiplication_dimensions() {
        let a: Matrix<3, 2> = Matrix::ones();
        let b: Matrix<2, 4> = Matrix::ones();
        let c: Matrix<3, 4> = a * b;

        // Each element is sum of row_a * col_b = 1*1 + 1*1 = 2
        assert_eq!(c.data[0][0], 2.0);
    }
}

6.2 Integration Tests

// tests/integration.rs
use const_matrix::Matrix;

#[test]
fn test_chain_operations() {
    let a: Matrix<2, 3> = Matrix::ones();
    let b: Matrix<3, 2> = Matrix::ones();
    let c: Matrix<2, 2> = Matrix::identity();

    let result = (a * b) * c;
    // (2x3) * (3x2) = (2x2), then (2x2) * (2x2) = (2x2)
    assert_eq!(result.data[0][0], 3.0);
}

6.3 Property-Based Testing

use proptest::prelude::*;

proptest! {
    #[test]
    fn matrix_addition_commutative(
        a00 in -1000.0..1000.0f64,
        a01 in -1000.0..1000.0f64,
        b00 in -1000.0..1000.0f64,
        b01 in -1000.0..1000.0f64,
    ) {
        let a = Matrix::new([[a00, a01]]);
        let b = Matrix::new([[b00, b01]]);
        prop_assert_eq!(a + b, b + a);
    }
}

6.4 Compile-Fail Tests

Create tests that verify invalid operations don’t compile:

// tests/ui/mul_mismatch.rs
use const_matrix::Matrix;

fn main() {
    let a: Matrix<3, 2> = Matrix::zeros();
    let b: Matrix<5, 4> = Matrix::zeros();
    let _ = a * b;  // ERROR: no impl for Mul<Matrix<5, 4>>
}

// tests/compile_fail.rs
#[test]
fn compile_fail_tests() {
    let t = trybuild::TestCases::new();
    t.compile_fail("tests/ui/*.rs");
}

7. Common Pitfalls & Debugging

Problem Symptom Root Cause Fix
Type inference failure “cannot infer R and C” Missing type annotation Use Matrix::<3,4>::zeros() or let m: Matrix<3,4> = ...
Stack overflow Crash on large matrix Too big for stack Box large matrices: Box::new(Matrix::zeros())
Confusing row/col Wrong results Index order confusion Remember: data[row][col]
Wrong output type Type mismatch in multiply Forgot output is Matrix<R1, C2> Check multiplication signature
Copy trait missing Cannot use matrix twice Forgot #[derive(Copy)] Add Copy bound or clone
Bounds check overhead Slower than expected Using get() in hot loop Use direct indexing after validation

Debugging Checklist

  1. Compilation errors about trait bounds: Check that your const generic parameters match exactly in trait impls
  2. “the trait bound is not satisfied”: Ensure your impl covers the right types
  3. Infinite monomorphization: Avoid recursive types with growing const params
  4. Stack overflow: Large matrices should use Box<Matrix<R,C>>

8. Extensions & Challenges

Extension 1: Gaussian Elimination

impl<const N: usize> Matrix<N, N> {
    /// Solve Ax = b using Gaussian elimination
    pub fn solve(&self, b: Matrix<N, 1>) -> Option<Matrix<N, 1>> {
        // Implementation: augment [A|b], forward elimination, back substitution
    }

    /// Calculate determinant
    pub fn determinant(&self) -> f64 {
        // Use row reduction, det = product of pivots
    }

    /// Calculate inverse (if exists)
    pub fn inverse(&self) -> Option<Self> {
        // Use Gauss-Jordan with identity augmented
    }
}

Extension 2: Generic Element Type

use num_traits::{Zero, One, NumOps};

pub struct Matrix<T, const R: usize, const C: usize>
where
    T: Copy + Zero + One + NumOps,
{
    data: [[T; C]; R],
}

// Now works with f32, f64, i32, complex numbers, etc.

Extension 3: SIMD-Optimized Operations

#![feature(portable_simd)]

impl<const R: usize, const C: usize> Matrix<R, C>
where
    [(); C % 4]:,  // Ensure C is multiple of 4
{
    pub fn add_simd(self, rhs: Self) -> Self {
        // Use SIMD for 4 f64s at once
    }
}

Extension 4: Strassen Algorithm

Implement Strassen’s matrix multiplication for O(n^2.807) complexity instead of O(n^3).

9. Real-World Connections

nalgebra

The industry-standard linear algebra library for Rust uses const generics:

use nalgebra::{SMatrix, SVector};

// Fixed-size matrix (compile-time dimensions)
type Mat3x3 = SMatrix<f64, 3, 3>;

// Your project teaches the fundamentals that nalgebra builds on

Graphics Programming (glam, cgmath)

Game engines need fast, fixed-size math:

use glam::Mat4;  // 4x4 matrix, size known at compile time

// Transformations in 3D graphics:
let model = Mat4::from_rotation_y(angle);
let view = Mat4::look_at_rh(eye, target, up);
let projection = Mat4::perspective_rh(fov, aspect, near, far);
let mvp = projection * view * model;  // Type-safe multiplication!

Machine Learning (burn, candle)

Neural network frameworks use const generics for tensor shapes:

// In modern ML frameworks:
type LayerWeights = Matrix<784, 128>;  // MNIST input -> hidden layer
type LayerBias = Matrix<128, 1>;

// Shape mismatches caught at compile time!

10. Resources

10.1 Documentation

10.2 Crates

10.3 Articles & Talks

11. Self-Assessment Checklist

  • I can explain what const generics are and why they’re useful
  • I can implement struct Matrix<const R: usize, const C: usize>
  • I understand monomorphization and its tradeoffs
  • I can implement operator traits (Add, Mul) with type constraints
  • I understand why Mul<Matrix<C1, C2>> returns Matrix<R1, C2>
  • I can write compile-fail tests using trybuild
  • I can explain the difference between Type State and const generics
  • My implementation compiles with no warnings
  • All tests pass (including compile-fail tests)
  • I can answer all the interview questions confidently

12. Submission / Completion Criteria

Your project is complete when:

  1. Core Implementation
    • Matrix<R, C> struct with const generic parameters
    • new(), zeros(), ones() constructors
    • Add, Sub, Mul trait implementations
    • identity() for square matrices
    • transpose() returning Matrix<C, R>
  2. Testing
    • Unit tests for all operations
    • Compile-fail tests proving dimension mismatches don’t compile
    • Property-based tests for mathematical properties
  3. Documentation
    • All public items documented
    • Example in examples/demo.rs
    • README with usage instructions
  4. Verification
    • cargo test passes
    • cargo clippy has no warnings
    • cargo doc generates clean documentation

Final Deliverable: A working library that catches matrix dimension errors at compile time, with comprehensive tests proving both correct operations and rejected invalid operations.