Close Issue 268 (#269)

* close #268; 
* add faster vander_solve function for one case

Author: jverzani

src/common.jl

@@ -75,16 +75,9 @@ function fit(P::Type{<:AbstractPolynomial},
              x::AbstractVector{T},
              y::AbstractVector{T},
              deg::Integer = length(x) - 1;
-    weights = nothing,
-    var = :x,) where {T}
-    x = mapdomain(P, x)
-    vand = vander(P, x, deg)
-    if weights !== nothing
-        coeffs = _wlstsq(vand, y, weights)
-    else
-        coeffs = pinv(vand) * y
-    end
-    return P(T.(coeffs), var)
+             weights = nothing,
+             var = :x,) where {T}
+    _fit(P, x, y, deg, weights, var)
 end
 
 fit(P::Type{<:AbstractPolynomial},
@@ -108,9 +101,27 @@ fit(x::AbstractVector,
     weights = nothing,
     var = :x,) = fit(Polynomial, x, y, deg; weights = weights, var = var)
 
+function _fit(P::Type{<:AbstractPolynomial},
+             x::AbstractVector{T},
+             y::AbstractVector{T},
+             deg::Integer = length(x) - 1;
+             weights = nothing,
+             var = :x,) where {T}
+    x = mapdomain(P, x)
+    vand = vander(P, x, deg)
+    if weights !== nothing
+        coeffs = _wlstsq(vand, y, weights)
+    else
+        coeffs = pinv(vand) * y
+    end
+    R = float(T)
+    return P(R.(coeffs), var)
+end
+
+
 # Weighted linear least squares
 _wlstsq(vand, y, W::Number) = _wlstsq(vand, y, fill!(similar(y), W))
-_wlstsq(vand, y, W::AbstractVector) = _wlstsq(vand, y, diagm(0 => W))
+_wlstsq(vand, y, W::AbstractVector) = _wlstsq(vand, y, Diagonal(W))
 _wlstsq(vand, y, W::AbstractMatrix) = (vand' * W * vand) \ (vand' * W * y)
 
 """

src/polynomials/standard-basis.jl

@@ -301,6 +301,51 @@ function vander(P::Type{<:StandardBasisPolynomial}, x::AbstractVector{T}, n::Int
 end
 
 
+## as noted at https://github.com/jishnub/PolyFit.jl, using method from SpecialMatrices is faster
+## https://github.com/JuliaMatrices/SpecialMatrices.jl/blob/master/src/vandermonde.jl
+## This is Algorithm 2 of https://www.maths.manchester.ac.uk/~higham/narep/narep108.pdf
+## Solve V(αs)⋅x = y where V is (1+n) × (1+n) Vandermonde matrix (Vᵀ in the paper)
+function solve_vander!(ys, αs)
+    n = length(ys) - 1
+    for k in 0:n-1
+        for j in n:-1:k+1
+            ys[1+j] = (ys[1+j] - ys[1+j-1])/(αs[1+j] - αs[1+j-k-1])
+        end
+    end
+    for k in n-1:-1:0
+        for j in k:n-1
+            ys[1+j] = ys[1+j] - αs[1+k] * ys[1 + j + 1]
+        end
+    end
+
+    nothing
+end
+
+# intercept one (typical) case for a faster variant
+function fit(P::Type{<:StandardBasisPolynomial},
+             x::AbstractVector{T},
+             y::AbstractVector{T},
+             deg::Integer = length(x) - 1;
+             weights = nothing,
+             var = :x,) where {T}
+
+    if deg == length(x) -1 && weights == nothing
+        _polynomial_fit(P, x, y; var=var)
+    else
+        _fit(P, x, y, deg; weights=weights, var=var)
+    end
+end
+
+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
+
+
+
 ## --------------------------------------------------
 """
     compensated_horner(p::P, x)

test/StandardBasis.jl

@@ -306,6 +306,10 @@ end
         # issue #214 --  should error
         @test_throws MethodError fit(Polynomial, rand(2,2), rand(2,2))
 
+        # issue #268 -- inexacterror
+        @test fit(P, 1:4, 1:4, var=:x) ≈ variable(P{Float64}, :x)
+        @test fit(P, 1:4, 1:4, 1, var=:x) ≈ variable(P{Float64}, :x)
+
     end
 
     # test default   (issue  #228)