diff --git a/src/__tests__/basis.spec.ts b/src/__tests__/basis.spec.ts
index f99892f..1541076 100644
--- a/src/__tests__/basis.spec.ts
+++ b/src/__tests__/basis.spec.ts
@@ -1,11 +1,22 @@
 import {describe, expect, it} from 'vitest';
 import {
+  det,
+  eye,
+  fractionalDet,
+  fractionalEye,
   fractionalGram,
+  fractionalInv,
   fractionalLenstraLenstraLovasz,
+  fractionalMatmul,
+  fractionalTranspose,
   gram,
+  inv,
   lenstraLenstraLovasz,
+  matmul,
+  transpose,
 } from '../basis';
 import {
+  FractionalMonzo,
   applyWeights,
   fractionalDot,
   fractionalMonzosEqual,
@@ -322,3 +333,255 @@ 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('computes another 3x3 inverse', () => {
+    const mat = [
+      [-2, -1, 2],
+      [2, 1, 4],
+      [-3, 3, -1],
+    ];
+    const inverse = inv(mat).map(row => row.map(x => Math.round(100000 * x)));
+    expect(inverse).toEqual([
+      [-24074, 9259, -11111],
+      [-18519, 14815, 22222],
+      [16667, 16667, 0],
+    ]);
+  });
+
+  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('computes another 3x3 inverse with fractional result', () => {
+    const mat = [
+      [-2, -1, 2],
+      [2, 1, 4],
+      [-3, 3, -1],
+    ];
+    const inverse = fractionalInv(mat).map(row => row.map(x => x.toFraction()));
+    expect(inverse).toEqual([
+      ['-13/54', '5/54', '-1/9'],
+      ['-5/27', '4/27', '2/9'],
+      ['1/6', '1/6', '0'],
+    ]);
+  });
+
+  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));
+      }
+    }
+  });
+});
+
+describe('Determinant', () => {
+  it('computes the determinant of a 3x3 matrix', () => {
+    const mat = [
+      [-2, -1, 2],
+      [2, 1, 4],
+      [-3, 3, -1],
+    ];
+    const determinant = det(mat);
+    expect(determinant).toBe(54);
+  });
+
+  it('computes 0 for the origin', () => {
+    const mat = [[0, 0], []];
+    const determinant = det(mat);
+    expect(determinant).toBe(0);
+  });
+
+  it('computes the area of a square', () => {
+    const sides = [
+      [1, 1],
+      [-1, 1],
+    ];
+    const determinant = det(sides);
+    expect(determinant).toBe(2);
+  });
+
+  it('computes the volume of a skew-aligned box', () => {
+    const sides = [
+      [1, 1, 1],
+      [-Math.SQRT2, Math.SQRT1_2, Math.SQRT1_2],
+      [0, 1, -1],
+    ];
+    const determinant = det(sides);
+    expect(determinant).toBe(-3 * Math.SQRT2);
+  });
+
+  it('computes the determinant of a 3x3 matrix with fractional entries', () => {
+    const mat = [
+      [-2, -1, 2],
+      [2, 0.5, 4],
+      [-3, 3, '-1/3'],
+    ];
+    const d = fractionalDet(mat);
+    expect(d.toFraction()).toBe('152/3');
+  });
+
+  // TODO: Fuzz against a naïve implementation.
+});
+
+describe('Transpose', () => {
+  it('transposes a 3x2 matrix', () => {
+    const mat = [[1, 2], [3], [4, 5]];
+    expect(transpose(mat)).toEqual([
+      [1, 3, 4],
+      [2, 0, 5],
+    ]);
+  });
+
+  it('transposes a 3x2 matrix with rational entries', () => {
+    const mat = [[1, 0.5], [3], ['2/7', 5]];
+    expect(
+      fractionalTranspose(mat).map(row => row.map(f => f.toFraction()))
+    ).toEqual([
+      ['1', '3', '2/7'],
+      ['1/2', '0', '5'],
+    ]);
+  });
+});
diff --git a/src/basis.ts b/src/basis.ts
index a8ed0f4..a8d30cb 100644
--- a/src/basis.ts
+++ b/src/basis.ts
@@ -218,3 +218,436 @@ export function fractionalLenstraLenstraLovasz(
     },
   };
 }
+
+/**
+ * Return a 2-D array with ones on the diagonal and zeros elsewhere.
+ * @param N Number of rows in the output.
+ * @param M Number of columns in the output.
+ * @param k Index of the diagonal.
+ * @returns An array where all elements are equal to zero, except for the `k`-th diagonal, whose values are equal to one.
+ */
+export function eye(N: number, M?: number, k = 0) {
+  M ??= N;
+  const result: number[][] = [];
+  for (let i = 0; i < N; ++i) {
+    result.push(Array(M).fill(0));
+    if (i + k < M) {
+      result[i][i + k] = 1;
+    }
+  }
+  return result;
+}
+
+/**
+ * Return a 2-D array with ones on the diagonal and zeros elsewhere.
+ * @param N Number of rows in the output.
+ * @param M Number of columns in the output.
+ * @param k Index of the diagonal.
+ * @returns An array where all elements are equal to zero, except for the `k`-th diagonal, whose values are equal to one.
+ */
+export function fractionalEye(N: number, M?: number, k = 0): FractionalMonzo[] {
+  M ??= N;
+  const result: FractionalMonzo[] = [];
+  for (let i = 0; i < N; ++i) {
+    const row: FractionalMonzo = [];
+    for (let j = 0; j < M; ++j) {
+      if (j === i + k) {
+        row.push(new Fraction(1));
+      } else {
+        row.push(new Fraction(0));
+      }
+    }
+    result.push(row);
+  }
+  return result;
+}
+
+// XXX: I'm sorry. This matrix inversion algorithm is not particularly sophisticated. Existing solutions just come with too much bloat.
+
+/**
+ * Compute the (multiplicative) inverse of a matrix.
+ * @param matrix Matrix to be inverted.
+ * @returns The multiplicative inverse.
+ * @throws An error if the matrix is not square or not invertible.
+ */
+export function inv(matrix: number[][]) {
+  let width = 0;
+  const height = matrix.length;
+  for (const row of matrix) {
+    width = Math.max(width, row.length);
+  }
+  if (width !== height) {
+    throw new Error('Non-square matrix');
+  }
+  const result: number[][] = [];
+  for (let i = 0; i < height; ++i) {
+    result.push(Array(width).fill(0));
+    result[i][i] = 1;
+  }
+  // Don't modify input
+  matrix = matrix.map(row => [...row]);
+  // Coerce missing entries to zeros
+  for (let y = 0; y < height; ++y) {
+    for (let x = matrix[y].length; x < width; ++x) {
+      matrix[y][x] = 0;
+    }
+  }
+  // Put ones along the diagonal, zeros in the lower triangle
+  for (let x = 0; x < width; ++x) {
+    let s = matrix[x][x];
+    if (!s) {
+      // Row echelon form
+      for (let y = x + 1; y < height; ++y) {
+        if (matrix[y][x]) {
+          let temp = matrix[y];
+          matrix[y] = matrix[x];
+          matrix[x] = temp;
+
+          temp = result[y];
+          result[y] = result[x];
+          result[x] = temp;
+          break;
+        }
+      }
+      s = matrix[x][x];
+      if (!s) {
+        throw new Error('Matrix is singular');
+      }
+    }
+    if (s !== 1) {
+      s = 1 / s;
+      matrix[x] = matrix[x].map(a => a * s);
+      result[x] = result[x].map(a => a * s);
+    }
+    for (let y = x + 1; y < height; ++y) {
+      s = matrix[y][x];
+      if (s) {
+        matrix[y] = matrix[y].map((a, i) => a - s * matrix[x][i]);
+        result[y] = result[y].map((a, i) => a - s * result[x][i]);
+      }
+    }
+  }
+  // Eliminate remaining entries in the upper triangle
+  for (let x = width - 1; x > 0; --x) {
+    for (let y = x - 1; y >= 0; --y) {
+      const s = matrix[y][x];
+      if (s) {
+        // No need to keep track of these entries anymore.
+        // matrix[y] = matrix[y].map((a, i) => a - s * matrix[x][i]);
+        result[y] = result[y].map((a, i) => a - s * result[x][i]);
+      }
+    }
+  }
+  return result;
+}
+
+/**
+ * Compute the (multiplicative) inverse of a matrix.
+ * @param matrix Matrix to be inverted.
+ * @returns The multiplicative inverse.
+ * @throws An error if the matrix is not square or not invertible.
+ */
+export function fractionalInv(matrix: ProtoFractionalMonzo[]) {
+  let width = 0;
+  const height = matrix.length;
+  for (const row of matrix) {
+    width = Math.max(width, row.length);
+  }
+  if (width !== height) {
+    throw new Error('Non-square matrix');
+  }
+  const result: FractionalMonzo[] = [];
+  for (let i = 0; i < height; ++i) {
+    const row: FractionalMonzo = [];
+    for (let j = 0; j < width; ++j) {
+      if (i === j) {
+        row.push(new Fraction(1));
+      } else {
+        row.push(new Fraction(0));
+      }
+    }
+    result.push(row);
+  }
+  // Don't modify input
+  const matrix_: FractionalMonzo[] = matrix.map(row =>
+    row.map(f => new Fraction(f))
+  );
+  // Coerce missing entries to zeros
+  for (let y = 0; y < height; ++y) {
+    for (let x = matrix_[y].length; x < width; ++x) {
+      matrix_[y][x] = new Fraction(0);
+    }
+  }
+  // Put ones along the diagonal, zeros in the lower triangle
+  for (let x = 0; x < width; ++x) {
+    let s = matrix_[x][x];
+    if (!s.n) {
+      // Row echelon form
+      for (let y = x + 1; y < height; ++y) {
+        if (matrix_[y][x].n) {
+          let temp = matrix_[y];
+          matrix_[y] = matrix_[x];
+          matrix_[x] = temp;
+
+          temp = result[y];
+          result[y] = result[x];
+          result[x] = temp;
+          break;
+        }
+      }
+      s = matrix_[x][x];
+      if (!s.n) {
+        throw new Error('Matrix is singular');
+      }
+    }
+    if (!s.isUnity()) {
+      matrix_[x] = matrix_[x].map(a => a.div(s));
+      result[x] = result[x].map(a => a.div(s));
+    }
+    for (let y = x + 1; y < height; ++y) {
+      s = matrix_[y][x];
+      if (s.n) {
+        matrix_[y] = matrix_[y].map((a, i) => a.sub(s.mul(matrix_[x][i])));
+        result[y] = result[y].map((a, i) => a.sub(s.mul(result[x][i])));
+      }
+    }
+  }
+  // Eliminate remaining entries in the upper triangle
+  for (let x = width - 1; x > 0; --x) {
+    for (let y = x - 1; y >= 0; --y) {
+      const s = matrix_[y][x];
+      if (s.n) {
+        // No need to keep track of these entries anymore.
+        // matrix_[y] = matrix_[y].map(...);
+        result[y] = result[y].map((a, i) => a.sub(s.mul(result[x][i])));
+      }
+    }
+  }
+  return result;
+}
+
+/**
+ * Compute the matrix product of two arrays of arrays of numbers.
+ * @param A The left operand.
+ * @param B The right operand.
+ * @returns The matrix product of the operands.
+ */
+export function matmul(A: number[][], B: number[][]) {
+  const height = A.length;
+  let width = 0;
+  for (const row of B) {
+    width = Math.max(width, row.length);
+  }
+  let n = 0;
+  for (const row of A) {
+    n = Math.max(n, row.length);
+  }
+  B = [...B];
+  while (B.length < n) {
+    B.push([]);
+  }
+  const result: number[][] = [];
+  for (let i = 0; i < height; ++i) {
+    const row = Array(width).fill(0);
+    const rowA = A[i];
+    for (let j = 0; j < width; ++j) {
+      for (let k = 0; k < rowA.length; ++k) {
+        row[j] += rowA[k] * (B[k][j] ?? 0);
+      }
+    }
+    result.push(row);
+  }
+  return result;
+}
+
+/**
+ * Compute the matrix product of two arrays of arrays of fractions.
+ * @param A The left operand.
+ * @param B The right operand.
+ * @returns The matrix product of the operands.
+ */
+export function fractionalMatmul(
+  A: ProtoFractionalMonzo[],
+  B: ProtoFractionalMonzo[]
+): FractionalMonzo[] {
+  const height = A.length;
+  let width = 0;
+  for (const row of B) {
+    width = Math.max(width, row.length);
+  }
+  let n = 0;
+  for (const row of A) {
+    n = Math.max(n, row.length);
+  }
+  const B_ = B.map(row => row.map(f => new Fraction(f)));
+  while (B_.length < n) {
+    B_.push([]);
+  }
+  for (const row of B_) {
+    while (row.length < width) {
+      row.push(new Fraction(0));
+    }
+  }
+  const result: FractionalMonzo[] = [];
+  for (let i = 0; i < height; ++i) {
+    const row: FractionalMonzo = [];
+    while (row.length < width) {
+      row.push(new Fraction(0));
+    }
+    const rowA = A[i].map(f => new Fraction(f));
+    for (let j = 0; j < width; ++j) {
+      for (let k = 0; k < rowA.length; ++k) {
+        row[j] = row[j].add(rowA[k].mul(B_[k][j]));
+      }
+    }
+    result.push(row);
+  }
+  return result;
+}
+
+/**
+ * Compute the determinant of a matrix.
+ * @param matrix Array of arrays of numbers to calculate determinant for.
+ * @returns The determinant.
+ */
+export function det(matrix: number[][]) {
+  let width = 0;
+  const height = matrix.length;
+  for (const row of matrix) {
+    width = Math.max(width, row.length);
+  }
+  if (width !== height) {
+    throw new Error('Non-square matrix');
+  }
+  matrix = matrix.map(row => [...row]);
+  let result = 1;
+  for (let x = 0; x < width; ++x) {
+    let d = matrix[x][x];
+    if (!d) {
+      // Row echelon form
+      for (let y = x + 1; y < height; ++y) {
+        if (matrix[y][x]) {
+          const temp = matrix[y];
+          matrix[y] = matrix[x];
+          matrix[x] = temp;
+
+          result = -result;
+          break;
+        }
+      }
+      d = matrix[x][x];
+      if (!d) {
+        return 0;
+      }
+    }
+    result *= d;
+    d = 1 / d;
+    for (let y = x + 1; y < height; ++y) {
+      const row = matrix[y];
+      const s = row[x] * d;
+      if (s) {
+        // Skip over entries that are not used later.
+        const upperRow = matrix[x];
+        for (let i = x + 1; i < upperRow.length; ++i) {
+          row[i] -= s * upperRow[i];
+        }
+        // Full row operation for reference:
+        // matrix[y] = matrix[y].map((a, i) => a - s * (matrix[x][i] ?? 0));
+      }
+    }
+  }
+  return result;
+}
+
+/**
+ * Compute the determinant of a matrix with rational entries.
+ * @param matrix Array of arrays of fractions to calculate determinant for.
+ * @returns The determinant.
+ */
+export function fractionalDet(matrix: ProtoFractionalMonzo[]) {
+  let width = 0;
+  const height = matrix.length;
+  for (const row of matrix) {
+    width = Math.max(width, row.length);
+  }
+  if (width !== height) {
+    throw new Error('Non-square matrix');
+  }
+  const matrix_ = matrix.map(row => row.map(f => new Fraction(f)));
+  let result = new Fraction(1);
+  for (let x = 0; x < width; ++x) {
+    let d = matrix_[x][x];
+    if (d === undefined || !d.n) {
+      // Row echelon form
+      for (let y = x + 1; y < height; ++y) {
+        if (matrix_[y][x]) {
+          const temp = matrix_[y];
+          matrix_[y] = matrix_[x];
+          matrix_[x] = temp;
+
+          result = result.neg();
+          break;
+        }
+      }
+      d = matrix_[x][x];
+      if (d === undefined || !d.n) {
+        return new Fraction(0);
+      }
+    }
+    result = result.mul(d);
+    for (let y = x + 1; y < height; ++y) {
+      const row = matrix_[y];
+      const s = row[x].div(d);
+      if (s.n) {
+        const upperRow = matrix_[x];
+        for (let i = x + 1; i < upperRow.length; ++i) {
+          row[i] = row[i].sub(s.mul(upperRow[i]));
+        }
+        // matrix_[y] = matrix_[y].map((a, i) => a.sub(s.mul(matrix_[x][i] ?? 0)));
+      }
+    }
+  }
+  return result;
+}
+
+/**
+ * Transpose a 2-D matrix (swap rows and columns).
+ * @param matrix Matrix to transpose.
+ * @returns The transpose.
+ */
+export function transpose(matrix: number[][]) {
+  let width = 0;
+  for (const row of matrix) {
+    width = Math.max(row.length, width);
+  }
+  const result: number[][] = [];
+  for (let i = 0; i < width; ++i) {
+    const row: number[] = [];
+    for (let j = 0; j < matrix.length; ++j) {
+      row.push(matrix[j][i] ?? 0);
+    }
+    result.push(row);
+  }
+  return result;
+}
+
+/**
+ * Transpose a 2-D matrix with rational entries (swap rows and columns).
+ * @param matrix Matrix to transpose.
+ * @returns The transpose.
+ */
+export function fractionalTranspose(matrix: ProtoFractionalMonzo[]) {
+  let width = 0;
+  for (const row of matrix) {
+    width = Math.max(row.length, width);
+  }
+  const result: FractionalMonzo[] = [];
+  for (let i = 0; i < width; ++i) {
+    const row: FractionalMonzo = [];
+    for (let j = 0; j < matrix.length; ++j) {
+      row.push(new Fraction(matrix[j][i] ?? 0));
+    }
+    result.push(row);
+  }
+  return result;
+}