/***************************************************
*  Program to demonstrate Lagrange interpolation   *
*     of Function SIN(X) in double precision       *
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
*                                                  *
*               C++ version by J-P Moreau, Paris   *
*                      (www.jpmoreau.fr)           *
* ------------------------------------------------ *
* SAMPLE RUN:                                      *
*                                                  *
*  Lagrange interpolation of SIN(X):               *
*                                                  *
*    Input X: 0.5                                  *
*    Order of interpolation: 4                     *
*                                                  *
*    SIN(X) = 0.47942552 (exact value: 0.47942554) *
*                                                  *
***************************************************/
#include <stdio.h>
#include <math.h>

// Label: e100

double  XL[10];
double  X[15], Y[15];
double  xx,yy;
int     iv,n;


/*******************************************************
*          Lagrange interpolation subroutine           *
* ---------------------------------------------------- *
* n is the level of the interpolation ( Ex. n=2 is     *
* quadratic ). v is the total number of table values.  *
* X(i), Y(i) are the coordinate table values, Y(i)     *
* being the dependant variable. The X(i) may be arbi-  *
* trarily spaced.  x is the interpolation point which  *
* is assumed to be in the interval  with at least one  *
* table value to the left, and n to the right. If this *
* is violated, n will be set to zero. It is assumed    *
* that the table values are in ascending X(i) order.   *
*******************************************************/
void Interpol_Lagrange()  {
  //Labels: e100,e200,e300,e400
  int i,j,k;
  // Check to see if interpolation point is correct
  if (xx < X[1]) goto e100; 
  if (xx <= X[iv-n]) goto e200;
  // An error has been encountered
e100: n=0 ; return;
  // Find the relevant table interval
e200: i=0;
e300: i++;
  if (xx > X[i]) goto e300;
  i--;
  // Begin interpolation
  for (j=0; j<n+1; j++)  XL[j]=1.0;
  yy=0.0;
  for (k=0; k<n+1; k++) {
    for (j=0; j<n+1; j++) {
      if (j==k) goto e400;
      XL[k]=XL[k]*(xx-X[j+i])/(X[i+k]-X[j+i]);
e400:; }
    yy=yy+XL[k]*Y[i+k];
  }
}


void main()  {
  printf("\n Lagrange interpolation of SIN(X):\n");
  iv=14;
  /* Input sine table
  -----------------------------------------------------------------
   Sine table values from  Handbook of mathematical functions
   by M. Abramowitz and I.A. Stegun, NBS, june 1964
  ---------------------------------------------------------------*/
  X[1]=0.000; Y[1]=0.00000000; X[2]=0.125; Y[2]=0.12467473;
  X[3]=0.217; Y[3]=0.21530095; X[4]=0.299; Y[4]=0.29456472;
  X[5]=0.376; Y[5]=0.36720285; X[6]=0.450; Y[6]=0.43496553;
  X[7]=0.520; Y[7]=0.49688014; X[8]=0.589; Y[8]=0.55552980;
  X[9]=0.656; Y[9]=0.60995199; X[10]=0.721; Y[10]=0.66013615;
  X[11]=0.7853981634; Y[11]=0.7071067812;
  X[12]=0.849; Y[12]=0.75062005; X[13]=0.911; Y[13]=0.79011709;
  X[14]=0.972; Y[14]=0.82601466;
  //---------------------------------------------------------------

  // Input interpolation point
e100: ;

  printf("\n   Input X: "); scanf("%lf",&xx);
  printf("   Order of interpolation: "); scanf("%d",&n);

  Interpol_Lagrange();

  if (n!=0)  {
    printf("\n   SIN(X) = %10.8f  (exact value: %10.8f)\n", yy,sin(xx));
    goto e100;
	printf("\n");
  }
  printf("\n");
}


// End of file Lagrange.cpp