#include <stdio.h>
#include <math.h>
/* spline.c : 1999 3 9 by Oguni */
#define N 7
#define M 26
main(){
  double x[N],y[N],xnew[M],ynew[M],yy[M],df[2];
  double tol=0.000001;
  int i, n=7, m=26;
  df[0]=1.; df[1]=1.;
  for (i=0; i<n; i++){
    x[i]=(double)i-1.0; y[i]=sin(x[i]); }
  for (i=0; i<m; i++){
    xnew[i]=((double)i-1.0)*0.25;
    yy[i]=sin(xnew[i]); }
  spline(n,m,x,y,df,xnew,ynew,tol);
  printf("original and estimated\n");
  for (i=0; i<m-4; i++){
    printf("%f %f\n", ynew[i], yy[i]); }
}
/* function */
spline(n,m,x,y,df,xnew,ynew,tol)
double x[N],y[N],xnew[M],ynew[M],df[2],tol;
int n,m;
{
  double h[N],lamda[N],u[N],d[N],c[N][3],em[N+1],A[N][N],z[N];
  int i,j,k;
  double xk;
  for (i=0; i<n-1; i++){
    h[i]=x[i+1]-x[i]; }
  for (i=1; i<n-1; i++){
    lamda[i]=h[i]/(h[i-1]+h[i]); u[i]=1.0-lamda[i]; }
  lamda[0]=1.; u[0]=1.;
  for (i=1; i<n-1; i++){
    d[i]=6./(h[i-1]+h[i])*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]);}
  d[0]=6./h[0]*(y[1]-y[0]-df[0])/h[0];
  d[n-1]=-6./h[n-2]*(y[n-1]-y[n-2]-df[1])/h[n-2];
  for (i=1; i<n; i++){
    A[i][i]=2.; A[i][i-1]=u[i]; A[i-1][i]=lamda[i-1]; }
  A[0][0]=2.;
  gauss0(A,d,z,tol,n);
  for (i=0; i<n; i++){
    em[i]=z[i]; }
  em[n]=0.;
  c[0][2]=(em[1]-em[0])/6. /h[0];
  for (i=0; i<n-1; i++){
    c[i][1]=em[i]/2.; c[i][2]=(em[i+1]-em[i])/6. /h[i];
    c[i][0]=(y[i+1]-y[i])/h[i]-c[i][1]*h[i]-c[i][2]*h[i]*h[i]; }
  for (j=0; j<m; j++){
    k=1;
    while (k<n-1){
      if (xnew[j]>x[k] && xnew[j]<=x[k+1]){
        xk=xnew[j]-x[k];
        ynew[j]=y[k]+c[k][0]*xk+c[k][1]*xk*xk+c[k][2]*xk*xk*xk;
        break;}
       k++; }
  }
  return;
}
/* function */
gauss0(A, b, x, tol, n)
double A[N][N], b[N], x[N], tol;
int n;
{
  int i, j, k;
  double akk, aik, s;
  for (k=0; k<n-1; k++){
    if (fabs(A[k][k])<tol){
      printf("The pivot is too small at step # %d", k); return; }
    akk=1.0/A[k][k];
    for (i=k+1; i<n; i++){
      aik=-A[i][k]*akk;
      for (j=k+1; j<n; j++){
        A[i][j]+=aik*A[k][j]; }
      b[i]+=aik*b[k];
    }
    for (j=k+1; j<n; j++){
      A[k][j]*=akk; }
    b[k]*=akk;
  }
  x[n-1]=b[n-1]/A[n-1][n-1];
  for (k=n-2; k>=0; k--){
    s=0.0;
    for (j=k+1; j<n; j++){ s+=A[k][j]*x[j]; }
    x[k]=b[k]-s;
  }
}