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