Loosen limit n ≤ 4000 in gaussjacobi

[dlfivefifty]

This seems to be a hold over from MATLAB since

julia> n=5000;a=rand(n);b=rand(n-1);@time eig(SymTridiagonal(a,b));
  4.166163 seconds (276.39 k allocations: 203.424 MB, 1.36% gc time)

[hyrodium]

Is there any technical problem for this limitations?

Is it okay to replace

    elseif n <= 4000 && max(a,b)>=5.
        jacobi_gw(n, a, b)
    else
        error("gaussjacobi($n,$a,$b) is not implemented: n must be ≤ 4000 for max(a,b)≥5.")
    end

with

    else
        jacobi_gw(n, a, b)
    end

?

https://github.com/JuliaApproximation/FastGaussQuadrature.jl/blob/master/src/gaussjacobi.jl#L25

If there's a problem, I guess the convergence of Newton's method is suspicious.
Because the Newton iteration algorithm has hard-coded loop limitation.
https://github.com/JuliaApproximation/FastGaussQuadrature.jl/blob/master/src/gaussjacobi.jl#L88

[dlfivefifty]

Go for it! I think the limitation on n ≤ 4000 was just because Matlab was not using SymTridiagonal.

If there's a problem, I guess the convergence of Newton's method is suspicious

True, this could effect the large a and b case. Even worse, the initial guesses might not be sufficiently accurate to converge. @ajt60gaibb any thoughts?

[dlfivefifty]

If we add GenericLinearAlgebra.jl as a dependency we could just use Golub-Welsch with BigFloat for large a and b and then cast to Float64.