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

[Albert Chan]

Revisiting original quad(), where we interpolate from the center, we may be able to improve it.
double_exponential.py is doing the same way, so the patch might be useful.

This was the snippet from before, started with r = w = exp(t-.25/t); // = exp(sinh(j*h))
For tanh-sinh, r ← (r-1/r) / (r+1/r), which may suffer rounding errors > 1/2 ulp

(04-11-2021 03:05 AM)robve Wrote:  

Code:

        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))
        }

Instead, we calculate from the edge: r ← (r-1/r) / (r+1/r) = 1 - 2/(1+r*r)
This is the proposed patch, where updated r should be much more accurate, at the edges.

Code:

    if (mode == 0) {        // if Tanh-Sinh
        r = 2/(1+r*r);
        w = w * r;
        r = 1 - r;          // = tanh(sinh(j*h))
        w = w * w;          // = 1/cosh(sinh(j*h))^2
    }
    else {                  // if Sinh-Sinh
        r = (r-1/r) / 2;    // = sinh(sinh(j*h))
        w = (w+1/w) / 2;    // = cosh(sinh(j*h))
    }

(04-12-2021 01:50 PM)Albert Chan Wrote:  lua> Q = require 'quad'
lua> f = function(x) return 1/sqrt(1-x*x) end

lua> Q.qthsh(Q.count(f), -1, 1), Q.n
3.141592644894014       59
lua> Q.quad(Q.count(f), -1, 1), Q.n
3.1415926047697873     115

This is still not as good as interpolate from the edge, but it improved a lot.
Redoing the above example, with the proposed patch.

lua> O = require 'quad0' -- original quad(), interpolate from center, but more accurate r
lua> O.quad(O.count(f), -1, 1), O.n
3.1415926448906006       59