/***************************************************
*   Program to demonstrate Newton 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:                                      *
*                                                  *
*  Newton interpolation of SIN(X):                 *
*                                                  *
*    Input X: 0.5                                  *
*                                                  *
*    SIN(X) = 0.47961461 (exact value: 0.47942554) *
*    Error estimate: -0.00000028                   *
*                                                  *    
***************************************************/
#include <stdio.h>
#include <math.h>

//Label: 100

#define  SIZE  15

double  X[SIZE],Y[SIZE],Y1[SIZE],Y2[SIZE],Y3[SIZE];
double  e,xx,yy;
int     iv, n;


/*******************************************************
* Newton divided differences interpolation subroutine  *
* ---------------------------------------------------- *
* Calculates cubic interpolation for a given table.    *
* iv 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 arbitrarily spa- *
* ced. xx is the interpolation point which is assumed  *
* to be in the interval  with at least one table value *
* to the left, and 3 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. e is the error   *
* estimate. The returned value is yy.                  *
*******************************************************/
void Interpol_Newton()  {
//Labels: e100,e200,e300
int i; double a,b,c;
  // Check to see if interpolation point is correct
  n=1;
  if (xx>=X[1]) goto e100;
  n=0 ; return;
e100: if (xx<=X[iv-3]) goto e200;
  n=0 ; return;
  // Generate divided differences
e200: for (i=1 ; i<iv; i++)
    Y1[i]=(Y[i+1]-Y[i])/(X[i+1]-X[i]);
  for (i=1; i<iv-1; i++)
    Y2[i]=(Y1[i+1]-Y1[i])/(X[i+1]-X[i]);
  for (i=1; i<iv-2; i++)
    Y3[i]=(Y2[i+1]-Y2[i])/(X[i+1]-X[i]);
  // Find relevant table interval
  i=0;
e300: i++;
  if (xx>X[i]) goto e300;
  i--;
  // Begin interpolation
  a=xx-X[i];
  b=a*(xx-X[i+1]);
  c=b*(xx-X[i+2]);
  yy=Y[i]+a*Y1[i]+b*Y2[i]+c*Y3[i];
  //Calculate next term in the expansion for an error estimate
  e=c*(xx-X[i+3])*yy/24.0;
}


void main()  {
  iv=14;
  printf("\n Newton interpolation of SIN(X):\n");
  /* 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");

  printf("   Input X: "); scanf("%lf",&xx);

  Interpol_Newton();

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

}

// End of file Newton.cpp