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

[Albert Chan]

Some issues I noticed for Tanh-Sinh quadrature.

Let x = tanh(pi/2*sinh(t))

∫(f(x), x=0 .. 1) = pi/2 * ∫(f(tanh(pi/2*sinh(t))) * (cosh(t)/cosh(pi*sinh(t)/2)^2), t = 0 .. inf)

Numerically, t does not need to go too high (certainly not inf)

lua> function ts(t) return tanh(pi/2*sinh(t)) end
lua> ts(2.8), ts(3), ts(3.2)
0.9999999999866902       0.999999999999957       1

At t=3, with double precision, we are very close to evaluating the end-point, f(1)
If f is already u-transformed, say, from f0, t is even smaller (to avoid touching f0(1))

∫(f0(x), x=0..1) = ∫(6*u*(1-u) * f0(u*u*(3-2*u)), u=0..1)       // RHS integrand = f(u)

lua> function u(u) return u*u*(3-2*u) end
lua> u(ts(2.2)), u(ts(2.4)), u(ts(2.6))
0.9999999999917426       0.9999999999999855       1

So, what is the big deal ?
If f(1) = ∞, u-transformation might make the algorithm fail.

>>> from mpmath import *
>>> f = lambda x: 1/sqrt(-ln(x))
>>> quadts(f, [0,1], error=True)
(mpf('1.7724538503750329'), mpf('1.0e-10'))
>>>
>>> uf = lambda u: 6*u*(1-u) * f(u*u*(3-2*u))
>>> quadts(uf, [0,1], error=True)
Traceback (most recent call last):
...
ZeroDivisionError