Solve performance regresion (#103)

* Use `fixorder` again, to avoid the if's used in `getindex`

The use of those if's is at the heart of the problem reported in [#102].
The idea was to avoid fixorder by setting the coefficients that are
not defined to zero. This involves if's that killed performance.

Also, fixorder was corrected to return a copy of the coefficients to
avoid side effects.

`max_order` is not anylonger needed, so it is deleted.

* Fix tests and fateman.jl

* Tiny correction to fixorder

* Improvements on allocations

* Inline mul! and div! methods

* Use a modified power_by_squaring, and inline the (mutating) recursion functions

* Inline fixorder

* Inlining div!

perf/fateman.jl

@@ -10,9 +10,9 @@ function fateman1(degree::Int)
     s = TaylorN( [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree )
     s = s^degree
     # s is converted to order 2*ndeg
-    s = TaylorN(s, 2*degree)
+    s = TaylorN(s.coeffs, 2*degree)
 
-    s * ( s+TaylorN(oneH, 2*degree) )
+    s * ( s+1 )
 end
 
 function fateman2(degree::Int)
@@ -22,7 +22,7 @@ function fateman2(degree::Int)
     s = TaylorN( [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree )
     s = s^degree
     # s is converted to order 2*ndeg
-    s = TaylorN(s, 2*degree)
+    s = TaylorN(s.coeffs, 2*degree)
     return s^2 + s
 end
 

src/TaylorSeries.jl

@@ -32,7 +32,7 @@ import Base: zero, one, zeros, ones, isinf, isnan,
     rem, mod, mod2pi, abs, gradient,
     sqrt, exp, log, sin, cos, tan,
     asin, acos, atan, sinh, cosh, tanh,
-    A_mul_B!,
+    A_mul_B!, power_by_squaring, 
     getindex, setindex!, endof
 
 export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries

src/arithmetic.jl

@@ -15,18 +15,10 @@ for T in (:Taylor1, :TaylorN)
         =={T<:Number,S<:Number}(a::$T{T}, b::$T{S}) = ==(promote(a,b)...)
 
         function =={T<:Number}(a::$T{T}, b::$T{T})
-            la = a.order+1
-            lb = b.order+1
-            if a.order == b.order
-                return all( a.coeffs .== b.coeffs )
-            elseif a.order < b.order
-                res1 = all( a[1:la] .== b[1:la] )
-                res2 = iszero(b[la+1:lb])
-            else a.order > b.order
-                res1 = all( a[1:lb] .== b[1:lb] )
-                res2 = iszero(a[lb+1:la])
+            if a.order != b.order
+                a, b = fixorder(a, b)
             end
-            return res1 && res2
+            return a.coeffs == b.coeffs
         end
     end
 end
@@ -40,17 +32,11 @@ end
 
 ## zero and one ##
 for T in (:Taylor1, :TaylorN), f in (:zero, :one)
-    @eval begin
-        ($f)(a::$T) = $T([($f)(a[1])], a.order)
-
-        ($f)(a::$T, order) = $T([($f)(a[1])], order)
-    end
+    @eval ($f)(a::$T) = $T(($f)(a[1]), a.order)
 end
 
 function zero{T<:Number}(a::HomogeneousPolynomial{T})
-    a.order == 0 && return HomogeneousPolynomial([zero(T)], 0)
-    v = Array{T}( size_table[a.order+1] )
-    @__dot__ v = zero(a.coeffs)
+    v = zeros(a.coeffs)
     return HomogeneousPolynomial(v, a.order)
 end
 
@@ -67,9 +53,7 @@ zeros{T<:Number}(::Type{HomogeneousPolynomial{T}}, order::Int) =
     zeros( HomogeneousPolynomial([zero(T)], 0), order)
 
 function one{T<:Number}(a::HomogeneousPolynomial{T})
-    a.order == 0 && return HomogeneousPolynomial([one(a[1])], 0)
-    v = Array{T}( size_table[a.order+1] )
-    @__dot__ v = one(a.coeffs)
+    v = ones(a.coeffs)
     return HomogeneousPolynomial{T}(v, a.order)
 end
 
@@ -96,19 +80,12 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
             ($f){T<:Number,S<:Number}(a::$T{T}, b::$T{S}) = $f(promote(a,b)...)
 
             function $f{T<:Number}(a::$T{T}, b::$T{T})
-                la = a.order+1
-                lb = b.order+1
-                v = similar(a.coeffs, max(la,lb))
-                if a.order == b.order
-                    @__dot__ v = $f(a.coeffs, b.coeffs)
-                elseif a.order < b.order
-                    @inbounds @__dot__ v[1:la] = $f(a[1:la], b[1:la])
-                    @inbounds @__dot__ v[la+1:lb] = $f(b[la+1:lb])
-                else
-                    @inbounds @__dot__ v[1:lb] = $f(a[1:lb], b[1:lb])
-                    @inbounds @__dot__ v[lb+1:la] = a[lb+1:la]
+                if a.order != b.order
+                    a, b = fixorder(a, b)
                 end
-                return $T(v)
+                v = similar(a.coeffs)
+                @__dot__ v = $f(a.coeffs, b.coeffs)
+                return $T(v, a.order)
             end
 
             function $f(a::$T)
@@ -120,8 +97,7 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!)))
             ($f){T<:Number,S<:Number}(a::$T{T}, b::S) = $f(promote(a,b)...)
 
             function $f{T<:Number}(a::$T{T}, b::T)
-                coeffs = similar(a.coeffs)
-                coeffs .= a.coeffs
+                coeffs = copy(a.coeffs)
                 @inbounds coeffs[1] = $f(a[1], b)
                 return $T(coeffs, a.order)
             end
@@ -253,10 +229,11 @@ for (T, W) in ((:Taylor1, :Number), (:TaylorN, :NumberNotSeriesN))
     @eval *{T<:$W, S<:$W}(a::$T{T}, b::$T{S}) = *(promote(a,b)...)
 
     @eval function *{T<:$W}(a::$T{T}, b::$T{T})
-        a == b && return square(a)
-        corder = max(a.order, b.order)
-        c = zero(a, corder)
-        for ord = 0:corder
+        if a.order != b.order
+            a, b = fixorder(a, b)
+        end
+        c = $T(zero(a[1]), a.order)
+        for ord = 0:c.order
             mul!(c, a, b, ord) # updates c[ord+1]
         end
         return c
@@ -268,7 +245,6 @@ end
     b::HomogeneousPolynomial{S}) = *(promote(a,b)...)
 
 function *{T<:NumberNotSeriesN}(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T})
-    a == b && return square(a)
 
     order = a.order + b.order
 
@@ -280,10 +256,11 @@ function *{T<:NumberNotSeriesN}(a::HomogeneousPolynomial{T}, b::HomogeneousPolyn
 end
 
 
-# Homogeneous coefficient for the multiplication
+# Internal multiplication functions
 for T in (:Taylor1, :TaylorN)
-    @eval function mul!(c::$T, a::$T, b::$T, k::Int)
+    @eval @inline function mul!(c::$T, a::$T, b::$T, k::Int)
 
+        # c[k+1] = zero( a[k+1] )
         @inbounds for i = 0:k
             if $T == Taylor1
                 c[k+1] += a[i+1] * b[k-i+1]
@@ -317,7 +294,7 @@ c_k = \sum_{j=0}^k a_j b_{k-j}.
 Return `c = a*b` with no allocation; all arguments are `HomogeneousPolynomial`.
 
 """
-function mul!(c::HomogeneousPolynomial, a::HomogeneousPolynomial,
+@inline function mul!(c::HomogeneousPolynomial, a::HomogeneousPolynomial,
         b::HomogeneousPolynomial)
 
     (iszero(b) || iszero(a)) && return nothing
@@ -394,13 +371,15 @@ end
 /{T<:Number,S<:Number}(a::Taylor1{T}, b::Taylor1{S}) = /(promote(a,b)...)
 
 function /{T<:Number}(a::Taylor1{T}, b::Taylor1{T})
-    corder = max(a.order, b.order)
+    if a.order != b.order
+        a, b = fixorder(a, b)
+    end
 
     # order and coefficient of first factorized term
     ordfact, cdivfact = divfactorization(a, b)
 
-    c = Taylor1([cdivfact], corder)
-    for ord = 1:corder-ordfact
+    c = Taylor1(cdivfact, a.order)
+    for ord = 1:a.order-ordfact
         div!(c, a, b, ord, ordfact) # updates c[ord+1]
     end
 
@@ -412,12 +391,15 @@ end
 
 function /{T<:NumberNotSeriesN}(a::TaylorN{T}, b::TaylorN{T})
     @assert !iszero(constant_term(b))
-    corder = max(a.order, b.order)
+
+    if a.order != b.order
+        a, b = fixorder(a, b)
+    end
 
     # first coefficient
     @inbounds cdivfact = a[1] / constant_term(b)
-    c = TaylorN(cdivfact, corder)
-    for ord in 1:corder
+    c = TaylorN(cdivfact, a.order)
+    for ord in 1:a.order
         div!(c, a, b, ord) # updates c[ord+1]
     end
 
@@ -465,7 +447,7 @@ c_k =  \frac{1}{b_0} (a_k - \sum_{j=0}^{k-1} c_j b_{k-j}).
 For `Taylor1` polynomials, `ordfact` is the order of the factorized
 term of the denominator.
 """
-function div!(c::Taylor1, a::Taylor1, b::Taylor1, k::Int, ordfact::Int=0)
+@inline function div!(c::Taylor1, a::Taylor1, b::Taylor1, k::Int, ordfact::Int=0)
     if k == 0
         @inbounds c[1] = a[ordfact+1] / b[ordfact+1]
         return nothing
@@ -478,7 +460,7 @@ function div!(c::Taylor1, a::Taylor1, b::Taylor1, k::Int, ordfact::Int=0)
     return nothing
 end
 
-function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
+@inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int)
     if k==0
         @inbounds c[1] = a[1] / constant_term(b)
         return nothing

src/auxiliary.jl

@@ -18,7 +18,8 @@ function resize_coeffs1!{T<:Number}(coeffs::Array{T,1}, order::Int)
     lencoef = length(coeffs)
     order ≤ lencoef-1 && return nothing
     resize!(coeffs, order+1)
-    coeffs[lencoef+1:end] .= zero(coeffs[1])
+    z = zero(coeffs[1])
+    @__dot__ coeffs[lencoef+1:order+1] = z
     return nothing
 end
 
@@ -35,7 +36,8 @@ function resize_coeffsHP!{T<:Number}(coeffs::Array{T,1}, order::Int)
     @assert order ≤ get_order() && lencoef ≤ num_coeffs
     num_coeffs == lencoef && return nothing
     resize!(coeffs, num_coeffs)
-    coeffs[lencoef+1:end] .= zero(coeffs[1])
+    z = zero(coeffs[1])
+    @__dot__ coeffs[lencoef+1:num_coeffs] = z
     return nothing
 end
 
@@ -71,10 +73,7 @@ Return the coefficient of order `n::Int` of a `a::Taylor1` polynomial.
 get_coeff(a::Taylor1, n::Int) = (@assert 0 ≤ n ≤ a.order+1;
     return a[n+1])
 
-function getindex(a::Taylor1, n::Int)
-    (1 ≤ n ≤ length(a.coeffs)) && return a.coeffs[n]
-    return zero(a.coeffs[1])
-end
+getindex(a::Taylor1, n::Int) = a.coeffs[n]
 getindex(a::Taylor1, n::UnitRange) = view(a.coeffs, n)
 getindex(a::Taylor1, c::Colon) = view(a.coeffs, c)
 setindex!{T<:Number}(a::Taylor1{T}, x::T, n::Int) = a.coeffs[n] = x
@@ -124,11 +123,7 @@ function get_coeff(a::TaylorN, v::Array{Int,1})
     get_coeff(a[order+1], v)
 end
 
-function getindex(a::TaylorN, n::Int)
-    1 ≤ n ≤ length(a.coeffs) && return a.coeffs[n]
-    n ≤ get_order()+1 && return zero_korder(a, n-1)
-    throw(BoundsError(a.coeffs,n))
-end
+getindex(a::TaylorN, n::Int) = a.coeffs[n]
 getindex(a::TaylorN, n::UnitRange) = view(a.coeffs, n)
 getindex(a::TaylorN, c::Colon) = view(a.coeffs, c)
 setindex!{T<:Number}(a::TaylorN{T}, x::HomogeneousPolynomial{T}, n::Int) =
@@ -167,17 +162,22 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
 end
 
 
-"""
-    max_order(a)
-
-Returns the maximum order of a `Taylor1` or `TaylorN`.
-
-For `a::Taylor1` returns `a.order` while for `a::TaylorN` returns
-`get_order()`.
-"""
-max_order(a::Taylor1) = a.order
+## fixorder ##
+for T in (:Taylor1, :TaylorN)
+    @eval begin
+        @inline function fixorder(a::$T, b::$T)
+            a.order == b.order && return a, b
+            a.order < b.order && return $T(copy(a.coeffs), b.order),
+                $T(copy(b.coeffs), b.order)
+            return $T(copy(a.coeffs), a.order), $T(copy(b.coeffs), a.order)
+        end
+    end
+end
 
-max_order(::TaylorN) = get_order()
+function fixorder(a::HomogeneousPolynomial, b::HomogeneousPolynomial)
+    @assert a.order == b.order
+    return a, b
+end
 
 
 """

src/constructors.jl

@@ -31,16 +31,19 @@ immutable Taylor1{T<:Number} <: AbstractSeries{T}
     ## Inner constructor ##
     function (::Type{Taylor1{T}}){T}(coeffs::Array{T,1}, order::Int)
         resize_coeffs1!(coeffs, order)
-        return new{T}(coeffs, length(coeffs)-1)
+        return new{T}(coeffs, order)
     end
 end
 
 ## Outer constructors ##
-Taylor1{T<:Number}(x::Taylor1{T}, order::Int) = Taylor1{T}(x.coeffs, order)
 Taylor1{T<:Number}(x::Taylor1{T}) = x
 Taylor1{T<:Number}(coeffs::Array{T,1}, order::Int) = Taylor1{T}(coeffs, order)
-Taylor1{T<:Number}(coeffs::Array{T,1}) = Taylor1{T}(coeffs, length(coeffs)-1)
-Taylor1{T<:Number}(x::T, order::Int) = Taylor1{T}([x], order)
+Taylor1{T<:Number}(coeffs::Array{T,1}) = Taylor1(coeffs, length(coeffs)-1)
+function Taylor1{T<:Number}(x::T, order::Int)
+    v = zeros(T, order+1)
+    v[1] = x
+    return Taylor1(v, order)
+end
 
 # Shortcut to define Taylor1 independent variables
 doc"""
@@ -85,8 +88,6 @@ immutable HomogeneousPolynomial{T<:Number} <: AbstractSeries{T}
     end
 end
 
-HomogeneousPolynomial{T<:Number}(x::HomogeneousPolynomial{T}, order::Int) =
-    HomogeneousPolynomial{T}(x.coeffs, order)
 HomogeneousPolynomial{T<:Number}(x::HomogeneousPolynomial{T}) = x
 HomogeneousPolynomial{T<:Number}(coeffs::Array{T,1}, order::Int) =
     HomogeneousPolynomial{T}(coeffs, order)
@@ -150,7 +151,6 @@ immutable TaylorN{T<:Number} <: AbstractSeries{T}
     end
 end
 
-TaylorN{T<:Number}(x::TaylorN{T}, order::Int) = TaylorN{T}(x.coeffs, order)
 TaylorN{T<:Number}(x::TaylorN{T}) = x
 TaylorN{T<:Number}(x::Array{HomogeneousPolynomial{T},1}, order::Int) =
     TaylorN{T}(x, order)

src/evaluate.jl

@@ -79,9 +79,11 @@ Substitute `x::Taylor1` as independent variable in a `a::Taylor1` polynomial.
 evaluate{T<:Number,S<:Number}(a::Taylor1{T}, x::Taylor1{S}) =
     evaluate(promote(a,x)...)
 function evaluate{T<:Number}(a::Taylor1{T}, x::Taylor1{T})
-    maxord = max(a.order, x.order)
+    if a.order != x.order
+        a, x = fixorder(a, x)
+    end
     @inbounds suma = a[end]
-    @inbounds for k = maxord:-1:1
+    @inbounds for k = a.order:-1:1
         suma = suma*x + a[k]
     end
     suma