/*
 * example7.c --- DE 公式
 *
 *        ∞
 *  I = ∫   1/(1+x^2) dx = π
 *       -∞
 *
 *  コンパイル: cc -o example7 example7.c
 *
 */

#include <stdio.h>
#include <math.h>

typedef double ddfunction(double);

double pi, halfpi;

/* φ */
double phi2(double t)
{
  return sinh(halfpi * sinh(t));
}

/* 2乗 */
double sqr(double x) { return x * x; }

/* φ' */
double dphi2(double t)
{
  return halfpi * cosh(t) * cosh(halfpi * sinh(t));
}

/* DE 公式による (-∞,∞) における定積分の計算 */
double de(ddfunction f, double h, double N)
{
  int n;
  double t, S, dS;
  S = f(phi2(0.0)) * dphi2(0.0);
  for (n = 1; n <= N; n++) {
    t = n * h;
    dS = f(phi2(t)) * dphi2(t) + f(phi2(-t)) * dphi2(- t);
    S += dS;
  }
  return h * S;
}

/* テスト用の被積分関数 その1 */
double f(double x)
{
  return 1.0 / (1.0 + x * x);
}

int main(void)
{
  int m, N;
  double h, IhN, I;

  pi = 4.0 * atan(1.0); halfpi = pi / 2;

  printf("1 / (1 + x^2) の積分)\n");
  I = pi;
  /* |t|≦4 まで計算することにする */
  h = 1.0; N = 4;
  /* h を半分, N を倍にして double exponential formula で計算してゆく */
  for (m = 1; m <= 10; m++) {
    IhN = de(f, h, N);
    printf("h=%f, I_h=%25.16f, I_h-I=%e\n", h, IhN, IhN - I);
    h /= 2; N *= 2;
  }

  return 0;
}