GaussLaguerreIntegration may return NaN for large orders

GaussianQuadrature's constructor forms each weight as GaussianQuadrature's constructor forms each weight as

// ql/math/integrals/gaussianquadratures.cpp
w_[i] = mu_0*ev[0][i]*ev[0][i] / orthPoly.w(x_[i]);

For Gauss-Laguerre the weight function is w(x) = x^s exp(-x) and the largest node grows roughly as 4n, so beyond a certain order w(x_max) underflows to zero. The division then yields inf, and because the integrand at that node is itself of order exp(-x), operator() evaluates inf * 0 and the entire integral comes back NaN.

The quotient itself is finite and of order one at every node — the squared eigenvector component underflows just as fast as w(x) does. Only the intermediate w(x_i) leaves the representable range.

Reproducer

#include <ql/math/integrals/gaussianquadratures.hpp>
#include <cmath>
#include <iostream>
using namespace QuantLib;

int
main() {
       for (Size n = 184; n <= 232; n += 8) {
              GaussLaguerreIntegration q(n);
              // integral of exp(-x) over [0, inf) is exactly 1
              std::cout << "n = " << n << "   result = "
                  << q([](Real x) { return std::exp(-x); }) << std::endl;
       }
}

Expected 1 throughout. Observed: correct up to the threshold, then nan.

Threshold

The first non-finite weight appears at n = 200 (x_max = 767.81, where exp(-767.81) is below even the smallest subnormal). If subnormals are flushed to zero the threshold drops to n = 192, which is the order used by test-suite/hestonmodel.cpp:1774 and :2721 (AnalyticHestonEngine::Integration::gaussLaguerre(192)). AnalyticHestonEngine's default integrationOrder of 144 is well inside the safe range.

Below the threshold there is no accuracy penalty: the relative error of sum w_i exp(-x_i) x_i^k against k! is ~1e-15 either way, for every order where the current code produces finite weights. This is a range problem, not a precision one.

Shape of a fix

Forming the quotient in logarithms keeps every intermediate in range:

const Real wx = orthPoly.w(x_[i]);
if (wx > 0.0) {
       w_[i] = mu_0*ev[0][i]*ev[0][i] / wx;
} else {
       w_[i] = std::exp(std::log(mu_0)
           + 2.0*std::log(std::abs(ev[0][i]))
           - orthPoly.logW(x_[i]));
}

This needs a logW() alongside the existing w() on GaussianOrthogonalPolynomial — s*log(x) - x for Laguerre, 2*mu*log(|x|) - x*x for Hermite. With it the reproducer returns 1 at every order, and sub-threshold results are bit-identical to the current code.