/* Exercice V, feuille d'exercices 3 */

/* Note : il ne s'agit pas de la correction de la routine a rendre */
/* en raison de la structure de donnee utilisee pour les polynomes */
/* Cliquez ici pour acceder a la correction de la routine a rendre */


#include <stdio.h>
#include <assert.h>

#define alloc(a)      ((a*)malloc(sizeof(a)))
#define alloctab(a,n) ((a*)malloc((n)*sizeof(a)))


/* A polynomial of degre n is a n+2 array of doubles
   The first element represents de degree of the polynomial
   If the polynomial is 0, the degree will be 0 */


int CheckPolynomial (double ** Polynomial_Ptr)
     /* Get rid of possible zeros */
{
  int degre;

  degre=*Polynomial_Ptr[0];

  while (((*Polynomial_Ptr)[degre+1]==0)&&degre>0)
    {
      (*Polynomial_Ptr[0])--;
      free((*Polynomial_Ptr)[degre+1]);
      degre--;
    }

}
	

int GetPolynomial (double ** Polynomial_Ptr)
{
  int i;
  int degre;

  printf("Entrez le degre du polynome : ");
  scanf("%d",&degre);

  /* The degree is known, allocate the array accordingly */

  *Polynomial_Ptr = alloctab(double,degre+2);

  (*Polynomial_Ptr)[0]=degre;

  for (i=0 ; i<=degre ; i++)
    {
      if (i==0)
	printf("Donner la constante : ");
      else if (i==1)
	printf("Donner la constante devant x : ");
      else printf("Donner le coefficient devant x^%d : ",i);
      scanf("%lf",&((*Polynomial_Ptr)[i+1])); 
    }

  assert(i=degre+1);

  CheckPolynomial (Polynomial_Ptr);
}


int PrintPolynomial (double * Polynomial)
{
  int i;
  int degre;

  degre = Polynomial[0];

  if (degre>=1) 
    for (i=degre ; i>=1 ; i--)
      printf("%lf x^%d + ",Polynomial[i+1],i);
  printf("%lf\n",Polynomial[1]);
}


int Degree (double * Polynomial)
{
  return Polynomial[0];
}


int Valuation (double * Polynomial)
{
  int i;
  int degre;

  degre=Polynomial[0];

  i=0;
  while ((Polynomial[i+1]==0)&&(i<=degre)) 
    i++;

  return i;
}


int AddPolynomial (double * P1 ,
		   double * P2 ,
		   double ** P3_Ptr)
/* P3_Ptr=P1+P2 */
{
  int i;
  int degre1, degre2, degre3, max, min;
  double * largerP;
  
  degre1 = P1[0];
  degre2 = P2[0];

  if (degre1>degre2)
    {
      max=degre1;
      min=degre2;
      largerP=P1;
    }
  else
    {
      max=degre2;
      min=degre1;
      largerP=P2;
    }

  degre3=max;
  *P3_Ptr = alloctab(double,degre3+2);
  (*P3_Ptr)[0] = degre3;

  for (i=0 ; i<=min ; i++)
    (*P3_Ptr)[i+1]=P1[i+1]+P2[i+1];

  while (i<=max)
    {
      (*P3_Ptr)[i+1]=largerP[i+1];
      i++;
    }

  CheckPolynomial(P3_Ptr);
}


int ExtMultPolynomial(double * P1 ,
		      double lambda ,
		      double ** P3_Ptr)
     /* P3 = lambda P1 */
{
  int i;
  int degre;


  degre=P1[0];
  *P3_Ptr = alloctab(double,degre+2);
  (*P3_Ptr)[0] = degre;

  for (i=0 ; i<=degre ; i++)
    (*P3_Ptr)[i+1]=lambda*P1[i+1];

  CheckPolynomial(P3_Ptr);
}


int MonMultPolynomial(double * P1 ,
		      int n ,
		      double ** P3_Ptr)
     /* P3 = x^n P1 */
{
  int i;
  int degre;

  degre=P1[0];

  if ((degre==0)&&(P1[1]==0))
    {
      *P3_Ptr = alloctab(double,2);
      *P3_Ptr[0]=0;
      *P3_Ptr[1]=0;
    }
  else
    {
      *P3_Ptr = alloctab(double,degre+2+n);
      (*P3_Ptr)[0] = degre+n;
      
      for (i=0 ; i<n ; i++)
	(*P3_Ptr)[i+1]=0;
      for (i=n ; i<=n+degre ; i++)
	(*P3_Ptr)[i+1]=P1[i-n+1];
    }
}


int IntMultPolynomial(double * P1 ,
		   double * P2 ,
		   double ** P3_Ptr)
/* P3_Ptr=P1 P2 */
{
  int i;
  int degre1;
  double * Auxi0;
  double * Auxi1;
  double * Auxi2;
  double * Auxi3;

  degre1=P1[0];

  Auxi0 = alloctab(double,2);
  Auxi0[0]=0;
  Auxi0[1]=0;

  for (i=0 ; i<=degre1 ; i++)
    {
      MonMultPolynomial(P2,i,&Auxi1);
      ExtMultPolynomial(Auxi1,P1[i+1],&Auxi2);
      AddPolynomial(Auxi0,Auxi2,&Auxi3);
      free(Auxi0);
      free(Auxi1);
      free(Auxi2);
      Auxi0=Auxi3;
    }

  *P3_Ptr=Auxi0;
}


double EvalPolynomial(double * P1 ,
		      double x)
     /* returns P1(x) */
{
  int i;
  int degre;
  double res;

  degre = P1[0];

  res = P1[degre+1];

  for (i=degre-1 ; i>=0 ; i--)
      res = res*x+P1[i+1];
  
  return res;
}


int GetPointsXY (int * NumberOfPoints_Ptr ,
		 double ** ListOfPointsX_PPtr ,
		 double ** ListOfPointsY_PPtr )
{
  int i;

  printf("Donner le nombre de points a interpoler : ");
  scanf("%d",NumberOfPoints_Ptr);

  /* The number of points is known, allocate the array accordingly */
  *ListOfPointsX_PPtr = alloctab(double,*NumberOfPoints_Ptr);
  *ListOfPointsY_PPtr = alloctab(double,*NumberOfPoints_Ptr);

  for (i=1 ; i<=*NumberOfPoints_Ptr ; i++)
    {
      printf("Donner la valeur X du point %d : ",i);
      scanf("%lf",&((*ListOfPointsX_PPtr)[i-1])); 
      printf("Donner la valeur Y du point %d : ",i);
      scanf("%lf",&((*ListOfPointsY_PPtr)[i-1])); 
    }

  assert(i=*NumberOfPoints_Ptr+1);
  printf("%d points enregistres\n",*NumberOfPoints_Ptr);
}


int Aitken (int NumberOfPoints ,
	    double * ListOfPointsX_Ptr, 
	    double * ListOfPointsY_Ptr,
	    double ** Polynomial_Ptr)
{
  int i, j, k;
  double ** P;
  double ** Q;
  double * P1;
  double * P2;
  double * Auxi1;
  double * Auxi2;
  double * Auxi3;
  double * Auxi4;

  k=NumberOfPoints-1;

  P=alloctab(double*,NumberOfPoints);
  
  for (i=0 ; i<=k ; i++)
    {
      P[i]=alloctab(double,2);
      (P[i])[0]=0;
      (P[i])[1]=ListOfPointsY_Ptr[i];
    }
  
  for (j=1 ; j<=k ; j++)
    {
      Q=alloctab(double*,k-j+1);
      for (i=0 ; i<=k-j ; i++)
	{
	  /* P for superindex j-1, Q for superindex j */

	  P1=alloctab(double,3);
	  P1[0]=1;
	  P1[1]=ListOfPointsX_Ptr[i+j];
	  P1[2]=-1;
	  
	  P2=alloctab(double,3);
	  P2[0]=1;
	  P2[1]=-ListOfPointsX_Ptr[i];
	  P2[2]=+1;

	  IntMultPolynomial(P1,P[i],&Auxi1);
	  IntMultPolynomial(P2,P[i+1],&Auxi2);
	  AddPolynomial(Auxi1,Auxi2,&Auxi3);
	  ExtMultPolynomial(Auxi3,1/(ListOfPointsX_Ptr[i+j]-ListOfPointsX_Ptr[i]),&Auxi4);
	  
	  free(Auxi1);
	  free(Auxi2);
	  free(Auxi3);
	  free(P1);
	  free(P2);
	  
	  Q[i]=Auxi4;
	}
      free(P);
      P=Q;
    }
  
  *Polynomial_Ptr=P[0];
}


int main()
{
  int Number;
  double * Polynomial;
  double * Xs;
  double * Ys;

  GetPointsXY(&Number,&Xs,&Ys);
  Aitken(Number,Xs,Ys,&Polynomial);
  PrintPolynomial(Polynomial);
}