I was trying to use Polynomials for Base.Rational, and encountered an issue when running something like this:

xs :: Vector{Rational} = 1:21
ys :: Vector{Rational} = map(x->23//31*x^10-181//87*x^5+543//211, xs)
p = Polynomials.fit(Polynomial{Rational}, xs, ys)

The expected behavior is that a Polynomial with Rational coefficients is returned, and all steps are done with exact arithmetic.

Instead I get this error.

ERROR: InexactError: Rational(1.5394101672544472e-18)

Note that for simpler tasks, a Polynomial is returned, but inexact arithmetic is being done in the process, just there is a smaller floating point error, so it can be converted back.

Doing some investigation, the issue is with the following function in standard-basis.jl:

function _polynomial_fit(P::Type{<:StandardBasisPolynomial}, x::AbstractVector{T}, y; var=:x) where {T}
    R = float(T)
    coeffs = Vector{R}(undef, length(x))
    copyto!(coeffs, y)
    solve_vander!(coeffs, x)
    P(R.(coeffs), var)
end

Changing the first line to R = T makes my example work correctly.

julia> p = Polynomials.fit(Polynomial{Rational}, xs, ys)
Polynomial(543//211 - 181//87*x^5 + 23//31*x^10)

Of course that is not the right solution, I'd recommend something like

if T <: Rational
    R = T
else
    R = float(T)
end

not sure though if it should be more general than that though.