/*******************************************************
*  LEAST SQUARE APPROXIMATION OF A DISCRETE FUNCTION   *
*  USING ORTHOGONAL POLYNOMIALS                        *
* ---------------------------------------------------- *
* SAMPLE RUN:                                          *
* (Approximate SIN(X) defined by 60 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 [BIBLI 18].                         *
*                                                      *
*                   C++ Release By J-P Moreau, Paris.  *
*******************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  100
#define  KMAX   10


int     I,K,M,NW;
double  X[NMAX],Y[NMAX],SIGMA[NMAX];
double  S[KMAX],ECART[KMAX],ALPHA[KMAX],BETA[KMAX];
double  W1[NMAX],W2[NMAX],W3[NMAX],W4[NMAX];
double  PI,R,W,W0,DW;


    //Auxiliary function
    double PRD(double *X, double *Y, double *Z, int M) {
      double SUM; int I;
      SUM=0.0;
      for (I=1; I<=M; I++)
        SUM += (X[I]*Y[I]/(Z[I]*Z[I]));
      return SUM;
    }


    void MCARRE(int K, int M, double *X, double *Y, double *SIGMA, double *S,
		        double *ALPHA, double *BETA, double *ECART,
				double *P1, double *P2, double *P3, double *P4)  {
/*----------------------------------------------------------------------
!     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:
!     P1..P4: 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.
!---------------------------------------------------------------------*/
      int    I,L;
      double OMEGA,T,W,WPR;
      for (I=1; I<=M; I++) {
        P1[I]=0.0;
        P2[I]=1.0;
      }
      W=(double) 1.0*M;
      BETA[0]=0.0;
      for (I=0; I<K; I++) {
        OMEGA=PRD(Y,P2,SIGMA,M);
        S[I]=OMEGA/W;
        T=PRD(P2,P2,SIGMA,M)/(W*W);
        ECART[I]=sqrt(T);
        for (L=1; L<=M; L++)  P4[L]=X[L]*P2[L];
        ALPHA[I+1]=PRD(P4,P2,SIGMA,M)/W;
        for (L=1; L<=M; L++)
          P3[L]=(X[L]-ALPHA[I+1])*P2[L]-BETA[I]*P1[L];
        WPR=PRD(P3,P3,SIGMA,M);
        if (I+1 <= K-1)  BETA[I+1]=WPR/W;
        W=WPR;
        for (L=1; L<=M; L++) {
          P1[L]=P2[L];
          P2[L]=P3[L];
        }
      }
      OMEGA=PRD(Y,P2,SIGMA,M);
      S[K]=OMEGA/W;
      T=PRD(P2,P2,SIGMA,M)/(W*W);
      ECART[K]=sqrt(T);
    }


    double P(int K, double *S, double *ALPHA, double *BETA, double X) {
/*--------------------------------------------------------------------
!     THIS FUNCTIONE 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.
!-------------------------------------------------------------------*/
      double B,BPR,BSD;
      int I;
      B=S[K];
      BPR=S[K-1]+(X-ALPHA[K])*S[K];
      for (I=K-2; I>=0; I--) {
        BSD=S[I]+(X-ALPHA[I+1])*BPR-BETA[I+1]*B;
        B=BPR;
        BPR=BSD;
      }
      return BPR;
    }


void main()  {

      K=7;
      M=50;
      NW=6;

      PI = 4.0*atan(1.0);

      for (I=1; I<=M; I++) {
        X[I]=PI*(I-1)/(2*(M-1));
        Y[I]=sin(X[I]);
        SIGMA[I]=1.0;
      }

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

      printf("\n   I    COEFFICIENTS    STD DEVIATION\n");
      for (I= 0; I<=K; I++)
        printf("%4d  %12.8f      %12.8f\n", I, S[I], ECART[I]);
      printf("\n   I    VARIABLE          EXACT R         APPROX. R\n");
      W0=0.0;
      DW=PI/10.0;
      for (I=1; I<=NW; I++) {
        W=W0+(I-1)*DW;
        R=P(K,S,ALPHA,BETA,W);
        printf("%4d  %12.8f      %12.8f    %12.8f\n", I, W, sin(W), R);
      }
	  printf("\n");

} // of main program


// end of file approx.cpp