(PC-12xx~14xx) qthsh Tanh-Sinh quadrature

[robve]

Here is my latest contribution with integrated Tanh-Sinh, Exp-Sinh, Sinh-Sinh rules. The Tanh-Sinh rule has the same Michalski & Mosig abscissas and weights, but reformulated to combine the three quadratures into one routine (see the updated doc for the formulas):

Code:

// Tanh-Sinh bounded [a,b], Exp-Sinh (one inf either side), Sinh-Sinh [-inf,inf]
double quad(double (*f)(double), double a, double b, double n, double eps) {
  const double tol = 10*eps;
  double c = 0, d = 1, s, sign = 1, e, v, h = 2;
  int k = 0, mode = 0; // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
  if (b < a) { // swap bounds
    v = b;
    b = a;
    a = v;
    sign = -1;
  }
  if (isfinite(a) && isfinite(b)) {
    c = (a+b)/2;
    d = (b-a)/2;
    s = f(c);
  }
  else if (isfinite(a)) {
    mode = 1; // Exp-Sinh
    c = a;
    s = f(a + 1);
  }
  else if (isfinite(b)) {
    mode = 1; // Exp-Sinh
    c = b;
    s = f(b - 1);
    d = -d;
    sign = -sign;
  }
  else {
    mode = 2; // Sinh-Sinh
    s = f(0);
  }
  do {
    double p = 0, q, t, eh;
    h /= 2;
    eh = exp(h);
    t = eh/2;
    if (k > 0)
      eh *= eh;
    do {
      double r, w, x, y;
      q = 0;
      r = w = exp(t-.25/t); // = exp(sinh(j*h))
      if (mode != 1) {      // if Tanh-Sinh or Sinh-Sinh
        w += 1/w;           // = 2*cosh(sinh(j*h))
        r -= 1/r;           // = 2*sinh(sinh(j*h))
        if (mode == 0) {    // if Tanh-Sinh
          r /= w;           // = tanh(sinh(j*h))
          w = 4/(w*w);      // = 1/cosh(sinh(j*h))^2
        }
        else {              // if Sinh-Sinh
          r /= 2;           // = sinh(sinh(j*h))
          w /= 2;           // = cosh(sinh(j*h))
        }
        x = c - d*r;
        if (x > a) {
          y = f(x);
          if (isfinite(y))  // if f(x) is finite, add to local sum
            q += y*w;
        }
      }
      else {                // Exp-Sinh
        x = c + d/r;
        if (x > a) {
          y = f(x);
          if (isfinite(y))  // if f(x) is finite, add to local sum
            q += y/w;
        }
      }
      x = c + d*r;
      if (x < b) {
        y = f(x);
        if (isfinite(y))    // if f(x) is finite, add to local sum
          q += y*w;
      }
      q *= t+.25/t;         // q *= cosh(j*h)
      p += q;
      t *= eh;
    } while (fabs(q) > eps*fabs(p));
    v = 2*s;
    s += p;
    ++k;
  } while (fabs(v-s) > tol*fabs(s) && k <= n);
  e = fabs(v-s)/(fabs(s)+eps);
  return sign*d*s*h; // result with estimated relative error e
}

This was highly optimized while keeping the code compact. Optimizations include strength reductions, branch eliminations, and variable reuse when possible to reduce CPU register pressure.

The partially optimized code for the inner loop is perhaps a bit easier to comprehend:

Code:

    t = eh = exp(h);
    if (k > 0)
      eh *= eh;
    do {
      double ch, u, r, x, w;
      ch = (t+1/t)/2; // = cosh(j*h)
      w = r = u = exp((t-1/t)/2); // = exp(sinh(j*h))
      if (mode != 1) {
        u += 1/u; // = 2*cosh(sinh(j*h))
        r -= 1/r;   // = 2*sinh(sinh(j*h))
        if (mode == 0) {
          r /= u; // = tanh(sinh(j*h))
          w = 4/(u*u); // = 1/cosh(sinh(j*h))
        }
        else {
          r /= 2;   // = sinh(sinh(j*h))
          w = u/2; // = cosh(sinh(j*h))
        }
      }
      x = d*r;
      q = 0;
      if (c+x < b) {
        double y = f(c+x);
        if (isfinite(y))
          q = y*w;
      }
      if (mode == 1)
        x = -d/r;
      if (c-x > a) {
        double y = f(c-x);
        if (isfinite(y))
          q += mode == 1 ? y/w : y*w;
      }
      q *= ch;
      p += q;
      t *= eh;
    } while (fabs(q) > eps*fabs(p));

I have verified the optimized routine with several examples, but more testing will be necessary to determine its numerical properties.

- Rob

Edit: minor change to fix mirrored bound checks.