/* Solution: Spline This */
/* E. Jason Riedy, ejr@cise.ufl.edu */
/*
  This type of system of equations is called tri-diagonal.  Gaussian
  elimination can be greatly simplified; only the next (or previous)
  row must be examined, and then only for a known number of columns.
  You can also optimize the storage and keep only the elements you
  need.  Below is the basic layout of what's kept from Ax=b.

    diag               b = ypp
        00 right       0
        011            1
    left 122           2
          233          3
           344         4
            45         5

  The routine below destroys b and replaces it with the solution.
*/

#include<stdlib.h>
#include<stdio.h>

enum { MAX_POINTS = 20 };

/* Find the correct interval */
int table_find( int n, double* x, double where )
{
  int high = n-1;
  int low = 0;
  int j;

  while (high - low > 1) {
    j = (high + low) / 2;
    if (x[j] >= where)
      high = j;
    else
      low = j;
  }

  j = low;

  /*
  printf( "Table find: x[%d] (%.2f) <= where (%.2f) <= x[%d] (%.2f)\n",
	  j, x[j], where, j+1, x[j+1]);
  /**/

  return j;
}

double eval_cspline( int n, double* x, double* y, double* ypp, double where )
{
  double A, B, C, D, xint;
  double out;
  int j = table_find( n, x, where );

  xint = x[j+1] - x[j];
  A = (x[j+1] - where)/xint;
  B = 1-A;
  C = (A*A*A-A) * xint*xint / 6.0;
  D = (B*B*B-B) * xint*xint / 6.0;

  /*
  printf( "Coefficients for %.2f: A = %.4f, B = %.4f, C = %.4f, D = %.4f\n",
	  where, A, B, C, D );
  printf( "ypp[%d] = %.4f, ypp[%d] = %.4f\n", j, ypp[j], j+1, ypp[j+1] );
  /**/

  out = A * y[j] + B * y[j+1] + C * ypp[j] + D * ypp[j+1];
  return out;
}

/*----------------------------------------------------------------------*/

void fill_b( int n, double* x, double* y, double* b )
{
  int i;
  double oldfrac, newfrac;

  oldfrac = (y[1] - y[0]) / (x[1] - x[0]);
  for (i = 1; i < n-1; i++) {
    newfrac = (y[i+1] - y[i]) / (x[i+1] - x[i]);
    b[i] = newfrac - oldfrac;
    oldfrac = newfrac;
  }
  b[0] = 0.0;
  b[n-1] = 0.0;
}

void fill_dlr( int n, double* x, double* diag, double* left, double* right )
{
  int i;

  for (i = 0; i < n-1; i++) {
    diag[i+1] = (x[i+2] - x[i]) / 3.0;
    left[i] = (x[i+1] - x[i]) / 6.0;
    right[i] = (x[i+1] - x[i]) / 6.0;
  }
}

void solve_cspline( int n, double* x, double* y, double* ypp )
{
  int i;
  double diag[MAX_POINTS];
  double left[MAX_POINTS];
  double right[MAX_POINTS];
  double m;

  fill_dlr( n, x, diag, left, right );
  fill_b( n, x, y, ypp );

  /* Eliminate the elements in left going down... */
  for (i = 1; i < n-1; i++) {
    m = diag[i];
    right[i] /= m;
    ypp[i] /= m;
    diag[i+1] -= right[i] * left[i];
    ypp[i+1] -= ypp[i] * left[i];
  }

  /* and then eliminate those in right going up.  */
  ypp[n-1] = 0.0;
  for ( i = n-2; i >= 0; i-- )
    ypp[i] -= ypp[i+1] * right[i];

  /*
    Note that the values being eliminated are never actually changed.
    They are never read again, so why bother setting them to zero?
  */
}

void read_points( int* n, double* x, double* y ) 
{
  int i;
  
  printf("Number of points: ");
  scanf("%d", n);

  for (i = 0; i < *n; i++) {
    printf("Point %d: ", i+1);
    scanf("%lg %lg", &x[i], &y[i]);
  }
}

void read_eval_print( int n, double* x, double* y, double* ypp )
{
  double where;
  double out;

  putchar('\n');
  for (;;) {
    printf("Evaluate at: ");
    scanf("%lg", &where);

    if (feof(stdin))
      return;

    out = eval_cspline( n, x, y, ypp, where );
    printf("The curve passes through (%.2f, %.2f).\n", where, out);
  }
}

void main(void)
{
  int n;
  double x[MAX_POINTS], y[MAX_POINTS], ypp[MAX_POINTS];

  read_points( &n, x, y );
  solve_cspline( n, x, y, ypp );
  read_eval_print( n, x, y, ypp );
  putchar('\n');
}