Skip to content

Commit

Permalink
Throw errors for inverse trigonometric functions, with tests, and cha…
Browse files Browse the repository at this point in the history 
…nges in docs
  • Loading branch information
Luis Benet committed Sep 2, 2016
1 parent 6ffee07 commit 5e29918
Show file tree
Hide file tree
Showing 2 changed files with 58 additions and 37 deletions.
88 changes: 53 additions & 35 deletions src/Taylor1.jl
View file
Expand Up @@ -83,29 +83,34 @@ function convert{T<:Integer, S<:AbstractFloat}(::Type{Taylor1{Rational{T}}},
end
Taylor1(v)
end
convert{T<:Number}(::Type{Taylor1{T}}, b::Array{T,1}) = Taylor1(b, length(b)-1)
convert{T<:Number, S<:Number}(::Type{Taylor1{T}}, b::Array{S,1}) =
Taylor1(convert(Array{T,1},b))
convert{T<:Number, S<:Number}(::Type{Taylor1{T}}, b::S) = Taylor1([convert(T,b)], 0)
convert{T<:Number}(::Type{Taylor1{T}}, b::T) = Taylor1([b], 0)

promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{Taylor1{T}}) = Taylor1{T}
promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{Taylor1{S}}) =
Taylor1{promote_type(T, S)}
promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{Array{T,1}}) = Taylor1{T}
promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{Array{S,1}}) =
Taylor1{promote_type(T, S)}
promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{T}) = Taylor1{T}
promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{S}) =
Taylor1{promote_type(T, S)}

## Auxiliary function ##
function firstnonzero{T<:Number}(a::Taylor1{T})
nonzero::Int = a.order+1
for i in eachindex(a.coeffs)
if a.coeffs[i] != zero(T)
function firstnonzero{T<:Number}(ac::Vector{T})
nonzero::Int = length(ac)
for i in eachindex(ac)
if ac[i] != zero(T)
nonzero = i-1
break
end
end
nonzero
end
firstnonzero{T<:Number}(a::Taylor1{T}) = firstnonzero(a.coeffs)

function fixshape{T<:Number}(a::Taylor1{T}, b::Taylor1{T})
if a.order == b.order
Expand Down Expand Up @@ -703,7 +708,6 @@ Return the Taylor expansion of $\sin(a)$, of order `a.order`, for
For details on making the Taylor expansion, see [`TaylorSeries.sincosHomogCoef`](@ref).
"""

sin(a::Taylor1) = sincos(a)[1]

## Cos ## 
Expand All @@ -713,7 +717,7 @@ doc"""
Return the Taylor expansion of $\cos(a)$, of order `a.order`,
for `a::Taylor1` polynomial
For details on making the Taylor expansion, see [`TaylorSeries.sincosHomogCoef`](@ref)..
For details on making the Taylor expansion, see [`TaylorSeries.sincosHomogCoef`](@ref).
"""
cos(a::Taylor1) = sincos(a)[2]

Expand Down Expand Up @@ -821,7 +825,6 @@ end

### INVERSE TRIGONOMETRIC FUNCTIONS ### 


## Arcsin ##
doc"""
asin(a)
Expand All @@ -835,6 +838,8 @@ function asin(a::Taylor1)
@inbounds aux = asin( a.coeffs[1] )
T = typeof(aux)
ac = convert(Array{T,1}, a.coeffs)
ac[1]^2 == one(T) && throw(ArgumentError("""
Recursion formula diverges due to vanishing `sqrt`."""))
rc = sqrt(one(T) - a^2).coeffs
coeffs = zeros(ac)
@inbounds coeffs[1] = aux
Expand All @@ -846,27 +851,29 @@ end

# Homogeneous coefficients for arcsin
doc"""
asinHomogCoef(kcoef, ac, r_f,coeffs)
asinHomogCoef(kcoef, ac, rc, coeffs)
Compute the `k-th` expansion coefficient of $s = \arcsin(a)$ given by
\begin{equation*}
s_k = \frac{1}{ \sqrt{1 - a_0^2} } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} s_j \big),
s_k = \frac{1}{ \sqrt{r_0} } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} s_j \big),
\end{equation*}
with $a$ a `Taylor1` polynomial and $r = \sqrt{1 - a^2}$.
Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the already calculated expansion coefficients `r_f`
of $r$ (see above), and the previews coefficients `coeffs` of `asinHomogCoef`.
`ac` of `a`, the already calculated expansion coefficients `rc`
of $r$ (see above), and the already calculated expansion coefficients
`coeffs` of `asin(a)`.
"""
function asinHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, r_f::Array{T,1}, coeffs::Array{T,1})
function asinHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, rc::Array{T,1},
coeffs::Array{T,1})
kcoef == 0 && return asin( ac[1] )
coefhomog = zero(T)
@inbounds for i in 1:kcoef-1
coefhomog += (kcoef-i) * r_f[i+1] * coeffs[kcoef-i+1]
coefhomog += (kcoef-i) * rc[i+1] * coeffs[kcoef-i+1]
end
@inbounds coefhomog = (ac[kcoef+1] - coefhomog/kcoef) / r_f[1]
@inbounds coefhomog = (ac[kcoef+1] - coefhomog/kcoef) / rc[1]
coefhomog
end

Expand All @@ -883,39 +890,42 @@ function acos(a::Taylor1)
@inbounds aux = asin( a.coeffs[1] )
T = typeof(aux)
ac = convert(Array{T,1}, a.coeffs)
ac[1]^2 == one(T) && throw(ArgumentError("""
Recursion formula diverges due to vanishing `sqrt`."""))
rc = sqrt(one(T) - a^2).coeffs
coeffs = zeros(ac)
@inbounds coeffs[1] = aux
@inbounds coeffs[1] = aux
@inbounds for k in 1:a.order
coeffs[k+1] = asinHomogCoef(k , ac, rc, coeffs)
end
@inbounds coeffs[1] = -acos( a.coeffs[1] )
@inbounds coeffs[1] = -acos( a.coeffs[1] )
Taylor1( -coeffs, a.order )
end

# Homogeneous coefficients for arccos
doc"""
acosHomogCoef(kcoef, ac, r_f, coeffs)
acosHomogCoef(kcoef, ac, rc, coeffs)
Compute the `k-th` expansion coefficient of $c = \arccos(a)$ given by
\begin{equation*}
c_k = - \frac{1}{ \sqrt{1 - a_0^2} } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big),
c_k = - \frac{1}{ r_0 } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big),
\end{equation*}
with $a$ a `Taylor1` polynomial and $r = \sqrt{1 - a^2}$.
Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the already calculated expansion coefficients `r_f`
of $r$ (see above), and the previews coefficients `coeffs` of `acosHomogCoef`.
`ac` of `a`, the already calculated expansion coefficients `rc`
of $r$ (see above), and the already calculated expansion coefficients
`coeffs` of `acos(a)`.
"""
function acosHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, r_f::Array{T,1}, coeffs::Array{T,1})
function acosHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, rc::Array{T,1}, coeffs::Array{T,1})
kcoef == 0 && return asin( ac[1] )
coefhomog = zero(T)
@inbounds for i in 1:kcoef-1
coefhomog += (kcoef-i) * r_f[i+1] * coeffs[kcoef-i+1]
coefhomog += (kcoef-i) * rc[i+1] * coeffs[kcoef-i+1]
end
@inbounds coefhomog = -(ac[kcoef+1] - coefhomog/kcoef) / r_f[1]
@inbounds coefhomog = -(ac[kcoef+1] - coefhomog/kcoef) / rc[1]
coefhomog
end

Expand All @@ -933,6 +943,8 @@ function atan(a::Taylor1)
@inbounds aux = atan( a.coeffs[1] )
T = typeof(aux)
rc = (one(T) + a^2).coeffs
rc[1] == zero(T) && throw(ArgumentError("""
Recursion formula has a pole."""))
ac = convert(Array{T,1}, a.coeffs)
coeffs = zeros(ac)
@inbounds coeffs[1] = aux
Expand All @@ -944,27 +956,28 @@ end

# Homogeneous coefficients for arctan
doc"""
atanHomogCoef(kcoef, ac, r_f, coeffs)
atanHomogCoef(kcoef, ac, rc, coeffs)
Compute the `k-th` expansion coefficient of $c = \arctan(a)$ given by
\begin{equation*}
t_k = \frac{1}{1 + a_0^2 }(a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} t_j) ,
t_k = \frac{1}{r_0}(a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} t_j) ,
\end{equation*}
with $a$ a `Taylor1` polynomial and $r = 1 + a^2$.
Inputs are the `kcoef`-th coefficient, the vector of the expansion coefficients
`ac` of `a`, the already calculated expansion coefficients `r_f`
of $r$ (see above), and the previews coefficients `coeffs` of `atanHomogCoef`.
`ac` of `a`, the already calculated expansion coefficients `rc`
of $r$ (see above), and the already calculated expansion coefficients
`coeffs` of `asin(a)`.
"""
function atanHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, r_f::Array{T,1}, coeffs::Array{T,1})
function atanHomogCoef{T<:Number}(kcoef::Int, ac::Array{T,1}, rc::Array{T,1}, coeffs::Array{T,1})
kcoef == 0 && return atan( ac[1] )
coefhomog = zero(T)
@inbounds for i in 1:kcoef-1
coefhomog += (kcoef-i) * r_f[i+1] * coeffs[kcoef-i+1]
coefhomog += (kcoef-i) * rc[i+1] * coeffs[kcoef-i+1]
end
@inbounds coefhomog = (ac[kcoef+1] - coefhomog/kcoef) / r_f[1]
@inbounds coefhomog = (ac[kcoef+1] - coefhomog/kcoef) / rc[1]
coefhomog
end

Expand Down Expand Up @@ -996,8 +1009,7 @@ end

## Integrating ##
"""
integrate(a, x)
integrate(a)
integrate(a, [x])
Return the integral of `a::Taylor1`. The constant of integration
(0-th order coefficient) is set to `x`, which is zero if ommitted.
Expand All @@ -1015,12 +1027,18 @@ integrate{T<:Number}(a::Taylor1{T}) = integrate(a, zero(T))

## Evaluating ##
"""
evaluate(a, dx)
evaluate(a)
evaluate(a, [dx])
Evaluate a `Taylor1` polynomial using Horner's rule (hand coded). If `dx` is
ommitted, its value is considered as zero.
"""
function evaluate{T<:Number}(a::Taylor1{T}, dx::T)
@inbounds suma = a.coeffs[end]
@inbounds for k = a.order:-1:1
suma = suma*dx + a.coeffs[k]
end
suma
end
function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, dx::S)
R = promote_type(T,S)
@inbounds suma = convert(R, a.coeffs[end])
Expand Down Expand Up @@ -1061,7 +1079,7 @@ function A_mul_B!{T<:Number}(y::Vector{Taylor1{T}}, a::Union{Matrix{T},SparseMat
# determine the maximal order of b
order = maximum([b1.order for b1 in b])

# Use matrices of coefficients (of proper size) and A_mul_Bt!
# Use matrices of coefficients (of proper size) and A_mul_B!
B = zeros(T, k, order+1)
for i = 1:k
@inbounds ord = b[i].order
Expand Down
7 changes: 5 additions & 2 deletions test/runtests.jl
View file
Expand Up @@ -95,12 +95,15 @@ facts("Tests for Taylor1 expansions") do
@fact evaluate(exp(Taylor1([0,1],17)),1.0) == 1.0*e --> true
@fact evaluate(exp(Taylor1([0,1],1))) == 1.0 --> true
@fact evaluate(exp(t),t^2) == exp(t^2) --> true

@fact sin(asin(tsquare)) == tsquare --> true
@fact tan(atan(tsquare)) == tsquare --> true
@fact atan(tan(tsquare)) == tsquare --> true
@fact asin(t) + acos(t) == pi/2 --> true
@fact asin(t) + acos(t) == pi/2 --> true
@fact derivative(acos(t)) == - 1/sqrt(1-t^2) --> true
@fact_throws ArgumentError asin(ta(1.0))
@fact_throws ArgumentError acos(ta(big(1.0)))
@fact_throws ArgumentError atan(ta(1//1*im))

@fact derivative(5, exp(ta(1.0))) == exp(1.0) --> true
@fact derivative(3, exp(ta(1.0pi))) == exp(1.0pi) --> true
Expand Down

0 comments on commit 5e29918

Please sign in to comment.