import numpy as np

rng = np.random.default_rng(0)  #generator
#----------- Input --------------
n = 5 # dimensiunea matricei
A = rng.integers(1,11, size = (n,n))
A = np.triu(A,-1)
A = A.astype('float')
b =  rng.integers(-10,11, size = (n,1)).astype('float')

#-------------------------------------
print('-----------A------------')
print(A)
print('-----------b------------')
print(b)

def EGPP(A: np.ndarray):
    n = len(A)
    p = np.zeros(n, dtype = int)

    for k in range(n - 1):
        p[k] = k
        for i in range(k, n):
            if abs(A[i][k]) is max(abs(A[k : n][k])):
                p[k] = i

        for j in range(k, n):
            temp = A[k][j]
            A[k][j] = A[p[k]][j]
            A[p[k], j] = temp

        for i in range(k + 1, n):
            A[i][k] = A[i][k] / A[k][k]

        for i in range(k + 1, n):
            for j in range(k + 1, n):
                A[i][j] = A[i][j] - A[i][k] * A[k][j]
        return A, p

def S_SUP_TR(U: np.ndarray, b: np.ndarray):
    n = len(U)
    x = np.zeros(n, dtype = float)

    for i in range(n - 1, -1, -1):
        s = b[i]
        for j in range(i + 1, n):
             s = s - U[i][j] * x[j]
        x[i] = s / U[i][i]

    return x

def Solver(A: np.ndarray, b: np.ndarray):
    x = np.zeros(n, dtype = float)
    U, p = EGPP(A)
    for k in range(n - 1):
        b[k], b[p[k]] = b[p[k]], b[k]
        for i in range(k + 1, n):
            b[i] = b[i] - A[i][k] * b[k]
    x = S_SUP_TR(U, b)

    return x

print("Result")
result = Solver(A, b)
print(result)

if result is np.linalg.inv(A) * b:
    print("E corect")
else:
    print("E gresit")