import numpy as np


def gauss_jordan_elimination(matrix, constants):
    n = len(matrix)
    matrix = matrix.astype(float)  # Převést matici na float64
    constants = constants.astype(float)  # Převést vektor na float64
    for i in range(n):
        # Normalizace
        div = matrix[i, i]
        if div == 0:
            raise ValueError("Nulový prvek na diagonále, nelze pokračovat.")
        matrix[i] /= div
        constants[i] /= div

        # Eliminace
        for k in range(n):
            if k != i:
                m = matrix[k, i]
                matrix[k] -= m * matrix[i]
                constants[k] -= m * constants[i]

    return constants


# Příklad použití
M = np.array([[2, 5, 4, 1, 20],
              [1, 3, 2, 1, 11],
              [2, 10, 9, 7, 40],
              [3, 8, 9, 2, 37]
              ])

# split matrix into two parts
A = M[:, :-1]
B = M[:, -1]

solution = gauss_jordan_elimination(A, B)
print("Řešení soustavy je:")
for i, x in enumerate(solution):
    if x == 0:
        x = abs(x)
    print(f"x{i + 1} = {x:.1f}".rstrip('0').rstrip('.'))