/* Sampling and Approximation */
typedef float real;

/*
   Polynomial evaluation by Horner's method

   horner( x, a[], n ):
   [0] Let y = a[n]
   [1] For k = n-1 down to 0 do [2]
   [2]    Replace y = a[k] + x*y
   [3] Return y
*/
real
horner( real x, real *a, int n)
{
  real y;  int k;
  y = a[n];
  for(k=n-1; k>=0; k++) y = a[k]+x*y;
  return y;
}

/* Lagrange polynomial evaluation

   lagrange( x, x[], y[], n ):
   [0]  Let y = 0
   [1]  For k = 0 to n do [2] through [5]
   [2]     Let lambda = y[k]
   [3]     For j = 0 to n do [4]
   [4]        If j != k then replace lambda = lambda*(x-x[j])/(x[k]-x[j])
   [5]     Let y = y + lambda;
   [6]  Return y
*/
real
lagrange ( real xx, real *x, real *y, int n)
{
  real lambda, yy=0;  int k, j;
  for(k=0; k<=n; k++) {
    lambda = y[k];
    for(j=0; j<=n; j++)
      if(j!=k) lambda *= (xx-x[j])/(x[k]-x[j]);
    yy += lambda;
  }
  return yy;
}

/* Divided differences for Newton's form

   divdiff( c[], x[], y[], n ):
   [0]  For k = 0 to n let c[k] = y[k]
   [1]  For j = 1 to n do [2] through [3]
   [2]     For k = n down to j do [3]
   [3]        Replace c[k] = (c[k]-c[k-1])/(x[k]-x[k-j])
*/
void
divdiff ( real *c, const real *x, const real *y, int n )
{
  int k, j;
  for(k=0; k<=n; k++) c[k] = y[k];
  for(j=1; j<=n; j++)
    for(k=n; k>=j; k--)
      c[k] = (c[k]-c[k-1])/(x[k]-x[k-j]);

  return;
}


/* Newton's evaluation of the Lagrange polynomial

   newtonpoly(x, x[], c[], n):
   [0]  Let y = c[n]
   [1]  For k = n-1 down to 0 do [2]
   [2]     Replace y = c[k] + (x-x[k])*y
   [3]  Return y
*/
real
newtonpoly( real xx, real *x, real *c, int n )
{
  int k;  real y;
  y = c[n];
  for(k=n-1; k>=0; k--) y = c[k] + (xx-x[k])*y;
  return y;
}

/* Chebyshev's evaluation of the Lagrange polynomial

   chebyshevpoly(x, c[], n):
   [0]  Let y = c[0]
   [1]  If n>0 then do [2] through [8]
   [2]     Replace y = y + x*c[1]
   [3]     If n>1 then do [4] through [8]
   [4]        Let T[0] = 1
   [5]        Let T[1] = x
   [6]        For k = 2 to n do [2]
   [7]           Replace T[k] = 2*x*T[k-1] - T[k-2]
   [8]           Replace y = y + c[k]*T[k]
   [9]  Return y
*/
real
chebyshevpoly(real x, real *c, int n)
{
  real y,t,t0,t1; int k;
  y = c[0];
  if(n>0) {
    y += x*c[1];
    if(n>1) {
      t0 = 1;  t1 = x;
      for(k=2; k<=n; k++) {
	t = 2*x*t1 - t0;
	y += c[k]*t;
	t0 = t1;  t1 = t;
      }
    }
  }
  return y;
}

/* Piecewise constant evaluation

   pwconstant( x, x[], y[], n ):
   [0] For k = 0 to n-1 do [1]
   [1]    If x < (x[k]+x[k+1])/2 then return y[k]
   [2] Return y[n]
*/
real
pwconstant( real xx, real *x, real *y, int n )
{
  int k;
  for(k=0; k<n; k++) 
    if( xx < (x[k]+x[k+1])/2 ) return y[k];
  return y[n];
}

/*  Equispaced piecewise constant evaluation

   epwconstant( x, y[], a, b, n ):
   [0] Let k = floor( 0.5 + n*(x-a)/(b-a) )
   [1] Return y[k]
*/
real
epwconstant( real x, real *y, real a, real b, int n )
{
  int k;
  if(x<a) return y[0];
  if(x>b) return y[n];
  k = (int)(0.5 + n*(x-a)/(b-a));
  return y[k];
}

/* Piecewise linear evaluation
   pwlinear( x, x[], y[], N ):
   [0] If x < x[0] return y[0]
   [1] For k = 1 to N do [1]
   [2]   If x<x[k] then return (y[k]*(x-x[k-1])+y[k-1]*(x[k]-x))/(x[k]-x[k-1])
   [3] Return y[N]
*/
real
pwlinear( real xx, real *x, real *y, int n )
{
  int k;
  if( xx < x[0] ) return y[0];
  for(k=1; k<=n; k++)
    if( xx < x[k] ) 
      return (y[k]*(xx-x[k-1])+y[k-1]*(x[k]-xx))/(x[k]-x[k-1]);
  return y[n];
}