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

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

  //  I = 4.442882938158366101788487867452204227448; /* B(1/4,3/4) */
  /* Iに真値を代入。真値がわからないときは1行上を削除 */

  printf("%10s %5s %5s %20s %15s\n","h","m","n","integral","difference"); 
  s = interval(xa,xb,h);
  //  printf("[%5.3e,%5.3e]\n",*xa,*xb);

  
  count = 0;
  while(fabs(s-s0)  > EPS2){
    s0 = s;
    h *= 0.5;
    interval(xa,xb,h);
    s = trapezoidal(h,xa,xb,s);
    //    printf("%20.15f %15.8e\n",s,s-I); /* 真値がわかっているとき */
    printf("%20.15f %15.8e\n",s,s-s0); 
  }

  printf("関数値計算回数=%d\n",count);
  
  free(xa);
  free(xb);

  return 0;
}

double interval(double *xa, double *xb, double h)
{
  int m, n;
  double x, fx, fx0, fxp, w, 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);
  *xb = h * n;

  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);
  *xa = -h * m;

  if (h >= 1.0){
    printf("%10.8f %5d %5d %19.15f\n",h,m,n,s);
  } else {
    printf("%10.8f %5d %5d",h,m,n);
  }	
  
  return s;
}

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

  x = *xa;
  w = 0.0;

  while(x  <=  *xb){
    fx = func(x);
    count++;
    w += fx;
    x += h;
  }

  s = h * w;

  return s;
}

double func(double x)
{
  return M_PI_2*cosh(x)/(sqrt(exp(M_PI_2*sinh(x)))*cosh(M_PI_2*sinh(x)));
  /* 桁落ち防止変換 p.165 (7.32) */
}