/* Double exponential formula */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX 100
#define EPS1 1.0e-16
#define EPS2 1.0e-12
#define PI 3.1415926535897932
double interval(double *xa, double *xb);
double trapezoidal(double h, double *xa, double *xb, double s);
double func(double x);

int main(void)
{
  double *xa, *xb, s = 0.0, h = 1.0, s0 = 0.0, I; 

  xa = malloc(sizeof(double));
  xb = malloc(sizeof(double));

  I = 4.442882938158366101788487867452204227448; /* B(1/4,3/4) */

  s = interval(xa,xb);
  printf("[%5.3e,%5.3e]\n",*xa,*xb);
  printf("%15.8e %20.15f %15.8e %15.8e\n",h,s,s-I,s-s0);

  while(fabs(s-s0) > EPS2){
    s0 = s;
    h *= 0.5;
    s = trapezoidal(h,xa,xb,s);
    printf("%15.8e %20.15f %15.8e %15.8e\n",h,s,s-I,s-s0);
  }

  free(xa);
  free(xb);

  return 0;
}

double interval(double *xa, double *xb)
{
  int m, n;
  double x, fx, fx0, fxp, w, h = 1.0, s;

  x = 0.0;
  fx = func(x);
  fx0 = fx;
  fxp = fx;
  x += h;
  fx = func(x);
  n = 1;
  s = h * (fxp + fx);
  w = fabs(fxp) + fabs(fx);

  /* determination of upper boundary  */
  while(w>EPS1 && n<MAX){
    fxp = fx;
    x += h;
    fx = func(x);
    //    printf("%10.3f %10.3e %20.15g\n",x,fx,1.0-tanh((PI/2.0)*sinh(x)));
    n++;
    s += h * fx;
    w = fabs(fxp) + fabs(fx);
  }

  *xb = h * (n - 1);

  fxp = fx0;
  x = -h;
  fx = func(x);
  m = 1;
  s += h * fx;
  w = fabs(fxp) + fabs(fx);

  /* determination of lower boundary  */
  while(w>EPS1 && m<MAX){
    x -= h;
    fxp = fx;
    fx = func(x);
    //    printf("%10.3f %10.3e %20.15g\n",x,fx,tanh((PI/2.0)*sinh(x))+1.0);
    m++;
    s += h * fx;
    w = fabs(fxp) + fabs(fx);
  }

  *xa = -h * (m - 1);

  return s;
}

double trapezoidal(double h, double *xa, double *xb, double s)
{
  double x, w, fx;

  x = *xa + h;
  w = 0.0;

  /* midpoint rule */
  while(x < *xb){
    fx = func(x);
    w += fx;
    x += 2 * h;
  }

  s = 0.5 * s + h * w;
  return s;
}

double func(double x)
{
  return (PI/2.0)*cosh(x)/(sqrt(exp((PI/2.0)*sinh(x)))*cosh((PI/2.0)*sinh(x)));
}