import numpy as np
from fractions import Fraction

M = np.array([[2, 5, 4, 1, 20],
              [1, 3, 2, 1, 11],
              [2, 10, 9, 7, 40],
              [3, 8, 9, 2, 37]
              ])

# rozdělení soustavy na matici koeficientů a konstant
A = M[:, :-1]
B = M[:, -1]

# hlavní determinant
det_A = np.round(np.linalg.det(A))

# determinanty pro x1, x2, ...
for column_idx in range(A.shape[1]):
    C = A.copy()
    C[:, column_idx] = B
    new_det = np.round(np.linalg.det(C))

    # dostat x1, x2, ...
    answer = new_det / det_A
    answer = 0.0 if answer == -0.0 else answer

    if isinstance(answer, float):
        answer = Fraction(answer).limit_denominator()
    print("Řešení soustavy je:")
    print(f'x{column_idx+1} = {answer}')