// number of points and M.R. polynom degree
N = 14
M = 2
Q = 3

dim X(N)            // data points
X(0) = 1.47 : X(1) = 1.50 : X(2) = 1.52
X(3) = 1.55 : X(4) = 1.57 : X(5) = 1.60
X(6) = 1.63 : X(7) = 1.65 : X(8) = 1.68
X(9) = 1.70 : X(10) = 1.73 : X(11) = 1.75
X(12) = 1.78 : X(13) = 1.80 : X(14) = 1.83
dim Y(N)            // data points
Y(0) = 52.21 : Y(1) = 53.12 : Y(2) = 54.48
Y(3) = 55.84 : Y(4) = 57.20 : Y(5) = 58.57
Y(6) = 59.93 : Y(7) = 61.29 : Y(8) = 63.11
Y(9) = 64.47 : Y(10) = 66.28 : Y(11) = 68.10
Y(12) = 69.92 : Y(13) = 72.19 : Y(14) = 74.46

dim S(N), T(N)        // linear system coefficient
dim A(M, Q)           // sistem to be solved
dim p(M, Q)

for k = 0 to 2*M
    S(k) = 0 : T(k) = 0
    for i = 0 to N
        S(k) = S(k) + X(i) ^ k
        if k <= M  T(k) = T(k) + Y(i) * X(i) ^ k
    next i
next k

// build linear system
for fila = 0 to M
    for columna = 0 to M
        A(fila, columna) = S(fila+columna)
    next columna
    A(fila, columna) = T(fila)
next fila

print "Linear system coefficents:"
for i = 0 to M
    for j = 0 to M+1
        print A(i,j) using "#####.#";
    next j
    print
next i

for j = 0 to M
    for i = j to M
        if A(i,j) <> 0  break
    next i
    if i = M+1 then
        print "\nSINGULAR MATRIX '"
        end
    end if
    for k = 0 to M+1
		p(j,k) = A(i,k) : A(i,k) = p(j,k) : A(j,k) = A(i,k)
    next k
    z = 1 / A(j,j)
    for k = 0 to M+1
        A(j,k) = z * A(j,k)
    next k
    for i = 0 to M
        if i <> j then
            z = -A(i,j)
            for k = 0 to M+1
                A(i,k) = A(i,k) + z * A(j,k)
            next k
        end if
    next i
next j

print "\nSolutions:"
for i = 0 to M
    print A(i,M+1) using "#######.#######";
next i
print
end