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

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

typedef double ddfunction(double);

double pi, halfpi;

//  tanh(t)-1 の計算 (t≫1の場合用)
double tanhm1(double t)
{
  if (t <= 354)
    return - 2.0 / (1.0 + exp(2.0 * t));
  else {
    printf("tahnm1(): %g は大きすぎる\n", t);
    return 0;
  }
}

//  tanh(t)+1 の計算 (t≪-1の場合用)
double tanhp1(double t)
{
  if (t >= - 354)
    return 2.0 / (1.0 + exp(- 2.0 * t));
  else {
    printf("tahnp1(): %g は小さ過ぎる\n", t);
    return 0;
  }
}

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

/* ξ(t)=φ(t)-1 */
double xi(double t)
{
  return tanhm1(halfpi * sinh(t));
}

/* η(t)=φ(t)+1 */
double eta(double t)
{
  return tanhp1(halfpi * sinh(t));
}

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

/* φ' */
double dphi1(double t)
{
  // printf("dphi1: %g\n", halfpi * cosh(t) / sqr(cosh(halfpi * sinh(t))));
  return halfpi * cosh(t) / sqr(cosh(halfpi * sinh(t)));
}

/* DE 公式による (a,b) における定積分の計算, x=a,b で桁落ち対策 */
double de3(ddfunction f, ddfunction g, ddfunction H,
	   double a, double b, double h, double N)
{
  int n;
  double p, q;
  double t, S, dS;
  p = (b - a) / 2.0; q = (a + b) / 2.0;
  S = 0.0; 
  for (n = -N; n <= N; n++) {
    t = n * h;
    if (t > 0.5)
      dS = dphi1(t) * g(p * xi(t));
    else if (t < - 0.5)
      dS = dphi1(t) * H(p * eta(t));
    else
      dS = dphi1(t) * f(p * phi1(t) + q);
    S += dS;
  }
  return p * h * S;
}

/* テスト用の被積分関数 その2 */
double f(double x)
{
  return 1.0 / sqrt(1.0 - sqr(x));
}

double g(double y)
{
  return 1.0 / sqrt(- y * (y + 2.0));
}

double h(double y)
{
  return 1.0 / sqrt(y * (2.0 - y));
}

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

  /* |t|≦4 まで計算することにする */
  h = 1.0; N = 4;
  /* h を半分, N を倍にして double exponential formula で計算してゆく */
  for (m = 1; m <= 10; m++) {
    IhN = de3(f, g, H, - 1.0, 1.0, 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("test2 (1/sqrt(1-x^2) の積分)\n");
  test(f, g, h, pi);

  return 0;
}