Merge pull request #39 from pwl/master

Resolves #38 by adding a `trunctate` function

Author: jverzani

src/Polynomials.jl

@@ -9,10 +9,11 @@ export Poly, poly
 export degree, coeffs
 export polyval, polyint, polyder, roots, polyfit
 export Pade, padeval
+export truncate!
 
 import Base: length, endof, getindex, setindex!, copy, zero, one, convert, norm, gcd
 import Base: show, print, *, /, //, -, +, ==, divrem, div, rem, eltype
-import Base: promote_rule
+import Base: promote_rule, truncate
 if VERSION >= v"0.4"
     import Base.call
 end
@@ -67,7 +68,7 @@ immutable Poly{T<:Number}
             return new(zeros(T,1), symbol(var))
         else
             # determine the last nonzero element and truncate a accordingly
-            a_last = max(1,findlast( p->(abs(p) > 2*eps(T)), a))
+            a_last = max(1,findlast(a))
             new(a[1:a_last], symbol(var))
         end
     end
@@ -136,6 +137,18 @@ coeffs(p::Poly) = p.a
 
 """
 
+`truncate{T}(p::Poly{T}; reltol = eps(T), abstol = eps(T))`: returns a polynomial with coefficients a_i truncated to zero if |a_i| <= reltol*maxabs(a)+abstol
+
+"""
+function truncate{T}(p::Poly{T}; reltol = eps(T), abstol = eps(T))
+    a = coeffs(p)
+    amax = maxabs(a)
+    anew = map(ai -> abs(ai) <= amax*reltol+abstol ? zero(T) : ai, a)
+    return Poly(anew)
+end
+
+"""
+
 * `norm(q::Poly, [p])`: return `p` norm of polynomial `q`
 
 """
@@ -253,7 +266,7 @@ function divrem{T, S}(num::Poly{T}, den::Poly{S})
 
     return pQ, pR
 end
-                  
+
 @deprecate /(num::Poly, den::Poly)                div(num::Poly, den::Poly)
 div(num::Poly, den::Poly) = divrem(num, den)[1]
 rem(num::Poly, den::Poly) = divrem(num, den)[2]
@@ -266,7 +279,7 @@ function ==(p1::Poly, p2::Poly)
     end
 end
 
-"""                  
+"""
 * `polyval(p::Poly, x::Number)`: Evaluate the polynomial `p` at `x` using Horner's method.
 
 Example:
@@ -314,7 +327,7 @@ polyint(Poly([1, 0, -1]))     # Poly(x - 0.3333333333333333x^3)
 polyint(Poly([1, 0, -1]), 2)  # Poly(2.0 + x - 0.3333333333333333x^3)
 ```
 
-"""                  
+"""
 function polyint{T}(p::Poly{T}, k::Number=0)
     n = length(p)
     R = typeof(one(T)/1)
@@ -379,7 +392,7 @@ Examples:
 roots(Poly([1, 0, -1]))    # [-1.0, 1.0]
 roots(Poly([1, 0, 1]))     # [0.0+1.0im, 0.0-1.0im]
 roots(Poly([0, 0, 1]))     # [0.0, 0.0]
-roots(poly([1,2,3,4]))     # [1.0,2.0,3.0,4.0] 
+roots(poly([1,2,3,4]))     # [1.0,2.0,3.0,4.0]
 ```
 """
 function roots{T}(p::Poly{T})
@@ -461,12 +474,12 @@ Original by [ggggggggg](https://github.com/Keno/Polynomials.jl/issues/19)
 function polyfit(x, y, n::Int=length(x)-1, sym::Symbol=:x)
     length(x) == length(y) || throw(DomainError)
     1 <= n <= length(x) - 1 || throw(DomainError)
-    
+
     A = [ float(x[i])^p for i = 1:length(x), p = 0:n ]
     Poly(A \ y, sym)
 end
 polyfit(x,y,sym::Symbol) = polyfit(x,y,length(x)-1, sym)
-   
+
 ### Pull in others
 include("pade.jl")
 

test/runtests.jl

@@ -55,7 +55,7 @@ sprint(show, pNULL)
 if VERSION >= v"0.4"
     @test pN(-.125) == 276.9609375
     @test pN([0.1, 0.2, 0.3]) == polyval(pN, [0.1, 0.2, 0.3])
-end 
+end
 
 @test poly([-1,-1]) == p3
 @test roots(p0)==roots(p1)==roots(pNULL)==[]
@@ -141,3 +141,15 @@ p = polyfit(xs, ys, :t)
 p = polyfit(xs, ys, 2)
 @test maximum(abs(map(x->polyval(p, x), xs) - ys)) <= 0.03
 
+
+## truncation
+p1 = Poly([1,1]/10)
+p2 = Poly([1,2]/10)
+p3 = Poly([1,3]/10)
+psum = p1 + p2 - p3
+@test degree(psum) == 1         # will have wrong degree
+@test degree(truncate(psum)) == 0 # the degree should be correct after truncation
+
+@test truncate(Poly([2,1]),reltol=1/2,abstol=0) == Poly([2])
+@test truncate(Poly([2,1]),reltol=1,abstol=0)   == Poly([0])
+@test truncate(Poly([2,1]),reltol=0,abstol=1)   == Poly([2])