The following algorithm solves the equation Ax=b.

import math
 
 
def factorize(A):
    n = int(math.sqrt(len(A)))
    diagA = True
    L = []
 
    # check if no diagonal entries are zero
    for i in range(0, len(A), n + 1):
        if A[i] == 0:
            diagA = False
 
    # factorize with no row exchanges
    if diagA:
        for i in range(len(A)):
            r, c = i2rc(i, n)
 
            if r == c:  # diagonal entries
                L.append(1)
                continue
            elif r > c:  # lower triangular entries
                ml = A[i] / A[n * (c - 1) + c - 1]
                L.append(ml)
                r0 = A[n * (r - 1):n * r]
                r1 = A[n * (c - 1):n * c]
                A[n * (r - 1):n * r] = rewrite(r0, r1, ml)
                continue
            elif r < c:  # upper triangular entries
                L.append(0)
                continue
 
    U = A
 
    return (L, U)
 
 
def rewrite(r0, r1, ml):
    for i in range(len(r0)):
        r0[i] = r0[i] - ml * r1[i]
 
    return r0
 
 
def rc2i(r, c, n):
    i = n * (r - 1) + c - 1
 
    return i
 
 
def i2rc(i, n):
    r = i // n + 1
    c = i % n + 1
 
    return (r, c)
 
 
def solv_lw_tri(L, b):
    n = len(b)
    c = []
 
    for i in range(n):
        c.append(b[i])
 
        for j in range(i):
            c[i] = c[i] - L[rc2i(i + 1, j + 1, n)] * c[j]
 
        c[i] = c[i] / L[rc2i(i + 1, i + 1, n)]
 
    return c
 
 
def solv_up_tri(U, c):
    n = len(c)
    x = [0] * n
 
    for i in range(n - 1, -1, -1):
        x[i] = c[i]
 
        for j in range(i + 1, n):
            x[i] = x[i] - U[rc2i(i + 1, j + 1, n)] * x[j]
 
        x[i] = x[i] / U[rc2i(i + 1, i + 1, n)]
 
    return x
 
 
class Test(object):
    """Test methods"""
    A = [2, 1, 3, 4, 3, 10, -2, 1, 7]
    b = [11, 28, 3]
    L = [1, 0, 0, 2, 1, 0, -1, 2, 1]
    U = [2, 1, 3, 0, 1, 4, 0, 0, 2]
    c = [11, 6, 2]
    x = [3, 2, 1]
 
    def test_factorize(self, A):
        assert factorize(A) == (self.L, self.U)
 
    def test_solv_lw_tri(self, L, b):
        assert solv_lw_tri(L, b) == self.c
 
    def test_solv_up_tri(self, U, c):
        assert solv_up_tri(U, c) == self.x
 
 
if __name__ == "__main__":
    test = Test()
    test.test_factorize(test.A)
    test.test_solv_lw_tri(test.L, test.b)
    test.test_solv_up_tri(test.U, test.c)

import numpy as np
 
def polyfit(x, fx, w, n):
    """Calculate the coeffients of the polynomial that best approximate (x, y).
 
    :param x (list): indipendent variable
    :param fx (list): dipendent variable
    :param w (list): weight of each node
    :param n (int): degree of the polynomial
    :return a (list): coeffientes of the polynomial
    """
    A = []
    for i in range(len(x)):  # rows
        r = []
        for j in range(n, 0, -1):  # columns
            r.append(w[i]**.5 * x[i]**(n-j))
        A.append(r)
 
    b = []
    for i in range(len(x)):  # columns
        b.append([w[i]**.5 * fx[i]])
 
    A = np.array(A)
    b = np.array(b)
 
    a = np.linalg.lstsq(np.dot(A.T, A), np.dot(A.T, b))  # (A' * A) \ (A' * b)
    a = a[0].T.tolist()[0]
 
    return a