/*    ATK IV - Numerical Programming (2006)      */
/*     Excercise 7.3 -  Spline Integration       */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#define N 100

float f(float x){return exp(-x);}

void tabulate(float (*f)(float x),float x1,float x2,int n,float x[],float y[]);
void sint(float x[],float y[],float ydd[],float yi[],int n);

int main(void){
  float x1=0.0;
  float x2=3.0;
  float h=(x2-x1)/(N-1);
  float x[N+1],y[N+1],ydd[N+1],yi[N+1];
  int i;

  tabulate(f,x1,x2,N,x,y);
  spline(x,y,N-1,-f(x1),-f(x2),ydd);
  sint(x,y,ydd,yi,N);

  for(i=1;i<=N-1;i++){
    printf("%f %f %f\n",x[i],yi[i],-f(x[i])+f(x1));
  }

  return 0;
}

void tabulate(float (*f)(float x),float x1,float x2,int n,float x[],float y[]){
  int i;
  float h=(x2-x1)/(n-1);

  for(i=1;i<=n;i++){
    x[i]=x1+h*(i-1);
    y[i]=f(x[i]);
  }
}

void sint(float x[],float y[],float ydd[],float yi[],int n){
  int i;
  float dx;

  yi[1]=0.0;
  for(i=2;i<=n-1;i++){
    dx=x[i+1]-x[i];
    yi[i]=yi[i-1]+0.5*dx*(y[i]+y[i+1]-1.0/12.0*dx*dx*(ydd[i]+ydd[i+1]));
  }
}