/********************************************************
*    Approximation of a discrete real function F(x) by  *
*    least squares                                      *
* ----------------------------------------------------- *
* Ref.: "Méthodes de calcul numérique, Tome 2 by Claude *
*        Nowakowski, PSI Edition, 1984" [BIBLI 04].     *
* ----------------------------------------------------- *
* 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.069930070                               *
*  A(1) =     3.496503497                               *
*  A(2) =    -3.496503497                               *                      
*  A(3) =     0.000000000                               *
*                                                       *
* 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                               *
*                                                       *
*                    C++ version by J-P Moreau, Paris.  *
*                            (www.jpmoreau.fr)          *
********************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  25

int    i,ij,j,k,n,n1,m,m1,m2;
double C[SIZE][SIZE];
double A[SIZE],B[SIZE],X[SIZE],Xc[SIZE],Y[SIZE],Yx[SIZE];
double p,xx,s,yc;


  double IntPower(double x,int k) {
    if (x==0)  return 0;
    else return (exp(k*log(x)));
  }


void main()  {
  printf("\n Number of points    : "); scanf("%d", &n);
  n--;
  printf("\n Degree of polynomial: "); scanf("%d", &m);
  n1=n+1; m1=m+1; m2=m+2;
  printf("\n Function to approximate:\n");
  for (i=1; i<=n1; i++) {
    printf("  X(%d), Y(%d) = ", i, i); scanf("%lf %lf", &X[i], &Y[i]);
  }
  for (k=1; k<=m2; k++) {
    Xc[k]=0;
    for (i=1; i<=n1; i++)  Xc[k] += IntPower(X[i],k);
  }
  yc=0;
  for (i=1; i<=n1; i++)  yc += Y[i];
  for (k=1; k<=m; k++) {
    Yx[k]=0;
    for (i=1; i<=n1; i++)  Yx[k] += Y[i]*IntPower(X[i],k);
  }
  for (i=1; i<=m1; i++)
	for (j=1; j<=m1; j++) {
      ij=i+j-2;
      if (i==1 && j==1) C[1][1]= n1;
      else C[i][j]=Xc[ij];
    }
  B[1]=yc; for (i=2; i<=m1; i++)  B[i]=Yx[i-1];
  for (k=1; k<=m; k++)
    for (i=k+1; i<=m1; i++) {
      B[i] -= C[i][k]/C[k][k]*B[k];
      for (j=k+1; j<=m1; j++)
        C[i][j] -= C[i][k]/C[k][k]*C[k][j];
    }
  A[m1]=B[m1]/C[m1][m1];
  for (i=m; i>0; i--)  {
    s=0;
    for (k=i+1; k<=m1; k++)  s = s + C[i][k]*A[k];
    A[i] = (B[i]-s)/C[i][i];
  }
  printf("\n Polynomial approximation of degree %d (%d points)\n",m,n+1);
  printf(" Coefficients of polynomial:\n");
  for (i=1; i<=m1; i++)  printf("  A(%d) = %15.9f\n", i-1, A[i]);
  printf("\n Approximated function:\n");
  printf("        X           Y\n");
  for (i=1; i<=n1; i++) {
    xx=X[i]; p=0;
    for (k=1; k<=m1; k++)  p = p*xx + A[m1+1-k];
    printf(" %11.6f %11.6f\n", xx, p);
  }
  printf("\n\n");
}

// end of file approx.cpp