arraybase 0
N = 14: M = 2: Q = 3 # number of points and M.R. polynom degree

dim X(N+1) # 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+1) # 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+1): dim T(N+1) # linear system coefficient
dim A(M+1, Q+1)    # system to be solved
dim p(M+1, Q+1)

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 then 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 rjust(A[i, j], 14.1);
	next j
	print
next i

for j = 0 to M
	for i = j to M
		if A[i, j] <> 0 then exit for
	next i
	if i = M + 1 then print: print "SINGULAR MATRIX '" : end
	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: print "Solutions:"
for i = 0 to M
	print rjust(A[i, M + 1], 16.7);
next i
end