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

[robve]

I also noticed other issues in the VB code in the Article, which I have not mentioned in my writeup, such as computing the relative error with:
errval = Abs(2 * sslast / ss - 1)
However, the distributive law does not apply to floating point arithmetic. This should be:
errval = Abs((2 * sslast - ss) / ss)
When 2*sslast ~ ss this is important.

However, this is a bit pointless to say, because both pieces of code may cause division by zero when the updated quadrature sum ss is zero. A function such as f(x)=x-.5 over [0,1] fails to integrate. It can get worse, because a division by zero produces errval=NaN. NaN is then compared to a threshold value. Comparisons to NaN are always false, even the test x==x is false, which is a neat way to detect that x==NaN by the way. NaN comparisons in loop conditions can lead to early termination (probably good) or non-termination (probably bad). In the VB code there is a WHILE instead of an UNTIL so NaN is "good" instead of "bad":
Loop While errval >= tol And k <= maxlevel

These kind of issues happen more often than we realize (so I'm hopefully not embarrassing anyone to point this out, because it happens a lot). Because FP is generally hard to use, a while ago I wrote a couple of lecture notes on the subject for graduate students. It's the usual fare on FP and SIMD and there is much more that I would like to add, in particular the reading material "what every computer scientist should know about floating point arithmetic".

If you're in to this kind of low level detail and there is not much to do because of the pandemic, then perhaps this is a text to get excited about!

- Rob