#include <stdio.h>
/* simpson.c : 1999 3 3 by Y.Yamazaki and T.Oguni */
double simpson(double (*f)(), double a, double b, int n)
{
  int i;
  double x, h, h2, sum;
  h = (b-a)/n;  h2 = h+h;
  sum = (*f)(a)+(*f)(b);
  for (i=1, x=a+h; i<n; i+=2, x+=h2)
    sum += 4*(*f)(x);
  for (i=2, x=a+h2; i<n; i+=2, x+=h2)
    sum += 2*(*f)(x);
  return (h/3)*sum;
}
double func(double x)
{ return x*x*x*x*x*x*x; }
void main()
{ printf("f(x)dx = %f\n", simpson(func, 0., 1., 4)); }