A benchmark comparing sympy and a pure python implementation of the Bareiss algorithm

benchmark.py

import sys
import time
import hashlib
from fractions import Fraction
from contextlib import contextmanager

from sympy import Matrix

# https://gist.github.com/vxgmichel/080e9999a1020711f27cd60b5c2d14de
from bareiss import adjugate


def bitstring(i, n=256):
    name = '{}\0'.format(i).encode()
    h = hashlib.sha256(name).digest()
    i = int.from_bytes(h, 'big')
    g = map(int, '{:0256b}'.format(i))
    return list(g)[:n]


@contextmanager
def timeit(*labels):
    start = time.time()
    yield
    delta = time.time() - start
    label = '[' + ', '.join(labels) + ']'
    print('{:<25} >>> {:.3f} seconds'.format(label, delta))


def benchmark(n=256):
    a = [bitstring(i, n) for i in range(n)]
    m1 = list(map(list, zip(*a)))
    m2 = Matrix(m1)
    sympy_inverse = n <= 25
    sympy_adjugate = n <= 15

    print('Using n={}:'.format(n))
    for row in m1:
        print('│', *row, '│', sep='')

    # Bareiss

    with timeit('Bareiss', 'Determinant'):
        det, adj = adjugate(m1)
        r11 = det

    with timeit('Bareiss', 'Adjugate'):
        det, adj = adjugate(m1)
        r12 = adj

    with timeit('Bareiss', 'Inverse'):
        det, adj = adjugate(m1)
        r13 = [[Fraction(v, det) for v in row] for row in adj]

    # Sympy

    with timeit('Sympy', 'Determinant'):
        r21 = m2.det()

    if sympy_adjugate:
        with timeit('Sympy', 'Adjugate'):
            r22 = m2.adjugate()

    if sympy_inverse:
        with timeit('Sympy', 'Inverse'):
            r23 = m2.inv()

    # Check
    assert r11 == r21
    if sympy_adjugate:
        assert Matrix(r12) == r22
    if sympy_inverse:
        assert Matrix(r13) == r23

    # Determinant
    print('(determinant = ', r11, ')', sep='')


if __name__ == '__main__':
    try:
        n = int(sys.argv[1])
    except:
        benchmark()
    else:
        benchmark(n)

Related: Sympy issue #14840