import math
from decimal import Decimal,getcontext

getcontext().prec = 40

def gaussPivoteamento(B,c): # faz eliminação de Gauss com pivoteamento parcial
    A=B 
    b=c
    sol=[]
    for i in range(len(A)):
        pivo = math.fabs(A[i][i])
        linha_pivo = i

        for j in range(i+1,len(A)):
            if math.fabs(A[i][i]) > pivo:
                pivo = math.fabs(A[i][i])
                linha_pivo = j
        
        if linha_pivo != i:
            linhaAuxiliar = A[i]
            A[i] = A[linha_pivo]
            A[linha_pivo] = linhaAuxiliar

            bAuxiliar = b[i]
            b[i] = b[linha_pivo]
            b[linha_pivo] = bAuxiliar

        for m in range(i + 1,len(A)):
            mult = A[m][i] / A[i][i]
            for n in range(i,len(A)):
                A[m][n] -=  mult * A[i][n]
            b[m] -= mult * b[i]
        sol.append(A)
        sol.append(b)
    return (sol)
    

def calculaSolucao(A,b): # Resolve sistema linear triangular superior
    vetorSolucao = []
    for i in range(len(A)):
        vetorSolucao.append(0)
    linha = len(A) - 1
    while linha >= 0:
        xl = b[linha]
        coluna = len(A) - 1
        while coluna > linha:
            xl -= A[linha][coluna]*vetorSolucao[coluna]
            coluna -= 1
        xl /= A[linha][linha]
        linha -= 1
        vetorSolucao[coluna] = xl
    print('Vetor solução:')
    for j in range(len(vetorSolucao)):
        print("a"+str(j)+" = "+str(vetorSolucao[j]))
    return vetorSolucao

# Resolvendo problema particular y=p(x)=a_0+a_1z+...+a_nx^n, com p(s) polinômio interpolador

def crie_matriz(n_linhas, n_colunas, valor): # Criação de matriz genérica com entradas constantes	
	    matriz = [] # lista vazia
	    for i in range(n_linhas):
	        # cria a linha i
	        linha = [] # lista vazia
	        for j in range(n_colunas):
	            linha.append(valor)
	
	        # coloque linha na matriz
	        matriz.append(linha)
	
	    return matriz


def Gram(u):  # Cria matriz para interpolação polinomial
    solucaoDecimal = crie_matriz(len(u),len(u),0)
    for i in range(len(u)):
        mult = u[i]
        solucaoDecimal[i][0] = Decimal(1)
        solucaoDecimal[i][1] = Decimal(mult)
        for j in range(2,len(x)):
            solucaoDecimal[i][j] = Decimal(mult * mult)
            mult = solucaoDecimal[i][j]
    return solucaoDecimal

# Dados para interpolação 
x=[Decimal(1), Decimal(2),Decimal(3)]
y=[Decimal(1), Decimal(3),Decimal(3)]

A=Gram(x)
b=y

sol=gaussPivoteamento(A,b)
print('A=',sol[0],'b=',sol[1])

calculaSolucao(A,b)
print()

def poli(u):
    p=sol[1]
    q=p[0]
    for i in range(1,len(p)):
        q+=p[i]*pow(u,i)
    return q 
    
print(poli(1))