Pure python implementation of the Bareiss algorithm

bareiss.py

def adjugate(a, m=None):
    # Initialize
    sign = 1
    previous = pivot = 1

    # Bareiss formula
    def do_pivot(a, b, c, d, e):
        x = a * d - b * c
        if m is None:
            q, r = divmod(x, e)
            assert r == 0
            return q
        q = x * pow(e, -1, m)
        return q % m

    # Assert square matrix
    n = len(a)
    assert all(len(row) == n for row in a)

    # Build augmented matrix
    am = [
        list(map(int, row)) + [0] * i + [1] + [0] * (n - i - 1)
        for i, row in enumerate(a)
    ]

    # Iterate diagonally
    for c in range(n):

        # Test pivot
        if am[c][c] == 0:

            # Find pivot
            for i in range(c, n):
                if am[i][c] != 0:
                    break

            # Non-invertible detection
            else:
                return 0, [[0] * n for _ in range(n)]

            # Switch rows
            am[c], am[i] = am[i], am[c]
            sign *= -1

        # Get pivot
        previous, pivot = pivot, am[c][c]

        # Set zeros
        for i, row in enumerate(am):
            if i == c:
                continue
            am[i] = [
                do_pivot(pivot, y, row[c], x, previous) for x, y in zip(row, am[c])
            ]

    # Determinant
    determinant = sign * pivot
    if m is not None:
        determinant %= m

    # Return
    return determinant, [row[n:] for row in am]

Hi Michel,
Could you confirm that 'm' is the size of a finite field?
Thank you

@orah29 Yes, every time the matrix is updated the result is taken modulo m. I used it mostly with m=2 since the addition in Z/2 corresponds to the XOR operation over bits.

With A = np.eye(9) + np.arange(9) there is no pb when I test adjugate(A) (=equal to np.linalg.det(A) * np.linalg.inv(A)) but adjugate(A, 256) stop on an AssertionError (r != 0). I have to deal with the ff Z/2^8. Obviously, detA % 256 != 0. I translated your code in numpy, but the pb is the same. Could you help me? Thank you per advance.

Sorry I confused Z/p with Z/p^n ... I going to use an irreducible polynomial to calculate on Z/256

@orah29 Oh right I made a mistake.

When m is provided, the new values should be computed not as (a * d - b * c) // e but as (a * d - b * c) * modinv(e, m) where modinv is the modular multiplicative inverse. Still, it requires that e and m are coprime, meaning m should be prime.