feat: add determinant algorithm

Author: cosmicproc

maths/determinant.ts

@@ -0,0 +1,73 @@
+/**
+ * @description
+ * Computes the determinant of the given matrix using elimination.
+ * - Rounding errors may occur for some matrices.
+ * - Only handles 6 decimal places. Rounds thereafter.
+ * @Complexity_Analysis
+ * Time complexity: O(n^3)
+ * Space Complexity: O(n^2)
+ * @param {number[][]} m - A square matrix (2D array)
+ * @return {number} - The determinant
+ * @example det([[1,1],[1,1]]) = 0
+ */
+
+function interchange(m: number[][], from: number, to: number): number[][] {
+  ;[m[to], m[from]] = [m[from], m[to]]
+  return m
+}
+
+function addition(
+  m: number[][],
+  from: number,
+  to: number,
+  c: number
+): number[][] {
+  m[to] = m[to].map((e, i) => e + c * m[from][i])
+  return m
+}
+
+function diagProduct(m: number[][]): number {
+  let product = 1
+  for (let i = 0; i < m.length; i++) {
+    product *= m[i][i]
+  }
+  return product
+}
+
+export function det(m: number[][]): number {
+  if (m.some((r) => r.length != m.length)) {
+    throw new Error('only square matrices can have determinants')
+  }
+
+  const decPlaces = 6
+  const epsilon = 1e-6
+
+  // track the number of applied interchange operations
+  let appliedICs = 0
+  for (let i = 0; i < m[0].length; i++) {
+    // partial pivotting
+    let idealPivot = null
+    let maxValue = 0
+    for (let j = i; j < m.length; j++) {
+      if (Math.abs(m[j][i]) > maxValue) {
+        maxValue = Math.abs(m[j][i])
+        idealPivot = j
+      }
+    }
+    if (idealPivot === null) {
+      return 0
+    }
+    if (idealPivot != i) {
+      m = interchange(m, i, idealPivot)
+      appliedICs++
+    }
+    // eliminate entries under the pivot
+    for (let j = i + 1; j < m.length; j++) {
+      if (Math.abs(m[j][i]) > epsilon) {
+        m = addition(m, i, j, -m[j][i] / m[i][i])
+      }
+    }
+  }
+  const result = diagProduct(m) * (-1) ** appliedICs
+  return parseFloat(result.toFixed(decPlaces))
+}

maths/test/determinant.test.ts

@@ -0,0 +1,89 @@
+import { det } from '../determinant'
+
+describe('determinant', () => {
+  test.each([
+    [
+      [
+        [1, 2],
+        [3, 4, 5]
+      ]
+    ],
+    [
+      [
+        [1, 2, 3],
+        [4, 5, 6]
+      ]
+    ],
+    [
+      [
+        [1, 2],
+        [3, 4],
+        [5, 6]
+      ]
+    ]
+  ])('should throw an error for non square matrix %p', (matrix) => {
+    expect(() => det(matrix)).toThrow(
+      'only square matrices can have determinants'
+    )
+  })
+
+  test.each([
+    [
+      [
+        [1, 2],
+        [3, 4]
+      ],
+      -2
+    ],
+    [
+      [
+        [1, 1],
+        [1, 1]
+      ],
+      0
+    ],
+    [
+      [
+        [1, 2],
+        [0, 0]
+      ],
+      0
+    ],
+    [
+      [
+        [8, 1, 5],
+        [9, 3, 7],
+        [1, 4, 4]
+      ],
+      8
+    ],
+    [
+      [
+        [15, 85, 32],
+        [76, 83, 23],
+        [28, 56, 92]
+      ],
+      -382536
+    ],
+    [
+      [
+        [2, -1, 0, 3],
+        [4, 0, 1, 2],
+        [3, 2, -1, 1],
+        [1, 3, 2, -2]
+      ],
+      -42
+    ],
+    [
+      [
+        [0.75483643, 0.68517541, 0.53548329, 0.5931435],
+        [0.37031247, 0.80103707, 0.82563949, 0.91266224],
+        [0.39293451, 0.27228353, 0.54093836, 0.51963319],
+        [0.60997323, 0.40161682, 0.58330774, 0.17392144]
+      ],
+      -0.051073
+    ]
+  ])('determinant of %p should be %d', (matrix, expected) => {
+    expect(det(matrix)).toBe(expected)
+  })
+})