Implement matrix inversion

Implement matmul and eye to verify implementation.

src/__tests__/basis.spec.ts

@@ -1,11 +1,18 @@
 import {describe, expect, it} from 'vitest';
 import {
+  eye,
+  fractionalEye,
   fractionalGram,
+  fractionalInv,
   fractionalLenstraLenstraLovasz,
+  fractionalMatmul,
   gram,
+  inv,
   lenstraLenstraLovasz,
+  matmul,
 } from '../basis';
 import {
+  FractionalMonzo,
   applyWeights,
   fractionalDot,
   fractionalMonzosEqual,
@@ -322,3 +329,158 @@ describe('Precise LLL basis reduction', () => {
     }
   });
 });
+
+describe('Matrix multiplication', () => {
+  it('multiplies two matrices', () => {
+    const A = [
+      [1, 0, 1],
+      [2, 1, 1],
+      [0, 1, 1],
+      [1, 1, 2],
+    ];
+    const B = [
+      [1, 2, 1],
+      [2, 3, 1],
+      [4, 2, 2],
+    ];
+    const C = matmul(A, B);
+    expect(C).toEqual([
+      [5, 4, 3],
+      [8, 9, 5],
+      [6, 5, 3],
+      [11, 9, 6],
+    ]);
+  });
+
+  it('multiplies two matrices (fractions)', () => {
+    const A = [
+      [1, 0, 1],
+      [2, 1, 1],
+      [0, 0.5, 1],
+      [1, 1, 2],
+    ];
+    const B = [
+      [1, 2, 1],
+      [2, 3, 1],
+      ['4/3', 2, 2],
+    ];
+    const C = fractionalMatmul(A, B);
+    expect(C.map(row => row.map(f => f.toFraction()))).toEqual([
+      ['7/3', '4', '3'],
+      ['16/3', '9', '5'],
+      ['7/3', '7/2', '5/2'],
+      ['17/3', '9', '6'],
+    ]);
+  });
+});
+
+describe('Matrix inverse', () => {
+  it('computes a 3x3 inverse', () => {
+    const mat = [
+      [2, -1], // Missing entry interpreted as 0
+      [-1, 2, -1],
+      [0, -1, 2],
+    ];
+    const inverse = inv(mat).map(row => row.map(x => Math.round(4 * x) / 4));
+    expect(inverse).toEqual([
+      [0.75, 0.5, 0.25],
+      [0.5, 1, 0.5],
+      [0.25, 0.5, 0.75],
+    ]);
+  });
+
+  it('throws for non-square matrix', () => {
+    expect(() =>
+      inv([
+        [1, 2],
+        [3, 4],
+        [5, 6],
+      ])
+    ).toThrow('Non-square matrix');
+  });
+
+  it('throws for singular matrix', () => {
+    expect(() =>
+      inv([
+        [1, 0],
+        [0, 0],
+      ])
+    ).toThrow('Matrix is singular');
+  });
+
+  it('computes a 3x3 with fractional entries', () => {
+    const mat = [
+      [2, -1], // Missing entry interpreted as 0
+      [-1, '2/1', -1],
+      [0, -1, '2'],
+    ];
+    const inverse = fractionalInv(mat).map(row => row.map(x => x.toFraction()));
+    expect(inverse).toEqual([
+      ['3/4', '1/2', '1/4'],
+      ['1/2', '1', '1/2'],
+      ['1/4', '1/2', '3/4'],
+    ]);
+  });
+
+  it('throws for non-square matrix with fractional entries', () => {
+    expect(() =>
+      fractionalInv([
+        [1, 2],
+        [3, 4],
+        [5, 6],
+      ])
+    ).toThrow('Non-square matrix');
+  });
+
+  it('throws for singular matrix with fractional entries', () => {
+    expect(() =>
+      fractionalInv([
+        [1, 0],
+        [0, 0],
+      ])
+    ).toThrow('Matrix is singular');
+  });
+
+  it.runIf(FUZZ)('fuzzes for random inverses', () => {
+    for (let k = 0; k < 10000; ++k) {
+      const mat: number[][] = [];
+      const N = Math.ceil(1 + Math.random() * 10);
+      for (let i = 0; i < N; ++i) {
+        const row: number[] = [];
+        for (let j = 0; j < N; ++j) {
+          row.push(Math.random() * 10 - 5);
+        }
+        mat.push(row);
+      }
+      const inverse = inv(mat);
+      const I = matmul(mat, inverse).map(row =>
+        row.map(x => Math.round(x * 1024) / 1024 || 0)
+      );
+      expect(I).toEqual(eye(N));
+    }
+  });
+
+  it.runIf(FUZZ)('fuzzes for random inverses (fractional)', () => {
+    for (let k = 0; k < 1000; ++k) {
+      const mat: number[][] = [];
+      const N = Math.ceil(1 + Math.random() * 7);
+      for (let i = 0; i < N; ++i) {
+        const row: number[] = [];
+        for (let j = 0; j < N; ++j) {
+          row.push(Math.round(Math.random() * 10 - 5));
+        }
+        mat.push(row);
+      }
+      let inverse: FractionalMonzo[] | undefined;
+      try {
+        inverse = fractionalInv(mat);
+      } catch {
+        /** empty */
+      }
+      if (inverse) {
+        const I = fractionalMatmul(mat, inverse);
+        expect(I).toEqual(fractionalEye(N));
+      }
+    }
+  });
+});