/* de_half_infinite.c */
/* 半無限積分に対するDE公式 */
#include <stdio.h>
#include <math.h>
#define NMAX 1024
#define EPS 1.0e-15
double func(double x);
double func_de_half(double x, double a);

int main(void)
{
  int j,n,m;
  double a,s,integral,integral0,h,x,fx,fx0, f0,w;
  //  double I = M_PI_4; /* 真値 */

  printf("%10s %5s %5s %16s %15s\n","h", "m", "n", "integral","difference");
  a = 1.0;
  integral0 = 0.0;
  for(h = 1.0;h >= 0.0625; h *= 0.5){
    x = 0.0;
    n = 0;
    f0 = func_de_half(x,a);
    s = f0;
    fx0 = f0;
    x += h;
    n++;
    fx = func_de_half(x,a);
    s += fx;
    fx0 = fx;
    w = fabs(f0) + fabs(fx);
    
    while(w > EPS && n < NMAX){
      x += h;
      n++;
      fx0 = fx;
      fx = func_de_half(x,a);
      s += fx;
      w = fabs(fx0) + fabs(fx);
    }

    x = -h;
    m = 1;
    fx = func_de_half(x,a);
    s += fx;
    fx0 = fx;  
    x -= h;
    m++;
    fx0 = fx;
    fx = func_de_half(x,a);
    s += fx;
    w = fabs(fx0) + fabs(fx);
    
    while(w > EPS && m < NMAX){
      x -= h;
      m++;
      fx0 = fx;
      fx = func_de_half(x,a);
      s += fx;
      w = fabs(fx0) + fabs(fx);
    }

    integral = h * s;
    printf("%10.6f %5d %5d %20.15f %10.3e\n",h,m,n,integral,integral-integral0);    
    //    printf("h = %6.4f m = %3d n = %3d int = %20.15f error = %10.3e\n",h,m,n,h*s,I-h*s);
    integral0 = h * s;

  }
  
  return 0;
}

double func(double x)
{
  return 1.0 / (1.0 + x * x);

}

double func_de_half(double x,double a)
{
  return M_PI_2 * func(a + exp(M_PI_2 * sinh(x))) * cosh(x) * exp(M_PI_2 * sinh(x)); 
}