[JacobHass8]

I'm working on adding more Gaussian quadrature methods using Laguerre, Associated Laguerre and Jacobi polynomials. This is done by generalizing the orthogonal polynomial root finder in #1461 to any family of polynomials. I also added tests for finding the root of these polynomials in the process.

So far I've only implemented everything for Laguerre polynomials, but I'm going to start on Associate Laguerre polynomials next.

[NAThompson]

@JacobHass8 : I encourage you to take on the adaptive versions of these in the first pass-or at least give an error estimate. This truly increases the value of these methods to users.

[JacobHass8]

@JacobHass8 : I encourage you to take on the adaptive versions of these in the first pass-or at least give an error estimate. This truly increases the value of these methods to users.

I deleted my earlier comment because I realized adaptive methods are much, much harder for Hermite and Laguerre polynomials than simply using Gauss-Kronrod. I went back to your previous comment and took a look at the paper here: https://arxiv.org/html/1203.1554v2. Could I just use the Mathematica code provided to generate the abscissa and weights for fixed number of points? I can also work with Claude to translate the code but I'm not sure how long that would take. It also needs to solve a linear system of equations so will require Eigen or a linear solver.

As you noted earlier, an error estimate could just be provided by evaluating a higher order method and taking the difference. This is very minimal but might not be the best estimation.

Do you have any preference on which route to go?

[NAThompson]

@JacobHass8 : Hmm, I think you're right-the literature is considerably less developed for this problem.

What about tackling Genz-Kaister first?

[JacobHass8]

@JacobHass8 : Hmm, I think you're right-the literature is considerably less developed for this problem.

What about tackling Genz-Kaister first?

It looks like Genz-Keister will be difficult/impossible to compute on the fly. Apparently the linear system of equations can be ill-conditioned and return complex roots. I think that I'll just use the mathematica code in the arxiv paper to calculate the coefficients to higher precision

[NAThompson]

@JacobHass8 : Could you get your quadrature methods added to the nanobind layer so I can play around with them?

[JacobHass8]

@JacobHass8 : Could you get your quadrature methods added to the nanobind layer so I can play around with them?

Yes, I can work on that next. I've only really implemented Gauss-Laguerre though. I still need to do Associated Laguerre and Jacobi methods.

I also have weights and abscissa for the nested hermite method. I still need to test if they actually give higher accuracy than regular gauss hermite. I assumed this would go in a separate PR but I can commit it here if you want.

[NAThompson]

I assumed this would go in a separate PR but I can commit it here if you want.

Separate PR is fine.

[NAThompson]

@JacobHass8 : I have removed the protected from the member function and allow the nodes and weights to be returned as const. These integrators are not only useful for integrating functions, but if the nodes and weights are available they can be used to solve integral equations.

[NAThompson]

@JacobHass8 : If we integrate any polynomial ≤ 2d, then we should get error on the order of the unit roundoff

If we use the orthogonal polynomial family associated with the weight, then if we go beyond the degree then we should get error of order 1. This is what we observe.

[NAThompson]

Another one: If the Gauss estimate is an overestimate, the anti-Gauss should be an underestimate, and vice versa:

Note that this is a "symlog" scale-it's actually extremely accurate.

[JacobHass8]

@JacobHass8 : I have removed the protected from the member function and allow the nodes and weights to be returned as const. These integrators are not only useful for integrating functions, but if the nodes and weights are available they can be used to solve integral equations.

I think that's a good idea. I also added hermite_zeroes, laguerre_zeroes and jacobi_zeroes to the special functions namespace.

[NAThompson]

@JacobHass8 : I have also played around with the functions through the nanobind layer:

>>> import boost.math as bm, math
>>> bm.gauss_hermite(lambda x : x*x, points=10)
(0.8862269254527579, 3.935638150721656e-16, 0.8862269254527579)
>>> math.sqrt(math.pi)/2
0.8862269254527579

[JacobHass8]

@JacobHass8 : If we integrate any polynomial ≤ 2d, then we should get error on the order of the unit roundoff

If we use the orthogonal polynomial family associated with the weight, then if we go beyond the degree then we should get error of order 1. This is what we observe.

Could you also generate these plots for Laguerre and Jacobi? These are really nice!

[JacobHass8]

I wonder if our performance is any better than numpy https://numpy.org/doc/2.1/reference/generated/numpy.polynomial.hermite.hermgauss.html

[NAThompson]

@JacobHass8 : Yes, Laguerre and Jacobi are done too. Everything is in one HDF5 file: gauss_correctness.h5.zip. To view it, install the H5Web fork's VS Code extension (h5web-fork-4d2ac35.vsix; in VS Code, Extensions → … → Install from VSIX) and open the .h5.

What's covered, for N = 1..40 (Laguerre from N = 2):

Plot	Hermite	Laguerre	Jacobi
Discrete orthogonality |Σ wᵢ pⱼ(xᵢ)pₖ(xᵢ) − δⱼₖ|/eps, sliced by N	✓	α = 0, 0.5, −0.5	(α,β) = (0,0), (0.5,−0.5), (−0.5,−0.5), (2,1)
Exactness staircase vs degree d	✓	same	same
Gauss vs anti-Gauss mirror, f = e^{ωx}, ω = 3, 8	—	—	(0,0), (0.5,−0.5), (2,1)

The pⱼ are orthonormal polynomials from an independent recurrence evaluated in cpp_bin_float_50. Every family is exact (≤ 64 eps) up to degree 2N−1, and the first failure is exactly at d = 2N for every N. The anti-Gauss reference is a closed form via hypergeometric_1F1, checked against tanh-sinh to 1e-49. The anti-Gauss plots are Jacobi only: the Laguerre anti-Gauss nodes and weights are private, so they can't be reached.

I also pushed nanobind bindings for gauss_laguerre(f, points, alpha) and a gauss_jacobi(n, alpha, beta) class exposing abscissa, weights, anti_gauss_abscissa, anti_gauss_weights and anti_gauss_is_internal.

[NAThompson]

I wonder if our performance is any better than numpy https://numpy.org/doc/2.1/reference/generated/numpy.polynomial.hermite.hermgauss.html

@JacobHass8 : This turned out to be a fruitful question. . . one that admittedly made the MR much longer lol.

[JacobHass8]

I wonder if our performance is any better than numpy https://numpy.org/doc/2.1/reference/generated/numpy.polynomial.hermite.hermgauss.html

@JacobHass8 : This turned out to be a fruitful question. . . one that admittedly made the MR much longer lol.

I'm not familiar with MR... merge request?

[NAThompson]

@JacobHass8 : Apologies they call it "MR" in gitlab-got mixed up.

Before your question this was both slower and less accurate that scipy/numpy. Then we switched to this paper and now we're way faster and way more accurate. In addition, we can compute Gauss and anti-Gauss rules out to way higher N.

If you've installed my VSCode fork I've done this for every rule so you can click around if you wish.

[NAThompson]

@JacobHass8 : I'm a bit worried about the changes to special_functions/{hermite,laguerre,jacobi}.hpp. They don't change any existing function, but every user of hermite(), laguerre(), jacobi() (and special_functions.hpp, hypergeometric_2F0.hpp) would now also compile the root finder, plus gamma/digamma/beta. hermite.hpp already needed an NVRTC guard because of it. These headers are also used by MSVC's STL for std::hermite / std::laguerre (stl/src/special_math.cpp includes them via its boost-math submodule), so I'd rather not grow them at all.

I'd like to leave those three headers exactly as on develop, and move the *_family classes into a quadrature-side detail header that only the gauss_* headers include. Would you be OK with splitting hermite_zeros/laguerre_zeros/jacobi_zeros into a follow-up PR? Or we could keep them here in a new header of their own.

[JacobHass8]

@JacobHass8 : Apologies they call it "MR" in gitlab-got mixed up.

Before your question this was both slower and less accurate that scipy/numpy. Then we switched to this paper and now we're way faster and way more accurate. In addition, we can compute Gauss and anti-Gauss rules out to way higher N.

If you've installed my VSCode fork I've done this for every rule so you can click around if you wish.

Wow that's a substantial improvement! I was mainly doing this to have a C++ implementation but maybe the scipy/numpy maintainers would be interested now.

[JacobHass8]

I'd like to leave those three headers exactly as on develop, and move the *_family classes into a quadrature-side detail header that only the gauss_* headers include. Would you be OK with splitting hermite_zeros/laguerre_zeros/jacobi_zeros into a follow-up PR? Or we could keep them here in a new header of their own.

That sounds good to me! I wasn't aware that hermite.hpp was so huge. I was thinking it would be helpful now that we have a way of computing the zeros now. I can work on the follow up

[NAThompson]

@JacobHass8 : I now want to merge this.

Could you give it a final read?

[JacobHass8]

@JacobHass8 : I now want to merge this.

Could you give it a final read?

Sure thing! I'll take a look tonight