/* Integral of exp(-PI*x*x)^2 sinc(PI*x)^2 on the complement of [-1,1], done
   with Romberg integration */

#include <assert.h>
#include <math.h>
#include <stdio.h>

#define PI 3.14159265358979323
#define JMAX 15			/* Recursion depth limit. */

/* Function to integrate: */
double f(double t) {
  double x, y;
  if(t==0)
    y = 1.0;
  else {
    x = PI*t;
    y = exp(-t*t)*sin(x)/x;
  }
  return y*y;
}

/* Recursive composite trapezoidal rule increment on 2^J segments of [a,b]: */
double rcti(double a, double b, int j)
{
  double h, sum;  int n, k;
  assert(a<b);  assert(j>=0);  assert(j<30);
  n = 1<<j;  h = (b-a)/n;
  if(j==0) sum = (f(a)+f(b))/2;
  else { sum = 0; for(k=1; k<n; k+=2) sum += f(a+k*h);  }
  return h*sum;
}

main(void)
{
  double r[JMAX][JMAX], phi;
  int i, j, k;
  double a=-5.3, b=5.3;	/* Integrate f() on [a,b] */
  double err, tol=5.0e-14;

  r[0][0] = rcti(a,b,0);	/* Initial trapezoid rule value */
  for(j=1; j<JMAX; j++) {
    r[j][0]=0.5*r[j-1][0]+rcti(a,b,j); /* recursive trapezoid rule */
    for(k=1; k<=j; k++)	/* Romberg integration */
      r[j][k] = r[j][k-1] + (r[j][k-1] - r[j-1][k-1])/((1<<(2*k))-1.0);
    err = fabs( r[j][j] - r[j-1][j-1] );
    phi = r[j][j];		/* Value of the integral */
    printf("phi[%3.1f,%3.1f]=%20.14f    err=%8.2e, J=%d\n", 
	   a, b, phi, err, j );
    if(err<tol) break;	/* escape j-loop */
  }
  return 0;
}