import math
from decimal import Decimal,getcontext

getcontext().prec = 40

def gaussPivoteamento(A,b):

    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]
    calculaSolucao(A,b)

vetorSolucao = []
solucFinal = []

def calculaSolucao(A,b):
    
    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
    
    for j in range(len(vetorSolucao)):
        print("a"+str(j)+" = "+str(vetorSolucao[j]))

def crie_matriz(n_linhas, n_colunas, valor):	
	    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 solucao_final(u):
    soma = 0
    for i in range(len(u)):
        for j in range(len(u)):
            aux = math.pow(u[i],j)
            soma += Decimal(vetorSolucao[j]) * Decimal(aux)
        
        solucFinal.append(soma)
        soma =0

    for j in range(0,len(u)):
        print("P("+str(j)+") = "+str(solucFinal[j]))
        
        
# Resolvendo problema particular y=p(x)=a_0+a_1z+...+a_nx^n, com p(s) polinômio interpolador

x=[Decimal(1), Decimal(2),Decimal(3)]
y=[Decimal(1), Decimal(3),Decimal(3)]

solucaoDecimal = crie_matriz(len(x),len(x),0)

for i in range(len(x)):
    mult = x[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]

gaussPivoteamento(solucaoDecimal,y)
print()
solucao_final(x)

print(y)