Find inverse of matrix using integers and fractions

inverse_int_frac_matrix.py

from typing import Tuple
import numpy as np
from fractions import Fraction
import math


def invert_integer_matrix(
    matrix: Tuple[Tuple[int]],
) -> Tuple[Tuple[int], Tuple[Tuple[int]]]:
    """
    Given a matrix A with integer entries, this function computes the
    inverse of this matrix by gaussian elimination.
    
    # Input: 
        matrix: Tuple[Tuple[int]]
            Square matrix A of size (m, m) of integer values
            
    # Output:
        diagonal: Tuple[int]
            The final diagonal D after gaussian elimination, with values d_i
        inverse: Tuple[Tuple[int]]
            The final inversed matrix M = diag(D) * A^{-1}
    """
    side = len(matrix)
    inverse = np.eye(side, dtype="object")
    matrix = np.column_stack((matrix, inverse))
    for k in range(side):
        for i in range(k + 1, side):
            matrix[i] = matrix[i] * matrix[k, k] - matrix[k] * matrix[i, k]
            gdcline = math.gcd(*matrix[i])
            if gdcline != 1:
                matrix[i] = matrix[i] // gdcline
    for k in range(side - 1, 0, -1):
        for i in range(k - 1, -1, -1):
            matrix[i] = matrix[i] * matrix[k, k] - matrix[k] * matrix[i, k]
            gdcline = math.gcd(*matrix[i])
            if gdcline != 1:
                matrix[i] = matrix[i] // gdcline
    diagonal = np.diag(matrix[:, :side])
    inverse = matrix[:, side:]
    return diagonal, inverse


def invert_fraction_matrix(matrix: Tuple[Tuple[Fraction]]) -> Tuple[Tuple[Fraction]]:
    """
    Given a matrix A with fraction entries, this function computes the
    inverse of this matrix by gaussian elimination.
    
    # Input: 
        matrix: Tuple[Tuple[Fraction]]
            Square matrix A of size (m, m) of Fraction entities
            
    # Output:
        inverse: Tuple[Tuple[Fration]]
            The final inversed matrix M = A^{-1}
    """
    matrix = np.array(matrix, dtype="object")
    side = len(matrix)
    intdiag = [1] * side
    for i, line in enumerate(matrix):
        denoms = [frac.denominator for frac in line]
        intdiag[i] = math.lcm(*denoms)
        matrix[i] *= intdiag[i]
    intdiag = np.array(intdiag, dtype="int64")
    integermatrix = np.array(matrix, dtype="int64")
    diag, inverseintegermat = invert_integer_matrix(integermatrix)
    inverse = np.array(inverseintegermat, dtype="object")
    for i, line in enumerate(inverseintegermat):
        for j, elem in enumerate(line):
            inverse[i, j] = Fraction(elem * intdiag[j], diag[i])
    return inverse


if __name__ == "__main__":
    matrix = [[1, 0], [0, 1]]
    diagonal, inverse = invert_integer_matrix(matrix)
    np.testing.assert_allclose(diagonal, np.ones(len(matrix)))
    np.testing.assert_allclose(inverse, matrix)

    matrix = [[3, 4], [1, 2]]
    diagonal, inverse = invert_integer_matrix(matrix)
    testident = inverse @ matrix
    testident = testident / diagonal
    np.testing.assert_allclose(testident, np.eye(len(matrix)))

    matrix = [[1, 2], [2, 3]]
    diagonal, inverse = invert_integer_matrix(matrix)
    testident = inverse @ matrix
    testident = testident / diagonal
    np.testing.assert_allclose(testident, np.eye(len(matrix)))

    frac = Fraction
    matrix = [[frac(1), frac(-2)], [frac(-1, 2), frac(3, 2)]]
    matrix = np.array(matrix, dtype="object")
    inverse = invert_fraction_matrix(matrix)
    
    size = 3
    floatmatrix = np.random.uniform(-1, 1, (size, size))
    fracmatrix = np.empty((size, size), dtype="object")
    for i, line in enumerate(floatmatrix):
        for j, elem in enumerate(line):
            fracmatrix[i, j] = Fraction(elem).limit_denominator(10)
    fracmatrix += np.transpose(fracmatrix)
    invfracmatrix = invert_fraction_matrix(fracmatrix)
    product = fracmatrix @ invfracmatrix
    product = np.array(product, dtype="float64")
    np.testing.assert_allclose(product, np.eye(size))