spurious underflow in kronrod for large n

[stevengj]

As mentioned on discourse, the kronrod routine yields NaN for larger n, e.g.

kronrod(550,1e-11,40)

It looks like this is due to spurious underflow in the _kronrodjacobi function. It forms arrays s and t by calculations that involve repeatedly multiplying by the b matrix, corresponding to the squared off-diagonal elements of the Jacobi matrix. For the canonical interval (-1,1) with a constant weight function, the entries of the b matrix asymptote to 0.25, so this effectively divides by 4 over and over and eventually underflows.

The underflow is spurious because the final Kronrod–Jacobi matrix only depends on ratios of s and t entries, so we are free to rescale them by any overall constant factor that we want. This fact is actually mentioned in a somewhat obscure way by Laurie's paper:

Actually, as recommended by Reichel [23], the coefficients were scaled to the interval (−2, 2) before performing the computation, and the resulting points and weights (obtained by Gautschi’s routine gauss from [10]) were scaled back to (−1, 1). On a machine with a small exponent range, this subterfuge is necessary to avoid underflow of the modified moments.

Rescaling the interval to (-2,2) works because it multiplies the integrals by 2, and hence the squared Jacobi matrix entries get multiplied by 4, causing b to asymptote to 1 and removing the exponential shrinkage.

We could do the same trick, but it's not clear how to generalize this to arbitrary Jacobi matrices. Instead, it should suffice to normalize the s and t matrices (e.g. scaling by their infinity norm) after each iteration. This should only slow things down by a small constant factor (1–2% in my current experiments).

[stevengj]

Even with the rescaling, however, Laurie's algorithm seems to suffer from some severe accuracy issues in Float64 precision for n ≳ 520 or so. (Because of that, this underflow might be a blessing in disguise.)

Update: this seems to be fixed by #111 (i.e. it's my fault, not Laurie's).