/* de_fourier.c */
/* フーリエ型積分に対するDE公式 */
#include <stdio.h>
#include <math.h>
#define IMAX 5
#define NMAX 128
#define EPS 1.0e-12
double f(double x);
double g(int k,double h);
double M;
/* 被積分関数 f(x) * sin(x) */


int main(void)
{
  int i,n,m;
  double s,s0,h,x,t,gx,gx0,w;
  //  double I = M_PI_2; /* 真値 */

  printf("%10s %10s %5s %5s %16s %15s\n","M","h", "n", "m", "integral","difference");
  h = 1.0;
  for(i = 1; i <= IMAX;i++){
    s0 = s;
    h *= 0.5;
    M = M_PI / h;
    n = 0;
    gx = g(n,h);
    s = gx;
    gx0 = gx;
    n++;
    gx = g(n,h);
    s += gx;
    w = fabs(gx0) + fabs(gx);

    
    while(w > EPS && n < NMAX){
      n++;
      gx0 = gx;
      gx = g(n,h);
      s += gx;
      w = fabs(gx0) + fabs(gx);
    }

    m = 1;
    gx = g(-m,h);
    s += gx;
    gx0 = gx;  
    m++;
    gx0 = gx;
    gx = g(-m,h);
    s += gx;
    w = fabs(gx0) + fabs(gx);
    
    while(w > EPS && m < NMAX){
      t -= h;
      m++;
      gx0 = gx;
      gx = g(-m,h);
      s += gx;
      w = fabs(gx0) + fabs(gx);
    }

    printf("%10.6f %10.6f %5d %5d %20.15f %10.3e\n",M,h,n,m,s,s-s0);
  }
  
  return 0;
}

double f(double x)
{
  double y;

  if (x == 0)
    y = 1.0;
  else
    y = 1.0 / x;
  
  return y;
}


double g(int k,double h)
{
  double t,y;

  if (k == 0)
    y = M_PI * f(M/6) * sin(M/6) /2;
  else{
    t = k * h;
    y = f(M * t / (1.0 - exp(-6 * sinh(t))));
    y *= sin(k * M_PI /(1.0 - exp(-6 * sinh(t))));
    y *= (1-(1+6*t*cosh(t))*exp(-6 * sinh(t)));
    y /= ((1.0 - exp(-6 * sinh(t))) * (1.0 - exp(-6 * sinh(t))));
    y *= M_PI;
  }
  
  return y;
}