PROGRAM MULTIPLE_REGRESSION

!$DOUBLE

CONST N=14,M=2,Q=3 ! number of points and M.R. polynom degree

DIM X[N],Y[N]      ! data points
DIM S[N],T[N]      ! linear system coefficient
DIM A[M,Q]         ! sistem to be solved

BEGIN

   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)
   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 I%=0 TO N DO
     READ(X[I%])
   END FOR

   FOR I%=0 TO N DO
     READ(Y[I%])
   END FOR

   FOR K%=0 TO 2*M DO
      S[K%]=0  T[K%]=0
      FOR I%=0 TO N DO
         S[K%]=S[K%]+X[I%]^K%
         IF K%<=M THEN T[K%]=T[K%]+Y[I%]*X[I%]^K% END IF
      END FOR
   END FOR

! build linear system

   FOR ROW%=0 TO M DO
     FOR COL%=0 TO M DO
       A[ROW%,COL%]=S[ROW%+COL%]
     END FOR
     A[ROW%,COL%]=T[ROW%]
   END FOR

   PRINT("LINEAR SYSTEM COEFFICENTS") PRINT
   FOR I%=0 TO M DO
     FOR J%=0 TO M+1 DO
        WRITE(" ######.#";A[I%,J%];)
     END FOR
     PRINT
   END FOR
   PRINT

   FOR J%=0 TO M DO
         FOR I%=J% TO M DO
              EXIT IF A[I%,J%]<>0
         END FOR
         IF I%=M+1 THEN
             PRINT("SINGULAR MATRIX !")
             !$STOP
         END IF
         FOR K%=0 TO M+1 DO
             SWAP(A[J%,K%],A[I%,K%])
         END FOR
         Y=1/A[J%,J%]
         FOR K%=0 TO M+1 DO
             A[J%,K%]=Y*A[J%,K%]
         END FOR
         FOR I%=0 TO M DO
             IF I%<>J% THEN
                 Y=-A[I%,J%]
                 FOR K%=0 TO M+1 DO
                    A[I%,K%]=A[I%,K%]+Y*A[J%,K%]
                 END FOR
             END IF
         END FOR
   END FOR
   PRINT

   PRINT("SOLUTIONS") PRINT
   FOR I%=0 TO M DO
      PRINT("c";I%;"=";)
      WRITE("#####.#######";A[I%,M+1])
   END FOR

END PROGRAM