/*
 * nc.c --- Newton-Cotes の積分公式: 複合中点公式, 複合台形公式, 複合Simpson公式
 */

typedef double ddfunction(double);

double midpoint(ddfunction, double, double, int);
double trapezoidal(ddfunction, double, double, int);
double simpson(ddfunction, double, double, int);

/* 関数 f の [a,b] における積分の複合中点公式による数値積分 M_N */
double midpoint(ddfunction f, double a, double b, int N)
{
  int j;
  double h, M;
  h = (b - a) / N;
  M = 0.0;
  for (j = 1; j <= N; j++) M += f(a + (j - 0.5) * h);
  M *= h;
  return M;
}

/* 関数 f の [a,b] における積分の複合台形公式による数値積分 T_N */
double trapezoidal(ddfunction f, double a, double b, int N)
{
  int j;
  double h, T;
  h = (b - a)/ N;
  T = (f(a) + f(b)) / 2;
  for (j = 1; j < N; j++) T += f(a + j * h);
  T *= h;
  return T;
}

/* 関数 f の [a,b] における積分の複合 Simpson 公式による数値積分 S_{N} */
double simpson(ddfunction f, double a, double b, int N)
{
  int m = N / 2;
  return (trapezoidal(f, a, b, m) + 2 * midpoint(f, a, b, m)) / 3;
}