gauss incorrect for large n

[danielwe]

gauss returns incorrect values for n > 1032 or so. The error increases gradually with n and is way off by n = 1080. This is not Float64 specific: I've confirmed that the same incorrect values are produced by gauss(Double64, n) with Double64 from DoubleFloats.jl. Not sure if it could be related to the accuracy issues for kronrod mentioned here: #93 (comment)

Here's an MWE using FastGaussQuadratures as reference:

julia> using FastGaussQuadrature, QuadGK
julia> qx, qw = gauss(1080);
julia> fx, fw = gausslegendre(1080);
julia> norm(qx - fx)
1.2449820836798886

julia> norm(qw - fw)
0.005887119847361489

Scatter plot of gauss nodes against FastGaussQuadrature nodes:

julia> using CairoMakie
julia> scatter(fx, qx)

Using quadgk to confirm that the error is in gauss, not in FastGaussQuadrature:

julia> using LinearAlgebra
julia> dot(exp.(qx), qw)  # QuadGK.gauss
2.346053485734814

julia> dot(exp.(fx), fw)  # FastGaussQuadrature.gausslegendre
2.3504023872876028

julia> quadgk(exp, -1, 1)[1]  # Ground truth
2.3504023872876028

[stevengj]

Thanks for catching this; it looks like there was an overflow/underflow problem in my eigensolver for large n (as well as a bug where it was doing 1000x more iterations than needed). Should be fixed by #111.