'*********************************************************
'*   Approximation of a discrete real function F(x) by   *
'*   least squares                                       *
'* ----------------------------------------------------- *
'* Ref.: "Methodes de calcul numerique, Tome 2 by Claude *
'*        Nowakowski, PS1 Edition, 1984" [4].            *
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'*                                                       *
'* Number of points    : 11                              *
'*                                                       *
'* Degree of polynomial: 3                               *
'*                                                       *
'* Function to approximate:                              *
'*  X(1), Y(1) = 0,0                                     *
'*  X(2), Y(2) = 0.1,0.2                                 *
'*  X(3), Y(3) = 0.2,0.4                                 *
'*  X(4), Y(4) = 0.3,0.6                                 *
'*  X(5), Y(5) = 0.4,0.8                                 *
'*  X(6), Y(6) = 0.5,1                                   *
'*  X(7), Y(7) = 0.6,0.8                                 *
'*  X(8), Y(8) = 0.7,0.6                                 *
'*  X(9), Y(9) = 0.8,0.4                                 *
'*  X(10), Y(10) = 0.9,0.2                               *
'*  X(11), Y(11) = 1,0                                   *
'*                                                       *
'* Polynomial approximation of degree 3 (11 points)      *
'* Coefficients of polynomial:                           *
'*  A( 0 ) =   -0.069930                                 *
'*  A( 1 ) =    3.496503                                 *
'*  A( 2 ) =   -3.496503                                 *                  
'*  A( 3 ) =   -0.000000                                 *
'*                                                       *
'* Approximated function:                                *
'*      X          Y                                     *
'*   0.000000  -0.069930                                 *
'*   0.100000   0.244755                                 *
'*   0.200000   0.489510                                 *
'*   0.300000   0.664336                                 *
'*   0.400000   0.769231                                 *
'*   0.500000   0.804196                                 *
'*   0.600000   0.769231                                 *
'*   0.700000   0.664336                                 *
'*   0.800000   0.489510                                 *
'*   0.900000   0.244755                                 *
'*   1.000000  -0.069930                                 *
'*                                                       *
'*                  Basic Version By J-P Moreau, Paris.  *
'*                           (www.jpmoreau.fr)           *
'*********************************************************
'Program Approx
DEFINT I-N
DEFDBL A-H, O-Z

NSIZE = 25

DIM C(NSIZE, NSIZE)
DIM A(NSIZE), B(NSIZE), X(NSIZE), Xc(NSIZE), Y(NSIZE), Yx(NSIZE)

  CLS
  PRINT
  INPUT " Number of points    : ", n
  n = n - 1
  PRINT
  INPUT " Degree of polynomial: ", m
  n1 = n + 1: m1 = m + 1: m2 = m + 2
  PRINT
  PRINT " Function to approximate:"
  FOR i = 1 TO n1
    PRINT "  X("; i; "), Y("; i; ") = ";
    INPUT "", X(i), Y(i)
  NEXT i
  FOR k = 1 TO m2
    Xc(k) = 0#
    FOR i = 1 TO n1
      Xc(k) = Xc(k) + X(i) ^ k
    NEXT i
  NEXT k
  yc = 0#
  FOR i = 1 TO n1
    yc = yc + Y(i)
  NEXT i
  FOR k = 1 TO m
    Yx(k) = 0#
    FOR i = 1 TO n1
      Yx(k) = Yx(k) + Y(i) * X(i) ^ k
    NEXT i
  NEXT k
  FOR i = 1 TO m1
    FOR j = 1 TO m1
      ij = i + j - 2
      IF i = 1 AND j = 1 THEN
	C(1, 1) = n1
      ELSE
	C(i, j) = Xc(ij)
      END IF
    NEXT j
  NEXT i
  B(1) = yc
  FOR i = 2 TO m1
    B(i) = Yx(i - 1)
  NEXT i
  FOR k = 1 TO m
    FOR i = k + 1 TO m1
      B(i) = B(i) - C(i, k) / C(k, k) * B(k)
      FOR j = k + 1 TO m1
	C(i, j) = C(i, j) - C(i, k) / C(k, k) * C(k, j)
      NEXT j
    NEXT i
  NEXT k
  A(m1) = B(m1) / C(m1, m1)
  FOR i = m TO 1 STEP -1
    s = 0#
    FOR k = i + 1 TO m1
      s = s + C(i, k) * A(k)
    NEXT k
    A(i) = (B(i) - s) / C(i, i)
  NEXT i
  F$ = "####.######"
  PRINT " Polynomial approximation of degree "; m; " ("; n + 1; " points)"
  PRINT
  PRINT " Coefficients of polynomial:"
  FOR i = 1 TO m1
    PRINT "  A("; i - 1; ") = "; : PRINT USING F$; A(i)
  NEXT i
  PRINT
  PRINT " Approximated function:"
  PRINT "      X          Y"
  FOR i = 1 TO n1
    xx = X(i): p = 0#
    FOR k = 1 TO m1
      p = p * xx + A(m1 + 1 - k)
    NEXT k
    PRINT USING F$; xx; : PRINT USING F$; p
  NEXT i
  PRINT
  PRINT

END

' end of file approx.bas