Polynomials.fit fails with InexactError for Rational

[agostonsipos]

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.

[jverzani]

Hmm, maybe we can make

R = typeof(one(T)/1)

Which will make Integer go to Float, but keep Rational as is? The underlying solve_vander only needs to have the type R preserved under /.

Would you want to make a PR and see if all the tests run?

[agostonsipos]

I like the idea! Yes I will try doing a PR.