From 60af7a91e5878c6d2f68d7e8d094d7073b405f7f Mon Sep 17 00:00:00 2001
From: Sarthak Tayal <sarthaktayal2@gmail.com>
Date: Wed, 5 Aug 2026 16:53:35 -0400
Subject: [PATCH 1/6]

---
 src/matrix.js | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/matrix.js b/src/matrix.js
index 887b59b..74d2c39 100644
--- a/src/matrix.js
+++ b/src/matrix.js
@@ -1230,8 +1230,8 @@ export class AbstractMatrix {
       // Crop output to the desired size (undo dynamic padding).
       let result = AbstractMatrix.zeros(2 * c11.rows, 2 * c11.columns);
       result = result.setSubMatrix(c11, 0, 0);
-      result = result.setSubMatrix(c12, c11.rows, 0);
-      result = result.setSubMatrix(c21, 0, c11.columns);
+      result = result.setSubMatrix(c12, 0, c11.columns);
+      result = result.setSubMatrix(c21, c11.rows, 0);
       result = result.setSubMatrix(c22, c11.rows, c11.columns);
       return result.subMatrix(0, rows - 1, 0, cols - 1);
     }

From a572b2262052ecd12ba16cae13002a22c3058964 Mon Sep 17 00:00:00 2001
From: Sarthak Tayal <sarthaktayal2@gmail.com>
Date: Wed, 5 Aug 2026 16:54:26 -0400
Subject: [PATCH 2/6] fix: mmulStrassen returns the matrix product for every
 input shape

---
 src/matrix.js | 23 ++++++++++++++---------
 1 file changed, 14 insertions(+), 9 deletions(-)

diff --git a/src/matrix.js b/src/matrix.js
index 74d2c39..f06e4be 100644
--- a/src/matrix.js
+++ b/src/matrix.js
@@ -1135,6 +1135,9 @@ export class AbstractMatrix {
         `Multiplying ${r1} x ${c1} and ${r2} x ${c2} matrix: dimensions do not match.`,
       );
     }
+    if (r1 === 0 || c2 === 0) {
+      return new Matrix(r1, c2);
+    }
 
     // Put a matrix into the top left of a matrix of zeros.
     // `rows` and `cols` are the dimensions of the output matrix.
@@ -1150,14 +1153,12 @@ export class AbstractMatrix {
       }
     }
 
-    // Make sure both matrices are the same size.
-    // This is exclusively for simplicity:
-    // this algorithm can be implemented with matrices of different sizes.
-
-    let r = Math.max(r1, r2);
-    let c = Math.max(c1, c2);
-    x = embed(x, r, c);
-    y = embed(y, r, c);
+    // pad both operands into the same square so that the block split lines up.
+    // zeros never reach the top left r1 x c2 corner of the product, which is
+    // the part that gets returned.
+    let n = Math.max(r1, c1, r2, c2);
+    x = embed(x, n, n);
+    y = embed(y, n, n);
 
     // Our recursive multiplication function.
     function blockMult(a, b, rows, cols) {
@@ -1236,7 +1237,11 @@ export class AbstractMatrix {
       return result.subMatrix(0, rows - 1, 0, cols - 1);
     }
 
-    return blockMult(x, y, r, c);
+    const product = blockMult(x, y, n, n);
+    if (product.rows === r1 && product.columns === c2) {
+      return product;
+    }
+    return product.subMatrix(0, r1 - 1, 0, c2 - 1);
   }
 
   scaleRows(options = {}) {

From 8ad71537934602d501babff661016c4de28d4b1f Mon Sep 17 00:00:00 2001
From: Sarthak Tayal <sarthaktayal2@gmail.com>
Date: Wed, 5 Aug 2026 16:55:06 -0400
Subject: [PATCH 3/6]

---
 src/__tests__/matrix/utility.test.js | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/src/__tests__/matrix/utility.test.js b/src/__tests__/matrix/utility.test.js
index b1893a3..52fa6ee 100644
--- a/src/__tests__/matrix/utility.test.js
+++ b/src/__tests__/matrix/utility.test.js
@@ -374,14 +374,10 @@ describe('utility methods', () => {
 
   it('mmul strassen on empty matrices', () => {
     // https://github.com/mljs/matrix/issues/114
-    // while the mathematically correct result is 0x0, we assert a 2x2 padded result that the current implementation produces
-    // (this call is actually just delegated to standard multiplication in mmul())
-    expect(
-      new Matrix(0, 2).mmulStrassen(new Matrix(2, 0)).to2DArray(),
-    ).toStrictEqual([
-      [0, 0],
-      [0, 0],
-    ]);
+    const result = new Matrix(0, 2).mmulStrassen(new Matrix(2, 0));
+    expect(result.rows).toBe(0);
+    expect(result.columns).toBe(0);
+    expect(result.to2DArray()).toStrictEqual([]);
   });
 
   it('mmul 2x2 and 3x3', () => {

From 59bc63f56f2c8ae6eaa0f0d5106dad5319648d2f Mon Sep 17 00:00:00 2001
From: Sarthak Tayal <sarthaktayal2@gmail.com>
Date: Wed, 5 Aug 2026 16:56:10 -0400
Subject: [PATCH 4/6]

---
 src/__tests__/matrix/mmulStrassen.test.js | 125 ++++++++++++++++++++++
 1 file changed, 125 insertions(+)
 create mode 100644 src/__tests__/matrix/mmulStrassen.test.js

diff --git a/src/__tests__/matrix/mmulStrassen.test.js b/src/__tests__/matrix/mmulStrassen.test.js
new file mode 100644
index 0000000..98f2dce
--- /dev/null
+++ b/src/__tests__/matrix/mmulStrassen.test.js
@@ -0,0 +1,125 @@
+import { describe, it, expect } from 'vitest';
+
+import { Matrix } from '../..';
+
+// https://github.com/mljs/matrix/issues/114
+describe('mmulStrassen agrees with mmul', () => {
+  const shapes = [
+    [1, 1, 1, 1],
+    [1, 2, 2, 2],
+    [1, 3, 3, 1],
+    [3, 1, 1, 3],
+    [3, 2, 2, 4],
+    [4, 2, 2, 3],
+    [2, 5, 5, 2],
+    [5, 5, 5, 5],
+    [7, 3, 3, 6],
+  ];
+
+  for (const [r1, c1, r2, c2] of shapes) {
+    it(`${r1}x${c1} by ${r2}x${c2}`, () => {
+      const a = Matrix.randInt(r1, c1, { min: 0, max: 9 });
+      const b = Matrix.randInt(r2, c2, { min: 0, max: 9 });
+      const expected = a.mmul(b);
+      const result = a.mmulStrassen(b);
+      expect(result.rows).toBe(r1);
+      expect(result.columns).toBe(c2);
+      expect(result.to2DArray()).toStrictEqual(expected.to2DArray());
+    });
+  }
+
+  it('keeps the operands untouched', () => {
+    const a = new Matrix([
+      [1, 2],
+      [3, 4],
+    ]);
+    const b = new Matrix([
+      [5, 6],
+      [7, 8],
+    ]);
+    a.mmulStrassen(b);
+    expect(a.to2DArray()).toStrictEqual([
+      [1, 2],
+      [3, 4],
+    ]);
+    expect(b.to2DArray()).toStrictEqual([
+      [5, 6],
+      [7, 8],
+    ]);
+  });
+
+  it('accepts a 2D array', () => {
+    const a = new Matrix([[1, 2]]);
+    expect(
+      a
+        .mmulStrassen([
+          [1, 2],
+          [3, 4],
+        ])
+        .to2DArray(),
+    ).toStrictEqual([[7, 10]]);
+  });
+});
+
+describe('mmulStrassen above the recursion threshold', () => {
+  // the recursive path only runs when both dimensions exceed 512, so anything
+  // smaller was silently delegating to mmul and never exercised the block split
+  const sizes = [513, 514, 600];
+
+  for (const n of sizes) {
+    it(`${n}x${n} matches mmul exactly`, () => {
+      // small integers keep every intermediate exact in a float64
+      const a = Matrix.randInt(n, n, { min: 0, max: 3 });
+      const b = Matrix.randInt(n, n, { min: 0, max: 3 });
+      const expected = a.mmul(b);
+      const result = a.mmulStrassen(b);
+      expect(result.rows).toBe(n);
+      expect(result.columns).toBe(n);
+
+      let differing = 0;
+      for (let i = 0; i < n; i++) {
+        for (let j = 0; j < n; j++) {
+          if (result.get(i, j) !== expected.get(i, j)) differing++;
+        }
+      }
+      expect(differing).toBe(0);
+    });
+  }
+
+  it('handles a non square shape above the threshold', () => {
+    const a = Matrix.randInt(520, 600, { min: 0, max: 3 });
+    const b = Matrix.randInt(600, 530, { min: 0, max: 3 });
+    const expected = a.mmul(b);
+    const result = a.mmulStrassen(b);
+    expect(result.rows).toBe(520);
+    expect(result.columns).toBe(530);
+
+    let differing = 0;
+    for (let i = 0; i < 520; i++) {
+      for (let j = 0; j < 530; j++) {
+        if (result.get(i, j) !== expected.get(i, j)) differing++;
+      }
+    }
+    expect(differing).toBe(0);
+  });
+});
+
+describe('mmulStrassen with degenerate matrices', () => {
+  it('a matrix without rows', () => {
+    const result = new Matrix(0, 2).mmulStrassen(new Matrix(2, 3));
+    expect(result.rows).toBe(0);
+    expect(result.columns).toBe(3);
+  });
+
+  it('a matrix without columns', () => {
+    const result = new Matrix(2, 3).mmulStrassen(new Matrix(3, 0));
+    expect(result.rows).toBe(2);
+    expect(result.columns).toBe(0);
+  });
+
+  it('two 0x0 matrices', () => {
+    const result = new Matrix(0, 0).mmulStrassen(new Matrix(0, 0));
+    expect(result.rows).toBe(0);
+    expect(result.columns).toBe(0);
+  });
+});

From 829a842eae0d6c0b778ba18e8889f322cfa152de Mon Sep 17 00:00:00 2001
From: Sarthak Tayal <sarthaktayal2@gmail.com>
Date: Wed, 5 Aug 2026 16:56:23 -0400
Subject: [PATCH 5/6]

---
 matrix.d.ts                               |  6 +++
 src/__tests__/matrix/mmulStrassen.test.js | 47 ++++++++---------------
 2 files changed, 23 insertions(+), 30 deletions(-)

diff --git a/matrix.d.ts b/matrix.d.ts
index c6f6153..6e086dc 100644
--- a/matrix.d.ts
+++ b/matrix.d.ts
@@ -668,6 +668,12 @@ export abstract class AbstractMatrix {
 
   strassen3x3(other: MaybeMatrix): Matrix;
 
+  /**
+   * Returns the same product as {@link AbstractMatrix.mmul}, computed with the Strassen algorithm.
+   * The recursive path only kicks in once both dimensions of the padded operands pass 512,
+   * below that the call is handed to the plain multiplication.
+   * @param y - The right operand.
+   */
   mmulStrassen(y: MaybeMatrix): Matrix;
 
   /**
diff --git a/src/__tests__/matrix/mmulStrassen.test.js b/src/__tests__/matrix/mmulStrassen.test.js
index 98f2dce..ddab8a7 100644
--- a/src/__tests__/matrix/mmulStrassen.test.js
+++ b/src/__tests__/matrix/mmulStrassen.test.js
@@ -62,46 +62,33 @@ describe('mmulStrassen agrees with mmul', () => {
 });
 
 describe('mmulStrassen above the recursion threshold', () => {
-  // the recursive path only runs when both dimensions exceed 512, so anything
-  // smaller was silently delegating to mmul and never exercised the block split
-  const sizes = [513, 514, 600];
+  // the recursive path only runs once both dimensions pass 512, so anything
+  // smaller was delegating to mmul and never exercised the block split.
+  // these two cases cover the odd padded size along with the even one.
+  // small integers keep every intermediate exact in a float64.
+  const cases = [
+    { r1: 513, c1: 513, c2: 513 },
+    { r1: 514, c1: 520, c2: 516 },
+  ];
 
-  for (const n of sizes) {
-    it(`${n}x${n} matches mmul exactly`, () => {
-      // small integers keep every intermediate exact in a float64
-      const a = Matrix.randInt(n, n, { min: 0, max: 3 });
-      const b = Matrix.randInt(n, n, { min: 0, max: 3 });
+  for (const { r1, c1, c2 } of cases) {
+    it(`${r1}x${c1} by ${c1}x${c2} matches mmul exactly`, () => {
+      const a = Matrix.randInt(r1, c1, { min: 0, max: 3 });
+      const b = Matrix.randInt(c1, c2, { min: 0, max: 3 });
       const expected = a.mmul(b);
       const result = a.mmulStrassen(b);
-      expect(result.rows).toBe(n);
-      expect(result.columns).toBe(n);
+      expect(result.rows).toBe(r1);
+      expect(result.columns).toBe(c2);
 
       let differing = 0;
-      for (let i = 0; i < n; i++) {
-        for (let j = 0; j < n; j++) {
+      for (let i = 0; i < r1; i++) {
+        for (let j = 0; j < c2; j++) {
           if (result.get(i, j) !== expected.get(i, j)) differing++;
         }
       }
       expect(differing).toBe(0);
-    });
+    }, 120000);
   }
-
-  it('handles a non square shape above the threshold', () => {
-    const a = Matrix.randInt(520, 600, { min: 0, max: 3 });
-    const b = Matrix.randInt(600, 530, { min: 0, max: 3 });
-    const expected = a.mmul(b);
-    const result = a.mmulStrassen(b);
-    expect(result.rows).toBe(520);
-    expect(result.columns).toBe(530);
-
-    let differing = 0;
-    for (let i = 0; i < 520; i++) {
-      for (let j = 0; j < 530; j++) {
-        if (result.get(i, j) !== expected.get(i, j)) differing++;
-      }
-    }
-    expect(differing).toBe(0);
-  });
 });
 
 describe('mmulStrassen with degenerate matrices', () => {

From 47a2146a06eaa5f96989efce53513ac95297b2f8 Mon Sep 17 00:00:00 2001
From: Sarthak Tayal <sarthaktayal2@gmail.com>
Date: Thu, 6 Aug 2026 16:38:52 -0400
Subject: [PATCH 6/6] fix: avoid square padding rectangular Strassen products

---
 src/__tests__/matrix/mmulStrassen.test.js | 29 ++++++++++++++++++++---
 src/matrix.js                             |  8 +++++++
 2 files changed, 34 insertions(+), 3 deletions(-)

diff --git a/src/__tests__/matrix/mmulStrassen.test.js b/src/__tests__/matrix/mmulStrassen.test.js
index ddab8a7..754d2b2 100644
--- a/src/__tests__/matrix/mmulStrassen.test.js
+++ b/src/__tests__/matrix/mmulStrassen.test.js
@@ -1,6 +1,10 @@
-import { describe, it, expect } from 'vitest';
+import { afterEach, describe, it, expect, vi } from 'vitest';
 
-import { Matrix } from '../..';
+import { AbstractMatrix, Matrix } from '../..';
+
+afterEach(() => {
+  vi.restoreAllMocks();
+});
 
 // https://github.com/mljs/matrix/issues/114
 describe('mmulStrassen agrees with mmul', () => {
@@ -68,7 +72,7 @@ describe('mmulStrassen above the recursion threshold', () => {
   // small integers keep every intermediate exact in a float64.
   const cases = [
     { r1: 513, c1: 513, c2: 513 },
-    { r1: 514, c1: 520, c2: 516 },
+    { r1: 514, c1: 514, c2: 514 },
   ];
 
   for (const { r1, c1, c2 } of cases) {
@@ -89,6 +93,25 @@ describe('mmulStrassen above the recursion threshold', () => {
       expect(differing).toBe(0);
     }, 120000);
   }
+
+  const rectangularCases = [
+    { r1: 1, c1: 513, c2: 1 },
+    { r1: 10, c1: 600, c2: 10 },
+  ];
+
+  for (const { r1, c1, c2 } of rectangularCases) {
+    it(`${r1}x${c1} by ${c1}x${c2} delegates to mmul`, () => {
+      const a = Matrix.randInt(r1, c1, { min: 0, max: 3 });
+      const b = Matrix.randInt(c1, c2, { min: 0, max: 3 });
+      const expected = a.mmul(b);
+      const multiply = vi.spyOn(AbstractMatrix.prototype, 'mmul');
+
+      const result = a.mmulStrassen(b);
+
+      expect(multiply).toHaveBeenCalledTimes(1);
+      expect(result.to2DArray()).toStrictEqual(expected.to2DArray());
+    });
+  }
 });
 
 describe('mmulStrassen with degenerate matrices', () => {
diff --git a/src/matrix.js b/src/matrix.js
index f06e4be..acd90c3 100644
--- a/src/matrix.js
+++ b/src/matrix.js
@@ -1139,6 +1139,14 @@ export class AbstractMatrix {
       return new Matrix(r1, c2);
     }
 
+    // Padding a rectangular product into a square whose side is the largest
+    // dimension can turn a cheap multiplication into an enormous allocation.
+    // The recursive implementation is only useful for equally sized square
+    // operands; mmul already handles every other compatible shape directly.
+    if (c1 === r2 && (r1 !== c1 || r2 !== c2 || r1 !== r2)) {
+      return x.mmul(y);
+    }
+
     // Put a matrix into the top left of a matrix of zeros.
     // `rows` and `cols` are the dimensions of the output matrix.
     function embed(mat, rows, cols) {