'********************************************************
'*  LEAST SQUARE APPROXIMATION OF A DISCRETE FUNCTION   *
'*  USING ORTHOGONAL POLYNOMIALS                        *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'* (Approximate SIN(X) defined by 50 points, from x=0   *
'*  to x=PI/2, with degree K=7).                        *
'*                                                      *
'*  I     COEFFICIENTS       STD DEVIATION              *
'*  0       0.63383395         0.14142136               *
'*  1       0.66275983         0.30570229               *
'*  2      -0.37764682         0.73926175               *
'*  3      -0.11371820         1.82234391               *
'*  4       0.02861664         4.52658861               *
'*  5       0.00574936        11.29632598               *
'*  6      -0.00096142        28.29696331               *
'*  7      -0.00013770        71.13668785               *
'*                                                      * 
'*  I     VARIABLE       EXACT R     APPROX. R          *
'*  1     0.00000000   0.00000000  -0.00000003          *
'*  2     0.31415927   0.30901699   0.30901698          *
'*  3     0.62831853   0.58778525   0.58778525          *
'*  4     0.94247780   0.80901699   0.80901700          *
'*  5     1.25663706   0.95105652   0.95105651          *
'*  6     1.57079633   1.00000000   0.99999997          *
'*                                                      *
'* ---------------------------------------------------- *   
'* Ref.: From Numath Library By Tuan Dang Trong in      *
'*       Fortran 77 [18].                               *
'*                                                      *
'*                Basic Release By J-P Moreau, Paris.   *
'*                         (www.jpmoreau.fr)            *
'********************************************************
'PRORAM APPROX
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

      M = 50  'NUMBER OF COUPLES X(I), Y(I)
      K = 7   'ORDER OF POLYNOMIALS

DIM X(M), Y(M), SIGMA(M)
DIM S(K), ECART(K), ALPHA(K), BETA(K)
DIM W1(M), W2(M), W3(M), W4(M)

      NW = 6  'NUMBER OF APPROXIMATED POINTS
     
      PI = 4# * ATN(1#)

      FOR I = 1 TO M
        X(I) = PI * (I - 1) / (2 * (M - 1))
        Y(I) = SIN(X(I))
        SIGMA(I) = 1#
      NEXT I

      GOSUB 1000 'Call MCARRE(K,M,X,Y,SIGMA,S,ALPHA,BETA,ECART,W1,W2,W3,W4)

      F$ = "####.########"
      CLS
      PRINT
      PRINT "    I    COEFFICIENTS       STD DEVIATION"
      FOR I = 0 TO K
        PRINT "   "; I; "  ";
        PRINT USING F$; S(I);
        PRINT "      ";
        PRINT USING F$; ECART(I)
      NEXT I
      PRINT
      PRINT "    I    VARIABLE       EXACT R     APPROX. R"
      W0 = 0!
      DW = PI / 10#
      FOR I = 1 TO NW
        W = W0 + (I - 1) * DW
        GOSUB 2000 'R=P(K,S,ALPHA,BETA,W)
        PRINT "   "; I;
        PRINT USING F$; W;
        PRINT USING F$; SIN(W);
        PRINT USING F$; R
      NEXT I
      PRINT

END   'OF MAIN PROGRAM


    'Auxiliary functions
500   SUM = 0#
      FOR II = 1 TO M
        SUM = SUM + Y(II) * W2(II) / SIGMA(II) ^ 2
      NEXT II
      PRD = SUM
RETURN
501   SUM = 0#
      FOR II = 1 TO M
        SUM = SUM + W2(II) * W2(II) / SIGMA(II) ^ 2
      NEXT II
      PRD = SUM
RETURN
502   SUM = 0#
      FOR II = 1 TO M
        SUM = SUM + W4(II) * W2(II) / SIGMA(II) ^ 2
      NEXT II
      PRD = SUM
RETURN
503   SUM = 0#
      FOR II = 1 TO M
        SUM = SUM + W3(II) * W3(II) / SIGMA(II) ^ 2
      NEXT II
      PRD = SUM
RETURN


1000  'Procedure MCARRE(K,M,X,Y,SIGMA,S,ALPHA,BETA,ECART,W1,W2,W3,W4)
'-----------------------------------------------------------------------
'     LEAST SQUARES APPROXIMATION OF A FUNCTION F(X) DEFINED BY M POINTS
'     X(I), Y(I) BY USING ORTHOGONAL POLYNOMIALS
'-----------------------------------------------------------------------
'     INPUTS:
'     K   : DEGREE OF POLYNOMIALS
'     M   : NUMBER OF POINTS
'     X,Y : TABLES OF DIMENSION M TO STORE M ABSCISSAS AND
'           M ORDINATES OF GIVEN POINTS
'     SIGMA : TABLE OF DIMENSION M TO STORE THE STANDARD DEVIATIONS
'             OF VARIABLE Y
'     OUTPUTS:
'     S    : TABLE OF DIMENSION(0:K)
'     ALPHA: TABLE OF DIMENSION (K)
'     BETA : TABLE OF DIMENSION (0:K-1)
'     WORKING SPACE:
'     W1..W4: AUXILIARY TABLES OF DIMENSION (M)
'     NOTE:
'     COEFFICIENTS S,ALPHA,BETA ARE USED TO EVALUATE VALUE
'     AT POINT Z OF THE BEST POLYNOMIAL OF DEGREE K
'     BY USING FUNCTION P(K,S,ALPHA,BETA,Z) DESCRIBED BELOW.
'-----------------------------------------------------------------------
      FOR I = 1 TO M
        W1(I) = 0#
        W2(I) = 1#
      NEXT I
      W = 1! * M
      BETA(0) = 0#
      FOR I = 0 TO K - 1
        GOSUB 500
        OMEGA = PRD
        S(I) = OMEGA / W
        GOSUB 501
        T = PRD / W ^ 2
        ECART(I) = SQR(T)
        FOR L = 1 TO M
          W4(L) = X(L) * W2(L)
        NEXT L
        GOSUB 502
        ALPHA(I + 1) = PRD / W
        FOR L = 1 TO M
          W3(L) = (X(L) - ALPHA(I + 1)) * W2(L) - BETA(I) * W1(L)
        NEXT L
        GOSUB 503
        WPR = PRD
        IF I + 1 <= K - 1 THEN BETA(I + 1) = WPR / W
        W = WPR
        FOR L = 1 TO M
          W1(L) = W2(L)
          W2(L) = W3(L)
        NEXT L
      NEXT I
      GOSUB 500
      OMEGA = PRD
      S(K) = OMEGA / W
      GOSUB 501
      T = PRD / W ^ 2
      ECART(K) = SQR(T)
RETURN


2000  'SUBROUTINE P(K,S,ALPHA,BETA,W)
'---------------------------------------------------------------------
'     THIS SUBROUTINE ALLOWS EVALUATING VALUE AT POINT X OF A FUNCTION
'     F(X), APPROXIMATED BY A SYSTEM OF ORTHOGONAL POLYNOMIALS Pj(X)
'     THE COEFFICIENTS OF WHICH, ALPHA,BETA HAVE BEEN DETERMINED BY
'     LEAST SQUARES (Result in R).
'---------------------------------------------------------------------
      B = S(K)
      BPR = S(K - 1) + (W - ALPHA(K)) * S(K)
      FOR II = K - 2 TO 0 STEP -1
        BSD = S(II) + (W - ALPHA(II + 1)) * BPR - BETA(II + 1) * B
        B = BPR
        BPR = BSD
      NEXT II
      R = BPR
RETURN


'end of file approx.bas