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

[robve]

Albert, in the Tanh-Sinh programs that I based my C and BASIC implementation on, I found several issues that can be corrected and optimized. Love to get your feedback on those. I also "played around" a bit with the Python mpmath Tanh-Sinh implementation, which is terrible (slower than Gauss-Legendre for several test cases when it should be faster and inaccurate results, perhaps more about that later...)

a. the Tanh-Sinh code I looked at attempts to prevent evaluating the function at points a and b. However, as I've explained in my previous reply, these checks do not reliably work. I can think of at least two possible options to fix this, though it might also be fine to leave them out entirely, because it is known and understood that Tanh-Sinh may fail for certain functions with singularities at the endpoints.

Option 1: changing the inner loop termination condition from this:

      if (r == 0 || r == 1)
        break;
      x = d*r;
      if (a+x > a)
        fp = f(a+x);
      if (b-x < b)
        fm = f(b-x);

to the following lines of code that check if |r|<MachEps instead of r==0 and by removing the unnecessary ad-hoc guards that wrongly used the old fp and fm values (i.e. reused previous points), which I did not find any justification for and I seriously doubt is okay since it may affect the result and the error estimation:

      if (1+r == 1 || r == 1)
        break;
      x = d*r;
      fp = f(a+x);
      fm = f(b-x);

This prevents x from getting too close to the endpoints a and b. This change takes the relative distance of the proximity of x to a and b into account, i.e. scaled with the size of the integration domain d=(b-a)/2, note that we have x=d*r.

Option 2: the first option can still fail if d<<1 and both a and b are large, i.e when |a-b|<<1 and a,b>>1 (or when c>>1 center). So:

      if (d+r == d || r == 1)
        break;
      x = d*r;
      if (c+x == c)
        break;
      fp = f(a+x);
      fm = f(b-x);

These are two possible options to reduce the chances of hitting a and b too close.

b. The code I looked at lacks some optimizations to further speed up the code, especially when exp() is expensive to compute, e.g. in software on BASIC pocket computers and calculators.

We can apply two "loop strength reductions", which do not negatively affect precision. The first strength reduction I already applied to the BASIC program. It is a change from:

double h;
int k = 0;
...
do {
    double p = 0, q, v;
    int j = 1;
    h = pow(2, -k); // h=1,1/2,1/4,1/8,...
    do {
      double t = exp(h*j);
      ...
      j += 1+(k>0);
    } while (fabs(q) > fabs(p*eps));

to the strength reduced form to remove pow(2, -k):

double h = 2;
int k = 0;
...
do {
    double p = 0, q, v;
    h /= 2;
    do {
      double t = exp(h*j);
      ...
      j += 1+(k>0);
    } while (fabs(q) > fabs(p*eps));
...
} while (e >= tol && k <= n);
...
return d*s*h; // we need the final h here

Then we apply a second loop strength reduction by observing that exp(h*j)=exp(h)^j then strength-reduce the power ^j away to obtain:

double h = 2;
int k = 0;
...
do {
    double p = 0, q, v;
    double exph, exphj;
    h /= 2;
    exph = exp(h);
    exphj = exph;
    if (k > 0)
      exph *= exph
    do {
      double t = exphj;
      ...
      exphj *= exph;
    } while (fabs(q) > fabs(p*eps));
...
} while (e >= tol && k <= n);
...
return d*s*h; // we need the final h here

Note that exp(h) in the outer loop has up to k==6 or k==7 values, namely \( \{e, e^{1/2}, e^{1/4}, e^{1/8}, e^{1/16}, e^{1/32}, e^{1/64} \} \) and can be precomputed to speed this up even more.

I experimented with these changes to verify them and everything looks good. More testing is necessary and additional improvements might be possible. More to come...

- Rob