STaylor1: add division, abs, rad2deg, deg2rad; fix real: imag, conj

benchmark/benchmarks.jl

@@ -47,7 +47,10 @@ for i in dims
     S["STaylor1{$i,Float64} + STaylor1{$i,Float64}"] = @benchmarkable f(+, $q, $q2, $n)
     S["Taylor1{$i,Float64} - Taylor1{$i,Float64}"] = @benchmarkable f(-, $x, $x2, $n)
     S["STaylor1{$i,Float64} - STaylor1{$i,Float64}"] = @benchmarkable f(-, $q, $q2, $n)
-
+    S["Taylor1{$i,Float64} * Taylor1{$i,Float64}"] = @benchmarkable f(*, $x, $x2, $n)
+    S["STaylor1{$i,Float64} * STaylor1{$i,Float64}"] = @benchmarkable f(*, $q, $q2, $n)
+    S["Taylor1{$i,Float64} / Taylor1{$i,Float64}"] = @benchmarkable f(/, $x, $x2, $n)
+    S["STaylor1{$i,Float64} / STaylor1{$i,Float64}"] = @benchmarkable f(/, $q, $q2, $n)
 end
 
 S = SUITE["functions"] = BenchmarkGroup()
@@ -56,22 +59,9 @@ for i in dims
     x = Taylor1(r)
     q = STaylor1(r)
     y = rand()
-    S["exp(Taylor1{Float64}), len = $i"] = @benchmarkable f(exp, $x, $n)
-    S["exp(STaylor1{$i,Float64}),"] = @benchmarkable f(exp, $q, $n)
-    S["zero(Taylor1{Float64}), len = $i"] = @benchmarkable f(zero, $x, $n)
-    S["zero(STaylor1{$i,Float64}),"] = @benchmarkable f(zero, $q, $n)
-    S["one(Taylor1{Float64}), len = $i"] = @benchmarkable f(one, $x, $n)
-    S["one(STaylor1{$i,Float64}),"] = @benchmarkable f(one, $q, $n)
-    S["iszero(Taylor1{Float64}), len = $i"] = @benchmarkable f(iszero, $x, $n)
-    S["iszero(STaylor1{$i,Float64}),"] = @benchmarkable f(iszero, $q, $n)
-end
-
-#=
-S = SUITE["STaylor1 Arithmetic"] = BenchmarkGroup()
-for i in (5,10,5)
-    S["STaylor1{$i,Float64} Scalar Addition"]
-    S["STaylor1{$i,Float64} Scalar Multiplication"]
-    S["STaylor1{$i,Float64} Addition"]
-    S["STaylor1{$i,Float64} Multiplication"]
+    for g in (exp, abs, zero, one, real, imag, conj, adjoint, iszero, isnan,
+              isinf, deg2rad, rad2deg)
+        S["$(Symbol(g))(Taylor1{Float64}), len = $i"] = @benchmarkable f($g, $x, $n)
+        S["$(Symbol(g))(STaylor1{$i,Float64}),"] = @benchmarkable f($g, $q, $n)
+    end
 end
-=#

src/arithmetic.jl

@@ -520,6 +520,45 @@ function /(a::STaylor1{N,T}, b::T) where {N, T<:Number}
     STaylor1{N,T}(div_tuple_by_scalar(a.coeffs, b))
 end
 
+@generated function /(a::STaylor1{N,T}, b::STaylor1{N,T}) where {N,T<:Number}
+
+    ex_calc = quote end
+    append!(ex_calc.args, Any[nothing for i in 1:N])
+    syms = Symbol[Symbol("c$i") for i in 1:N]
+
+    # add error
+    ex_line = quote
+                  if iszero(b[0])
+                      throw(ArgumentError("""The 0th order STaylor1 coefficient
+                                             must be non-zero for b, (a/b)(x) is
+                                             not differentiable at x=0)."""))
+                  end
+              end
+    ex_calc.args[1] = ex_line
+
+    # add recursion relation
+    for j = 0:(N-1)
+        ex_line = :(a[$(j)])
+        for k = 1:j
+            sym = syms[j-k+1]
+            ex_line = :($ex_line - $sym*b[$k])
+        end
+        sym = syms[j+1]
+        ex_line = :($sym = ($ex_line)/b[0])
+        ex_calc.args[j+2] = ex_line
+    end
+
+    exout = :(($(syms[1]),))
+    for i = 2:N
+        push!(exout.args, syms[i])
+    end
+    return quote
+               Base.@_inline_meta
+               $ex_calc
+               return STaylor1{N,T}($exout)
+            end
+end
+
 /(a::TaylorN{T}, b::TaylorN{S}) where
     {T<:NumberNotSeriesN, S<:NumberNotSeriesN} = /(promote(a,b)...)
 

src/other_functions.jl

@@ -21,7 +21,7 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN)
     @eval isnan(a::$T) = any( isnan.(a.coeffs) )
 end
 for f in (:real, :imag, :conj)
-    @eval ($f)(a::STaylor1{N,T}) where {N,T<:Number} = STaylor1(($f).(a.coeffs), Val(N))
+    @eval ($f)(a::STaylor1{N,T}) where {N,T<:Number} = STaylor1{N,T}(($f).(a.coeffs))
 end
 
 adjoint(a::STaylor1) = conj(a)
@@ -101,6 +101,18 @@ for T in (:Taylor1, :TaylorN)
     end
 end
 
+function abs(a::STaylor1{N,T}) where {N,T<:Real}
+    if a[0] > zero(T)
+        return a
+    elseif a[0] < zero(T)
+        return -a
+    else
+        throw(ArgumentError(
+        """The 0th order Taylor1 coefficient must be non-zero
+        (abs(x) is not differentiable at x=0)."""))
+    end
+end
+
 function mod2pi(a::TaylorN{Taylor1{T}}) where {T<:Real}
     coeffs = copy(a.coeffs)
     @inbounds coeffs[1] = mod2pi( constant_term(a) )
@@ -285,6 +297,12 @@ for T in (:Taylor1, :TaylorN)
     @eval rad2deg(z::$T{T}) where {T<:Real} = z * (180 / convert(float(T), pi))
 end
 
+deg2rad(z::STaylor1{N, T}) where {N, T<:AbstractFloat} = z * (convert(T, pi) / 180)
+deg2rad(z::STaylor1{N, T}) where {N, T<:Real} = z * (convert(float(T), pi) / 180)
+
+rad2deg(z::STaylor1{N, T}) where {N, T<:AbstractFloat} = z * (180 / convert(T, pi))
+rad2deg(z::STaylor1{N, T}) where {N, T<:Real} = z * (180 / convert(float(T), pi))
+
 # Internal mutating deg2rad!, rad2deg! functions
 for T in (:Taylor1, :TaylorN)
     @eval @inline function deg2rad!(v::$T{T}, a::$T{T}, k::Int) where {T<:AbstractFloat}

test/onevariable.jl

@@ -580,11 +580,34 @@ end
     # check that STaylor1 and Taylor yeild same result
     t1 = STaylor1([1.1, 2.1, 3.1])
     t2 = Taylor1([1.1, 2.1, 3.1])
-    @test test_vs_Taylor1(exp(t1), exp(t2))
+    for f in (exp, abs)
+        @test test_vs_Taylor1(f(t1), f(t2))
+    end
+    t1a = STaylor1([2.1, 2.1, 3.1])
+    t2a = Taylor1([2.1, 2.1, 3.1])
+    @test isapprox((t1/t1a)[0], (t2/t2a)[0], atol=1E-10)
+    @test isapprox((t1/t1a)[1], (t2/t2a)[1], atol=1E-10)
+    @test isapprox((t1/t1a)[2], (t2/t2a)[2], atol=1E-10)
+
+    @test isapprox((t1*t1a)[0], (t2*t2a)[0], atol=1E-10)
+    @test isapprox((t1*t1a)[1], (t2*t2a)[1], atol=1E-10)
+    @test isapprox((t1*t1a)[2], (t2*t2a)[2], atol=1E-10)
 
     a = STaylor1([0.0, 1.2, 2.3, 4.5, 0.0])
     @test findfirst(a) == 1
     @test findlast(a) == 3
+
+    a = STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
+    @test isapprox(deg2rad(a)[0], 0.087266, atol=1E-5)
+    @test isapprox(deg2rad(a)[2], 0.040142, atol=1E-5)
+    @test isapprox(rad2deg(a)[0], 286.4788975, atol=1E-5)
+    @test isapprox(rad2deg(a)[2], 131.7802928, atol=1E-5)
+    @test real(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
+    @test imag(a) == STaylor1([0.0, 0.0, 0.0, 0.0, 0.0])
+    @test adjoint(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
+    @test conj(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
+    @test a == abs(a)
+    @test a == abs(-a)
 end
 
 #=