Merge pull request #143 from jverzani/v0.7

degree of zero poly; fix v0.7 warnings

Author: jverzani

src/Polynomials.jl

@@ -15,8 +15,8 @@ export Pade, padeval
 
 import Compat.lastindex
 
-import Base: start, next, done, length, size, eltype, collect, eachindex
-import Base:  getindex, setindex!, copy, zero, one, convert, gcd
+import Base: length, size, eltype, collect, eachindex
+import Base: getindex, setindex!, copy, zero, one, convert, gcd
 import Base: show, print, *, /, //, -, +, ==, isapprox, divrem, div, rem, eltype
 import Base: promote_rule, truncate, chop,  conj, transpose, hash
 import Base: isequal
@@ -146,10 +146,16 @@ eltype(::Type{Poly{T}}) where {T} = Poly{T}
 length(p::Poly) = length(coeffs(p))
 lastindex(p::Poly)  = length(p) - 1
 
-start(p::Poly)        = start(coeffs(p)) - 1
-next(p::Poly, state)  = (temp = fill!(similar(coeffs(p)), 0); temp[state+1] = p[state]; (Poly(temp), state+1))
-done(p::Poly, state)  = state > degree(p)
-
+if VERSION < v"0.7.0-DEV.5126"
+    import Base: start, next, done
+    start(p::Poly)        = start(coeffs(p)) - 1
+    next(p::Poly, state)  = (temp = fill!(similar(coeffs(p)), 0); temp[state+1] = p[state]; (Poly(temp), state+1))
+    done(p::Poly, state)  = state > degree(p)
+else
+    import Base: iterate
+    Base.iterate(p::Poly) = (p[0] * one(p), 1)
+    Base.iterate(p::Poly, state) = state <= degree(p) ? (p[state]*variable(p)^(state), state+1) : nothing
+end
 
 # shortcut for collect(eltype, collection)
 collect(p::Poly{T}) where {T} = collect(Poly{T}, p)
@@ -161,9 +167,9 @@ size(p::Poly, i::Integer) = size(p.a, i)
     degree(p::Poly)
 
 Return the degree of the polynomial `p`, i.e. the highest exponent in the polynomial that
-has a nonzero coefficient.
+has a nonzero coefficient. The degree of the zero polynomial is defined to be -1.
 """
-degree(p::Poly) = length(p) - 1
+degree(p::Poly) = iszero(p) ? -1 : length(p) - 1
 
 """
     coeffs(p::Poly)
@@ -244,7 +250,7 @@ the `p`-norm is
 ||q||_p = (|q_0|^p + \\ldots + |q_n|^p)^{1/p}
 ``
 """
-norm(q::Poly, args...) = norm(coeffs(q), args...)
+norm(q::Poly, p::Real=2) = norm(coeffs(q), p)
 
 
 """
@@ -270,7 +276,7 @@ function setindex!(p::Poly, value, idx::Int)
     n = length(p.a)
     if n ≤ idx
         resize!(p.a, idx+1)
-        p.a[n+1:idx] = 0
+        p.a[n+1:idx] .= 0
     end
     p.a[idx+1] = value
     return p
@@ -352,6 +358,9 @@ function *(p1::Poly{T}, p2::Poly{S}) where {T,S}
     Poly(a,p1.var)
 end
 
+# quiet this deprecation https://github.com/JuliaLang/julia/pull/23332
+import Base: ^
+^(p::Poly, n::Integer) = Base.power_by_squaring(p,n)
 
 # are any values NaN
 hasnan(p::Poly) = reduce(|, (isnan.(p.a)))
@@ -706,17 +715,17 @@ function polyfit(x, y, n::Int=length(x)-1, sym::Symbol=:x)
     # here unsure, whether similar(float(x[1]),...), or similar(x,...)
     # however similar may yield unwanted surprise in case of e.g. x::Int
     #
-    A=similar(float.(x[1:1]), length(x), n+1)
+    T = eltype(float(x[1]))
+    A = Array{T}(undef, length(x), n+1)
     #
     # TODO: add support for poly coef bitmap
     # (i.e. polynomial with some terms fixed to zero)
     #
-    A[:,1]=1
+    A[:,1] .= 1
     for i=1:n
-        A[:,i+1]=A[:,i] .* x   # cumulative product more precise than x.^n
+        A[:,i+1] .= A[:,i] .* x   # cumulative product more precise than x.^n
     end
-    Aqr=qrfact(A)   # returns QR object, not a matrix
-    p=Aqr\y         # least squares solution via QR
+    p = A \ float.(y)
     Poly(p, sym)
 end
 polyfit(x,y,sym::Symbol) = polyfit(x,y,length(x)-1, sym)

test/REQUIRE

@@ -0,0 +1 @@
+SpecialFunctions

test/runtests.jl

@@ -3,6 +3,7 @@ using Compat
 using Compat.Test
 using Compat.LinearAlgebra
 using Polynomials
+using SpecialFunctions
 
 import Compat.SparseArrays: sparse, speye, nnz
 
@@ -46,7 +47,7 @@ sprint(show, pNULL)
 @test pNULL^3 == pNULL
 @test pNULL*pNULL == pNULL
 
-@test map(degree, [pNULL,p0,p1,p2,p3,p4,p5,pN,pR,p1000]) == [0,0,0,1,2,3,4,4,2,999]
+@test map(degree, [pNULL,p0,p1,p2,p3,p4,p5,pN,pR,p1000]) == [-1,-1,0,1,2,3,4,4,2,999]
 
 @test polyval(poly(Int[]), 2.) == 1.
 @test polyval(pN, -.125) == 276.9609375
@@ -168,7 +169,9 @@ p = polyfit(xs, ys)
 p = polyfit(xs, ys, :t)
 p = polyfit(xs, ys, 2)
 @test maximum(map(abs,map(x->polyval(p, x), xs) - ys)) <= 0.03
-
+#https://stackoverflow.com/questions/50832823/error-with-polyfit-function-julia
+# relax type assumptions on x, y
+polyfit(Any[1,2,3], Any[2,3,1])
 
 ## truncation
 p1 = Poly([1,1]/10)
@@ -225,8 +228,8 @@ p1[0:1] = [7,8]
 ## conjugate of poly (issue #59)
 as = [im, 1, 2]
 bs = [1, 1, 2]
-@test conj(Poly(as)) == Poly(conj(as))
-@test conj(Poly(bs)) == Poly(conj(bs))
+@test conj(Poly(as)) == Poly(conj.(as))
+#@test conj(Poly(bs)) == Poly(conj(bs)) # conj gives warning as no defn on T
 ## and transpose
 @test transpose(Poly(as)) == Poly(as)
 
@@ -355,7 +358,7 @@ p2s = Poly([1], :s)
 # test size
 @test size(Poly([0.5, 0.2])) == (2,)
 @test size(Poly([0.5, 0.2]), 1) == 2
-@test size(Poly([0.5, 0.2]), 1, 2) == (2,1)
+# @test size(Poly([0.5, 0.2]), 1, 2) == (2,1) # deprecated in v0.7
 
 # test iteration
 p1 = Poly([1,2,0,3])
@@ -393,3 +396,9 @@ p1 = poly([1.,2.,3.])
 p2 = poly([1.,2.,6.])
 
 @test (res = roots(gcd(p1, p2)); 1. ∈ res && 2. ∈ res)
+
+## Getting error on passing Real arrays to polyfit #146
+xx = Real[20.0, 30.0, 40.0]
+yy = Real[15.7696, 21.4851, 28.2463]
+polyfit(xx,yy,2)
+