/*
 * example6.c --- DE 公式
 *
 *         1      2 (1/2)
 *  I1 = ∫   (1-x )     dx = π/2
 *        -1
 *
 *         1      2 (-1/2)
 *  I2 = ∫   (1-x )      dx = π
 *        -1
 *
 *  いずれも端点に特異性があって、古典的な数値積分公式はうまく行かない。
 *  double exponential formula (DE公式) ならばうまく計算できる。
 *
 *  コンパイル: cc -o example6 example6.c
 *
 */

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

typedef double ddfunction(double);

double pi, halfpi;

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

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

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

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

/* テスト用の被積分関数 その1 */
double f1(double x)
{
  return sqrt(1 - x * x);
}

/* テスト用の被積分関数 その2 */
double f2(double x)
{
  if (x >= 1.0 || x <= -1.0)
    return 0;
  else
      return 1 / sqrt(1 - x * x);
}

void test(ddfunction f, double I)
{
  int m, N;
  double h, IhN;

  /* |t|≦3 まで計算することにする */
  h = 1.0; N = 3;
  /* 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;
  }
}

int main(void)
{
  pi = 4.0 * atan(1.0); halfpi = pi / 2;

  printf("test1 (sqrt(1-x^2) の積分)\n");
  test(f1, halfpi);

  printf("test2 (1/sqrt(1-x^2) の積分)\n");
  test(f2, pi);

  return 0;
}