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

int main(void)
{
  int j,n,m;
  
  double s,h,x,fx,fx0, f0,w;
  double integral,integral0;
  //  double I = M_PI; /* Iに真値を代入。M_PIはmath.hで定義されている円周率 */

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

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

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

double func(double x)
{
  //  return 1.0 /(1.0 + x * x * x * x);
  return 1.0 /(1.0 + x * x);
  //  return (x*x + x + 1.0) * exp(-x * x);
  //  return cos(x) * exp(-x * x);
  //  return cos(10 * x) * exp(-x * x);
}

double func_de(double x)
{
  return 0.5 * M_PI * func(sinh(0.5 * M_PI * sinh(x))) * cosh(x)
    * cosh(0.5 * M_PI * sinh(x)); 
}