מטריצות הפוכות

function inverse(matrix, relativeTolerance=1e-12) {
  const n = matrix.length;
  if (!n || matrix.some((row) => row.length !== n)) {
    throw new Error("A must be a nonempty square matrix");
  }

  const numericMatrix = matrix.map((row) => row.map(Number));
  const rowScales = numericMatrix.map((row) => Math.max(...row.map(Math.abs)));
  if (rowScales.some((scale) => scale === 0)) return null;

  const augmented = numericMatrix.map((row, i) => [
    ...row,
    ...Array.from({ length: n }, (_, j) => Number(i === j)),
  ]);

  for (let col = 0; col < n; col += 1) {
    let pivot = col;
    for (let row = col + 1; row < n; row += 1) {
      const candidateRatio = Math.abs(augmented[row][col]) / rowScales[row];
      const pivotRatio = Math.abs(augmented[pivot][col]) / rowScales[pivot];
      if (candidateRatio > pivotRatio) {
        pivot = row;
      }
    }
    if (augmented[pivot][col] === 0) return null;
    if (Math.abs(augmented[pivot][col]) <= relativeTolerance*rowScales[pivot]) {
      throw new Error("nonzero pivot is numerically unsafe at this tolerance");
    }
    [augmented[col], augmented[pivot]] = [augmented[pivot], augmented[col]];
    [rowScales[col], rowScales[pivot]] = [rowScales[pivot], rowScales[col]];

    const scale = augmented[col][col];
    augmented[col] = augmented[col].map((value) => value/scale);
    for (let row = 0; row < n; row += 1) {
      if (row === col) continue;
      const factor = augmented[row][col];
      augmented[row] = augmented[row].map(
        (value, j) => value - factor*augmented[col][j],
      );
    }
  }
  return augmented.map((row) => row.slice(n));
}

function multiply(A, B) {
  return A.map((row) =>
    B[0].map((_, j) =>
      row.reduce((sum, value, k) => sum + value*B[k][j], 0),
    ),
  );
}

const A = [[2, 1, 0], [1, 2, 1], [0, 1, 2]];
const Ainv = inverse(A);
console.table(Ainv);
console.log({
  A_times_Ainv: multiply(A, Ainv),
  Ainv_times_A: multiply(Ainv, A),
});
console.log(inverse([[1, 2], [2, 4]]));
console.log(inverse([[1, 0], [0, 1e-13]]));

def inverse(matrix, relative_tolerance=1e-12):
    n = len(matrix)
    if not n or any(len(row) != n for row in matrix):
        raise ValueError("A must be a nonempty square matrix")

    numeric_matrix = [[float(value) for value in row] for row in matrix]
    row_scales = [max(abs(value) for value in row) for row in numeric_matrix]
    if any(scale == 0 for scale in row_scales):
        return None

    augmented = [
        row[:]
        + [float(i == j) for j in range(n)]
        for i, row in enumerate(numeric_matrix)
    ]

    for col in range(n):
        pivot = max(
            range(col, n),
            key=lambda row: abs(augmented[row][col])/row_scales[row],
        )
        if augmented[pivot][col] == 0:
            return None
        if abs(augmented[pivot][col]) <= relative_tolerance*row_scales[pivot]:
            raise ArithmeticError(
                "nonzero pivot is numerically unsafe at this tolerance"
            )
        augmented[col], augmented[pivot] = augmented[pivot], augmented[col]
        row_scales[col], row_scales[pivot] = row_scales[pivot], row_scales[col]

        scale = augmented[col][col]
        augmented[col] = [value/scale for value in augmented[col]]
        for row in range(n):
            if row == col:
                continue
            factor = augmented[row][col]
            augmented[row] = [
                value - factor*pivot_value
                for value, pivot_value
                in zip(augmented[row], augmented[col])
            ]
    return [row[n:] for row in augmented]

def multiply(A, B):
    return [
        [
            sum(a*b for a, b in zip(row, column))
            for column in zip(*B)
        ]
        for row in A
    ]

A = [[2, 1, 0], [1, 2, 1], [0, 1, 2]]
A_inv = inverse(A)
print(A_inv)
print({
    "A_times_A_inv": multiply(A, A_inv),
    "A_inv_times_A": multiply(A_inv, A),
})
print(inverse([[1, 2], [2, 4]]))
print(inverse([[1, 0], [0, 1e-13]]))