The precision of `gausslegendre` can be improved

[hyrodium]

MWE:

julia> using FastGaussQuadrature, LinearAlgebra, Plots

julia> errors = Float64[]
Float64[]

julia> for n in 6:100
           nodes, weights = gausslegendre(n)
           push!(errors, norm([dot(weights, nodes.^p) - 2/(1+p) for p in 2:2:10], Inf))
       end

julia> plot(6:100, errors)

The precision with $n\le60$ is worse than with $n>60$, and FastGaussQuadrature.rec can be improved.

FastGaussQuadrature.jl/src/gausslegendre.jl

Lines 46 to 52 in 03b6fb5

	elseif n ≤ 60
	# NEWTON'S METHOD WITH THREE-TERM RECURRENCE:
	return rec(n)
	else
	# USE ASYMPTOTIC EXPANSIONS:
	return asy(n)
	end

[hyrodium]

At least, it seems that the threshold 60 can be replaced with 40.

julia> using FastGaussQuadrature, LinearAlgebra, Plots

julia> errors = Float64[]
Float64[]

julia> for n in 6:100
           nodes, weights = gausslegendre(n)
           push!(errors, norm([dot(weights, nodes.^p) - 2/(1+p) for p in 2:2:10], Inf))
       end

julia> plot(6:100, errors)

julia> errors_asy = Float64[]
Float64[]

julia> for n in 6:100
           nodes, weights = FastGaussQuadrature.asy(n)
           push!(errors_asy, norm([dot(weights, nodes.^p) - 2/(1+p) for p in 2:2:10], Inf))
       end

julia> plot!(6:100, errors_asy, ylims=(0, 2e-15))