OPTION BASE 0
LET n = 14                        ! number of points and M.R. polynom degree
LET m = 2
LET q = 3
DIM x(0)                          ! data points
MAT REDIM x(n)
DATA 1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 1.73, 1.75, 1.78, 1.80, 1.83
FOR c = LBOUND(x) TO UBOUND(x)
    READ x(c)
NEXT c
DIM y(0)                          ! data points
MAT REDIM y(n)
DATA 52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 66.28, 68.10, 69.92, 72.19, 74.46
FOR c = LBOUND(y) TO UBOUND(y)
    READ y(c)
NEXT c
DIM s(0)                          ! linear system coefficient
MAT REDIM s(n)
DIM t(0)
MAT REDIM t(n)
DIM a(0,0)                        ! system to be solved
MAT REDIM a(m, q)
DIM p(0,0)
MAT REDIM p(m, q)

FOR k = 0 TO 2*m
    LET s(k) = 0
    LET t(k) = 0
    FOR i = 0 TO n
        LET s(k) = s(k)+x(i)^k
        IF k <= m THEN LET 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
        LET a(fila, columna) = s(fila+columna)
    NEXT columna
    LET a(fila, columna) = t(fila)
NEXT fila
PRINT "Linear system coefficents:"
FOR i = 0 TO m
    FOR j = 0 TO m+1
        PRINT  USING "######.#": a(i, j);
    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 '"
       STOP
    END IF
    FOR k = 0 TO m+1
        LET p(j, k) = a(i, k)
        LET a(i, k) = p(j, k)
        LET a(j, k) = a(i, k)
    NEXT k
    LET z = 1/a(j, j)
    FOR k = 0 TO m+1
        LET a(j, k) = z*a(j, k)
    NEXT k
    FOR i = 0 TO m
        IF i <> j THEN
           LET z = -a(i, j)
           FOR k = 0 TO m+1
               LET 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  USING "  #####.#######": a(i, m+1);
NEXT i
END