/*    ATK IV - Numerical Programming (2006)      */
/*      Excercise 7.1 -  Trapezoidal Rule        */

#include<stdio.h>
#include<math.h>
#define N 10000

double f(double x){return sin(x)*sin(x);}

double fi(double x){return 0.5*x-0.25*sin(2.0*x);}

double maxi(double x, double y){
   if(x>y) return x;
   else return y;
}

void trapz(double y[],double h,int n,double yi[]);
void tabulate(double (*g)(double),double x1,double x2,int n,double y[]);

int main(void){
  double x1 = 0.0;
  double x2 = 2.0*acos(-1.0);
  double h = (x2-x1)/(N-1);
  double y[N],yi[N],ye[N];
  double err,maxerr;
  int i;

  tabulate(f,x1,x2,N,y);
  tabulate(fi,x1,x2,N,ye);
  trapz(y,h,N,yi);

  for(i=0;i<=N-1;i++){
    err = fabs(yi[i]-ye[i]);
    maxerr = maxi(maxerr,err);
    printf("%f %f %f %f\n",y[i],yi[i],ye[i],err);
  }
  printf("\nMax. error = %LG\n",maxerr);
  return 0;
}

void trapz(double y[],double h,int n,double yi[]){
  int i;

  yi[0]=h*0.5*y[0];

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

  yi[n-1]=yi[n-2]+h*0.5*y[n-1];
}

void tabulate(double (*g)(double x),double x1,double x2,int n,double y[]){
  int i;

  double h=(x2-x1)/(n-1);
  for(i=0;i<=n-1;i++){
    y[i]=g(x1+h*(i));
  }
}