/*******************************************************************
*       BEST APPROXIMATION OF A DISCRETE REAL FUNCTION F(X)        *
*             BY STIEFEL-REMES'S POLYNOMIAL METHOD                 *
* ---------------------------------------------------------------- *
* SAMPLE RUN:                                                      *
* (Approximate Exp(x) from 0.0 to 1.0).                            *
*                                                                  *
* AFTER 7 ITERATIONS                                               *
* POLYNOMIAL COEFFICIENTS ARE:                                     *
*   0   0.99999891                                                 *
*   1   1.00008075                                                 *
*   2   0.49908690                                                 *
*   3   0.17042178                                                 *
*   4   0.03478373                                                 *
*   5   0.01390866                                                 *
*           X          Y         YAPP       DIFF                   *
*   1   0.000000   1.000000   0.999999  -0.00000109                *   
*   2   0.100000   1.105171   1.105172   0.00000097                *
*   3   0.200000   1.221403   1.221402  -0.00000074                * 
*   4   0.300000   1.349859   1.349858  -0.00000092                *
*   5   0.400000   1.491825   1.491825   0.00000030                *
*   6   0.500000   1.648721   1.648722   0.00000109                *
*   7   0.600000   1.822119   1.822119   0.00000046                *
*   8   0.700000   2.013753   2.013752  -0.00000082                *
*   9   0.800000   2.225541   2.225540  -0.00000084                * 
*  10   0.900000   2.459603   2.459604   0.00000088                *
*  11   1.000000   2.718282   2.718281  -0.00000109                *
*                                                                  *
* ---------------------------------------------------------------- *
* Reference: From Numath Library By Tuan Dang Trong in Fortran 77  *
*            [BIBLI 18].                                           * 
*                                                                  *
*                               C++ Release By J-P Moreau, Paris   *
*                                      (www.jpmoreau.fr)           *
*******************************************************************/ 
#include <stdio.h>
#include <math.h>

#define  NMAX  100


  double X[NMAX],Y[NMAX],A[NMAX];
  int    I,J,M,N,N2,NZ,IT,ITMAX,NW;
  double V,Z,Z0,DZ;


    void CHEFIT(int M, int N, double *X, double *Y, double *A, double *T, 
	            double *AX, double *AY, double *AH, double *BY, double *BH, 
	  		    int *IN0, int *IT, int ITMAX);

    void CHEAPP(int M, int N, double *X, double *Y, int ITMAX,
	                 int *IT, double *A)  {
/*==================================================================
!     BEST APPROXIMATION OF A CONTINUOUS FUNCTION WITH REAL VARIABLE
!     IN THE SENSE OF MAXIMUM NORM.
!     THE STIEFEL-REMES'S METHOD IS USED.
!----- CALLING MODE:
!          CHEAPP(M,N,X,Y,ITMAX,IT,A,W,IW);
!----- INPUTS:
!     M:  NUMBER OF POINTS OF THE FUNCTION
!     N:  DEGREE OF APPROXIMATION POLYNOMIAL
!     X:  TABLE OF SIZE (M) STORING THE ABSCISSAS
!     Y:  TABLE OF SIZE (M) STORING THE ORDINATES
!     ITMAX: MAXIMUM NUMBER OF ITERATIONS
!     IT:    ACTUAL NUMBER OF ITERATIONS MADE
!----- OUTPUTS:
!     A:  TABLE OF SIZE (0:N) STORING THE COEFFICIENTS
!         OF APPROXIMATION POLYNOMIAL
!----- WORKING ZONES:
!     W : TABLE OF SIZE (M+5N+10)
!     IW: TABLE OF INTEGERS OF SIZE (N+2)
!----- NOTE:
!     ABSCISSAS X(I) MUST BE IN INCREASING ORDER 
!----- REFERENCE:
!     PROCEDURES ALGOL EN ANALYSE NUMERIQUE
!     EDITIONS DU CNRS, 1967
!==================================================================*/
      double W1[NMAX],W2[NMAX],W3[NMAX],W4[NMAX],W5[NMAX],W6[NMAX];
      int IW[NMAX];

      CHEFIT(M,N,X,Y,A,W1,W2,W3,W4,W5,W6,IW,IT,ITMAX);

	}

    void CHEFIT(int M, int N, double *X, double *Y, double *A, double *T, 
	            double *AX, double *AY, double *AH, double *BY, double *BH, 
				int *IN0, int *IT, int ITMAX)  {
      int I,IMAX,J,K,L;
      double DX,H,P,TMAX,W;
      K=(M-1) / (N+1);
      for (I=1; I<=N+1; I++)  IN0[I]=(I-1)*K+1;
      IN0[N+2]=M;
      *IT=0;
e1:   *IT=*IT+1;
      for (I=1; I<=N+2; I++) {
        AX[I]=X[IN0[I]];
        AY[I]=Y[IN0[I]];
        AH[I]=pow(-1.0,I-1);
      }
      for (I=2; I<=N+2; I++) {
		for (J=I-1; J<=N+2; J++) {
          BY[J]=AY[J];
          BH[J]=AH[J];
        }
        for (J=I; J<=N+2; J++) {
          DX=AX[J]-AX[J-I+1];
		  if (fabs(DX)>1e-12) {
            AY[J]=(BY[J]-BY[J-1])/DX;
            AH[J]=(BH[J]-BH[J-1])/DX;
          }
		  else {
            printf(" CHEFIT Divide by zero error.\n");
			return;
          }
        }
      }
      H=-AY[N+2]/AH[N+2];
      for (I=0; I<=N; I++) {
        A[I]=AY[I+1]+AH[I+1]*H;
        BY[I+1]=0.0;
      }
      BY[1]=1.0;
      TMAX=fabs(H);
      IMAX=IN0[1];
      for (I=1; I<=N; I++)
		for (J=0; J<I; J++) {
          W=BY[I+1-J]-BY[I-J]*X[IN0[I]];
          A[J]=A[J]+A[I]*W;
          BY[I+1-J]=W;
        }
	  for (I=1; I<=M; I++) {
        T[I]=A[N];
        for (J=1; J<=N; J++)
          T[I]=T[I]*X[I]+A[N-J];
        T[I]=T[I]-Y[I];
        if (fabs(T[I]) > TMAX) {
          TMAX=fabs(T[I]);
          IMAX=I;
        }
      }
      for (I=1; I<=N+2; I++) {
        L=I;
        if (IMAX < IN0[L]) goto e2;
        if (IMAX == IN0[L] || *IT >= ITMAX) return; //normal return
      }
      L=N+2;
e2:   P=T[IMAX]*T[IN0[L]];
      if (P < 0.0) goto e3;
      IN0[L]=IMAX;
      goto e1;

e3:   if (IN0[1] < IMAX) goto e4;
      for (I=1; I<=N+1; I++)
        IN0[N+3-I]=IN0[N+2-I];
      IN0[1]=IMAX;
      goto e1;

e4:   if (IN0[N+2] <= IMAX) goto e5;
      IN0[L-1]=IMAX;
      goto e1;

e5:   for (I=1; I<=N+1; I++)  IN0[I]=IN0[I+1];
      IN0[N+2]=IMAX;
      goto e1;
	}


void main()  {

  M=21;         //NUMBER OF POINTS OF THE FUNCTION
  N=5;          //DEGREE OF APPROXIMATION POLYNOMIAL
  NZ=11;        //NUMBER OF POINTS FOR VERIFICATION
  ITMAX=20;     //MAXIMUM NUMBER OF ITERATIONS

/*  X: TABLE OF SIZE (M) STORING THE ABSCISSAS
    Y: TABLE OF SIZE (M) STORING THE ORDINATES
    A: TABLE OF SIZE (0:N) STORING THE COEFFICIENTS
       OF APPROXIMATION POLYNOMIAL
  
  DEFINE CURVE (M PTS)  */
  for (I=1; I<=M; I++) {
    X[I]=1.0*(I-1)/(1.0*(M-1));
    Y[I]=exp(X[I]);
  }

  CHEAPP(M,N,X,Y,ITMAX,&IT,A);

//PRINT RESULTS
  printf("\n AFTER %d ITERATIONS\n", IT);
  printf(" POLYNOMIAL COEFFICIENTS ARE:\n");
  for (I=0; I<=N; I++) printf("%4d%13.8f\n", I, A[I]);

//VERIFICATION (NZ PTS)
  printf("           X          Y         YAPP       DIFF\n");
  Z0=0.0;
  DZ=0.1;
  for (I=1; I<=NZ; I++) {
    Z=Z0+(I-1)*DZ;
    V=A[N];
    for (J=N-1; J>-1; J--) V=A[J]+V*Z;
    printf("%4d%11.6f%11.6f%11.6f%13.8f\n", I,Z,exp(Z),V,V-exp(Z));
  }
  printf("\n");
}

// end of file tcheapp.cpp