diff --git a/.travis.yml b/.travis.yml index b608514e..9fc5106c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,13 +5,16 @@ os: - osx julia: - - 0.5 - 0.6 - nightly notifications: email: false +matrix: + allow_failures: + - julia: nightly + script: - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi - julia --check-bounds=yes -e 'Pkg.clone(pwd()); Pkg.build("TaylorSeries"); Pkg.test("TaylorSeries"; coverage=true)' diff --git a/README.md b/README.md index 02909687..154e11f0 100644 --- a/README.md +++ b/README.md @@ -3,8 +3,8 @@ A [Julia](http://julialang.org) package for Taylor expansions in one or more independent variables. -[![TaylorSeries](http://pkg.julialang.org/badges/TaylorSeries_0.5.svg)](http://pkg.julialang.org/?pkg=TaylorSeries) [![TaylorSeries](http://pkg.julialang.org/badges/TaylorSeries_0.6.svg)](http://pkg.julialang.org/?pkg=TaylorSeries) +[![TaylorSeries](http://pkg.julialang.org/badges/TaylorSeries_0.7.svg)](http://pkg.julialang.org/?pkg=TaylorSeries) [![Coverage Status](https://coveralls.io/repos/JuliaDiff/TaylorSeries.jl/badge.svg?branch=master)](https://coveralls.io/github/JuliaDiff/TaylorSeries.jl?branch=master) [![](https://img.shields.io/badge/docs-stable-blue.svg)](http://www.juliadiff.org/TaylorSeries.jl/stable) diff --git a/REQUIRE b/REQUIRE index 2652857b..9927ee10 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,2 @@ -julia 0.5 -Compat 0.17.0 +julia 0.6 +Compat 0.54.0 diff --git a/docs/REQUIRE b/docs/REQUIRE new file mode 100644 index 00000000..0a6d78af --- /dev/null +++ b/docs/REQUIRE @@ -0,0 +1 @@ +Documenter 0.16.1 0.16.1+ diff --git a/docs/make.jl b/docs/make.jl index fbae6db3..8391d3c5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,7 +16,7 @@ makedocs( deploydocs( repo = "github.com/JuliaDiff/TaylorSeries.jl.git", target = "build", - julia = "0.5", + julia = "0.6", osname = "linux", deps = nothing, make = nothing diff --git a/docs/src/api.md b/docs/src/api.md index 31d9d87d..121a9590 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -6,6 +6,7 @@ CurrentModule = TaylorSeries ``` +## Module ```@docs TaylorSeries ``` @@ -17,20 +18,23 @@ Taylor1 HomogeneousPolynomial TaylorN AbstractSeries -ParamsTaylorN ``` ## Functions and methods ```@docs -Taylor1([::Type{Float64}], [order::Int64=1]) -HomogeneousPolynomial{T<:Number}(::Type{T}, ::Int) -TaylorN{T<:Number}(::Type{T}, nv::Int; [order::Int=get_order()]) +Taylor1(::Type{T}, ::Int64=1) where {T<:Number} +HomogeneousPolynomial(::Type{T}, ::Int) where {T<:Number} +TaylorN(::Type{T}, ::Int; ::Int=get_order()) where {T<:Number} set_variables +get_variables show_params_TaylorN -get_coeff +show_monomials +getcoeff evaluate evaluate! +taylor_expand +update! derivative integrate gradient @@ -39,27 +43,33 @@ jacobian! hessian hessian! inverse +abs +norm +isapprox +isfinite +displayBigO ``` ## Internals ```@docs +ParamsTaylorN +_InternalMutFuncs generate_tables generate_index_vectors in_base make_inverse_dict -order_posTb resize_coeffs1! resize_coeffsHP! -zero_korder constant_term mul! -mul!(c::HomogeneousPolynomial, a::HomogeneousPolynomial, b::HomogeneousPolynomial) +mul!(::HomogeneousPolynomial, ::HomogeneousPolynomial, ::HomogeneousPolynomial) +mul!(::Vector{Taylor1{T}}, ::Union{Matrix{T},SparseMatrixCSC{T}},::Vector{Taylor1{T}}) where {T<:Number} div! pow! square sqr! -sqr!(c::HomogeneousPolynomial, a::HomogeneousPolynomial) +sqr!(::HomogeneousPolynomial{T}, ::HomogeneousPolynomial{T}) where {T<:NumberNotSeriesN} sqrt! exp! log! @@ -70,7 +80,12 @@ acos! atan! sinhcosh! tanh! -A_mul_B! +derivative! +_internalmutfunc_call +_dict_unary_ops +_dict_binary_calls +_dict_unary_calls +_dict_binary_ops ``` ## Index diff --git a/docs/src/examples.md b/docs/src/examples.md index c46aa9a3..85ed4fd5 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -20,7 +20,7 @@ The first example shows that the four-square identity holds: \end{eqnarray} ``` which was originally proved by Euler. The code can also be found in -[this test](../../test/identities_Euler.jl) of the package. +[this test](https://github.com/JuliaDiff/TaylorSeries.jl/blob/master/test/identities_Euler.jl) of the package. First, we reset the maximum degree of the polynomial to 4, since the RHS of the equation has *a priori* terms of fourth order, and define the 8 @@ -112,7 +112,7 @@ We note that the above functions use expansions in `Int128`. This is actually required, since some coefficients are larger than `typemax(Int)`: ```@repl fateman -get_coeff(f2, [1,6,7,20]) # coefficient of x y^6 z^7 w^{20} +getcoeff(f2, (1,6,7,20)) # coefficient of x y^6 z^7 w^{20} ans > typemax(Int) length(f2) sum(TaylorSeries.size_table) @@ -127,7 +127,8 @@ monomials in 4 variables. ### Bechmarks The functions described above have been compared against Mathematica v11.1. -The relevant files used for benchmarking can be found [here](../../perf/). +The relevant files used for benchmarking can be found +[here](https://github.com/JuliaDiff/TaylorSeries.jl/tree/master/perf). Running on a MacPro with Intel-Xeon processors 2.7GHz, we obtain that Mathematica requires on average (5 runs) 3.075957 seconds for the computation, while for `fateman1` and `fateman2` above we obtain 2.811391 and 1.490256, diff --git a/docs/src/userguide.md b/docs/src/userguide.md index d89181a6..e31b0cae 100644 --- a/docs/src/userguide.md +++ b/docs/src/userguide.md @@ -31,22 +31,31 @@ order considered for the expansion (fieldname `order`). The coefficients are arranged in ascending order with respect to the degree of the monomial, so that `coeffs[1]` is the constant term, `coeffs[2]` gives the first order term (`t^1`), -etc. This is a dense representation of the polynomial. +etc. Yet, it is possible to have the natural ordering with respect +to the degree; see below. This is a dense representation of the polynomial. The order of the polynomial can be omitted in the constructor, which is then fixed by the length of the -vector of coefficients; otherwise, it defines the maximum -degree of the polynomial, and may be used to truncate it. +vector of coefficients. If the length of the vector does not correspond with +the `order`, `order` is used, which effectively truncates polynomial to degree `order`. ```@repl userguide Taylor1([1, 2, 3],4) # Polynomial of order 4 with coefficients 1, 2, 3 Taylor1([0.0, 1im]) # Also works with complex numbers -Taylor1(ones(8), 2) # Polynomial of order 2 +Taylor1(ones(8), 2) # Polynomial truncated to order 2 shift_taylor(a) = a + Taylor1(typeof(a),5) ## a + taylor-polynomial of order 5 t = shift_taylor(0.0) # Independent variable `t` ``` Note that the information about the maximum order considered is displayed -using a big-𝒪 notation. +using a big-𝒪 notation. In some cases, it is desirable to not display +the big-𝒪 notation. The function [`displayBigO`](@ref) allows to +control whether it is displayed or not. +```@repl userguide +displayBigO(false) # turn-off displaying big O notation +t +displayBigO(true) # turn it on +t +``` The definition of `shift_taylor(a)` uses the method [`Taylor1([::Type{Float64}], [order::Int64=1])`](@ref), which is a @@ -67,7 +76,7 @@ t*(t^2-4)/(t+2) tI = im*t (1-t)^3.2 (1+t)^t -t^6 # order is 5 +t^6 # t is of order 5 ``` If no valid Taylor expansion can be computed, an error is thrown, for instance @@ -107,6 +116,16 @@ sqrt(t) log(t) ``` +To obtain a specific coefficient, [`getcoeff`](@ref) can be used. Another +alternative is to request the specific degree using the vector notation, +where the index corresponds to the degree of the term. + +```@repl userguide +expon = exp(t) +getcoeff(expon, 0) == expon[0] +rationalize(expon[3]) +``` + Differentiating and integrating is straightforward for polynomial expansions in one variable, using [`derivative`](@ref) and [`integrate`](@ref). These functions return the corresponding [`Taylor1`](@ref) expansions. @@ -123,7 +142,7 @@ integrate(exp(t)) integrate( exp(t), 1.0) integrate( derivative( exp(-t)), 1.0 ) == exp(-t) derivative(1, exp(shift_taylor(1.0))) == exp(1.0) -derivative(5, exp(shift_taylor(1.0))) == exp(1.0) # Fifth derivative of `exp(1+t)` +derivative(5, exp(shift_taylor(1.0))) == exp(1.0) # 5-th derivative of `exp(1+t)` ``` To evaluate a Taylor series at a given point, Horner's rule is used via the @@ -142,6 +161,36 @@ eBig = evaluate( exp(tBig), one(BigFloat) ) e - eBig ``` +Another way to obtain the value of a `Taylor1` polynomial `p` at a given value `x`, is to call `p` as if it were a function, i.e., `p(x)`: + +```@repl userguide +t = Taylor1(15) +p = sin(t) +evaluate(p, pi/2) # value of p at pi/2 using `evaluate` +p(pi/2) # value of p at pi/2 by evaluating p as a function +p(pi/2) == evaluate(p, pi/2) +p(0.0) +p() == p(0.0) # p() is a shortcut to obtain the 0-th order coefficient of `p` +``` + +Note that the syntax `p(x)` is equivalent to `evaluate(p, x)`, whereas `p()` is +equivalent to `evaluate(p)`. For more details about function-like behavior for a +given type in Julia, see the [Function-like objects](https://docs.julialang.org/en/stable/manual/methods/#Function-like-objects-1) +section of the Julia manual. + +Useful shortcuts are [`taylor_expand`](@ref) are [`update!`](@ref). +The former returns +the expansion of a function around a given value `t0`. In turn, `update!` +provides an in-place update of a given Taylor polynomial, that is, it shifts +it further by the provided amount. + +```@repl userguide +p = taylor_expand( x -> sin(x), pi/2, order=16) # 16-th order expansion of sin(t) around pi/2 +update!(p, 0.025) # updates the expansion given by p, by shifting it further by 0.025 +p +``` + + ## Many variables A polynomial in ``N>1`` variables can be represented in (at least) two ways: @@ -155,10 +204,10 @@ was discussed [here](https://groups.google.com/forum/#!msg/julia-users/AkK_UdST The structure [`TaylorN`](@ref) is constructed as a vector of parameterized homogeneous polynomials defined by the type [`HomogeneousPolynomial`](@ref), which in turn is a vector of -coefficients of given order (degree). This implementation imposes that the user -has to specify the (maximum) order considered and the number of independent -variables, which is done using the [`set_variables`](@ref) function. -A vector of the resulting Taylor variables is returned: +coefficients of given order (degree). This implementation imposes the user +to specify the (maximum) order considered and the number of independent +variables at the beginning, which can be conveniently done using +[`set_variables`](@ref). A vector of the resulting Taylor variables is returned: ```@repl userguide x, y = set_variables("x y") @@ -168,19 +217,32 @@ x.coeffs ``` As shown, the resulting objects are of `TaylorN{Float64}` type. -There is an optional `order` keyword argument in [`set_variables`](@ref): +There is an optional `order` keyword argument in [`set_variables`](@ref), +used to specify the maximum order of the `TaylorN` polynomials. Note that +one can specify the variables using a vector of symbols. ```@repl userguide -set_variables("x y", order=10) +set_variables([:x, :y], order=10) ``` -Numbered variables are also available by specifying a single +Similarly, numbered variables are also available by specifying a single variable name and the optional keyword argument `numvars`: ```@repl userguide set_variables("α", numvars=3) ``` +Alternatively to `set_variables`, [`get_variables`](@ref) can be used if one +doesn't want to change internal dictionaries. `get_variables` returns a vector +of `TaylorN` independent variables of a desired `order` +(lesser than `get_order` so the +internals doesn't have to change) with the length and variable names defined +by `set_variables` initially. + +```@repl userguide +get_variables(order=2) #vector of independent variables of order 2 +``` + The function [`show_params_TaylorN`](@ref) displays the current values of the parameters, in an info block. @@ -188,12 +250,6 @@ parameters, in an info block. show_params_TaylorN() ``` - julia> show_params_TaylorN() - INFO: Parameters for `TaylorN` and `HomogeneousPolynomial`: - Maximum order = 10 - Number of variables = 3 - Variable names = UTF8String["α₁","α₂","α₃"] - Internally, changing the parameters (maximum order and number of variables) redefines the hash-tables that translate the index of the coefficients of a [`HomogeneousPolynomial`](@ref) @@ -220,7 +276,7 @@ using [`HomogeneousPolynomial`](@ref) objects directly, which is uncomfortable. ```@repl userguide -set_variables("x", numvars=2); +set_variables(:x, numvars=2); # symbols can be used HomogeneousPolynomial([1,-1]) TaylorN([HomogeneousPolynomial([1,0]), HomogeneousPolynomial([1,2,3])],4) ``` @@ -246,25 +302,40 @@ x, y = set_variables("x y", order=10); exy = exp(x+y) ``` -The function [`get_coeff`](@ref) +The function [`getcoeff`](@ref) gives the normalized coefficient of the polynomial that corresponds to the -monomial specified by a vector `v` containing the powers. For instance, for +monomial specified by the tuple or vector `v` containing the powers. +For instance, for the polynomial `exy` above, the coefficient of the monomial ``x^3 y^5`` is - +obtained using `getcoeff(exy, (3,5))` or `getcoeff(exy, [3,5])`. ```@repl userguide -get_coeff(exy, [3,5]) +getcoeff(exy, (3,5)) rationalize(ans) ``` +Similar to `Taylor1`, vector notation can be used to request specific +coefficients of `HomogeneousPolynomial` or `TaylorN` objects. For `TaylorN` +objects, the index refers to the degree of the `HomogeneousPolynomial`. +In the case of `HomogeneousPolynomial` the index refers to the position +of the hash table. The function [`show_monomials`](@ref) can be used to +obtain the coefficient a specific monomial, given the degree of the +`HomogeneousPolynomial`. + +```@repl userguide +exy[8] # get the 8th order term +show_monomials(8) +exy[8][6] # get the 6th coeff of the 8th order term +``` + Partial differentiation is also implemented for [`TaylorN`](@ref) objects, through the function [`derivative`](@ref), specifying the number -of the variable as the second argument; integration is yet to be implemented. +of the variable, or its symbol, as the second argument. ```@repl userguide p = x^3 + 2x^2 * y - 7x + 2 q = y - x^4 derivative( p, 1 ) # partial derivative with respect to 1st variable -derivative( q, 2 ) +derivative( q, :y ) # partial derivative with respect to :y ``` If we ask for the partial derivative with respect to a non-defined variable, @@ -274,13 +345,64 @@ an error is thrown. derivative( q, 3 ) # error, since we are dealing with 2 variables ``` +To obtain more specific partial derivatives we have two specialized methods +that involve a tuple, which represents the number of derivatives with +respect to each variable (so the tuple's length has to be the +same as the actual number of variables). These methods either return +the `TaylorN` object in question, or the coefficient corresponding to +the specified tuple, normalized by the factorials defined by the tuple. +The latter is in essence the 0-th order coefficient of the former. + +```@repl userguide +derivative(p, (2,1)) # two derivatives on :x and one on :y +derivative((2,1), p) # 0-th order coefficient of the previous expression +derivative(p, (1,1)) # one derivative on :x and one on :y +derivative((1,1), p) # 0-th order coefficient of the previous expression +``` + +Integration with respect to the `r`-th variable for +`HomogeneousPolynomial`s and `TaylorN` objects is obtained +using [`integrate`](@ref). Note that `integrate` for `TaylorN` +objects allows to specify a constant of integration, which must +be independent from the integrated variable. Again, the integration +variable may be specified by its symbol. + +```@repl userguide +integrate( derivative( p, 1 ), 1) # integrate with respect to the first variable +integrate( derivative( p, 1 ), :x, 2) # integration with respect to :x, constant of integration is 2 +integrate( derivative( q, 2 ), :y, -x^4) == q +integrate( derivative( q, 2 ), 2, y) +``` + [`evaluate`](@ref) can also be used for [`TaylorN`](@ref) objects, using it on vectors of numbers (`Real` or `Complex`); the length of the vector must coincide with the -number of independent variables. +number of independent variables. [`evaluate`](@ref) also allows to specify only +one variable and a value. ```@repl userguide -evaluate(exy, [.1,.02]) == e^0.12 +evaluate(exy, [.1,.02]) == exp(0.12) +evaluate(exy, :x, 0.0) == exp(y) # evaluate `exy` for :x -> 0 +``` + +Analogously to `Taylor1`, another way to obtain the value of a `TaylorN` +polynomial `p` at a given point `x`, is to call it as if it were a function: +the syntax `p(x)` for `p::TaylorN` is equivalent to `evaluate(p,x)`, and +`p()` is equivalent to `evaluate(p)`. + +```@repl userguide +exy([.1,.02]) == exp(0.12) +exy(:x, 0.0) +``` + +The functions `taylor_expand` and `update!` work as well for `TaylorN`. + +```@repl userguide +xysq = x^2 + y^2 +update!(xysq, [1.0, -2.0]) # expand around (1,-2) +xysq +update!(xysq, [-1.0, 2.0]) # shift-back +xysq == x^2 + y^2 ``` Functions to compute the gradient, Jacobian and @@ -300,19 +422,20 @@ jacobian([p,q], [2,1]) hessian(r, [1.0,1.0]) ``` -Other specific applications are described in the next [section](examples). +Other specific applications are described in the +[Examples](@ref). ## Mixtures As mentioned above, `Taylor1{T}`, `HomogeneousPolynomial{T}` and `TaylorN{T}` -are parameterized structures structures such that `T<:AbstractSeries`, the latter +are parameterized structures such that `T<:AbstractSeries`, the latter is a subtype of `Number`. Then, we may actually define Taylor expansions in ``N+1`` variables, where one of the variables (the `Taylor1` variable) is somewhat special. ```@repl userguide x, y = set_variables("x y", order=3) -t = Taylor1([zero(x), one(x)], 5) +t1N = Taylor1([zero(x), one(x)], 5) ``` The last line defines a `Taylor1{TaylorN{Float64}}` variable, which is of order @@ -320,7 +443,7 @@ The last line defines a `Taylor1{TaylorN{Float64}}` variable, which is of order such polynomials: ```@repl userguide -cos(2.1+x+t) +cos(2.1+x+t1N) ``` This kind of expansions are of interest when studying the dependence of @@ -329,3 +452,7 @@ the dependence of the solution of a differential equation on the initial conditi around a given solution. In this case, `x` and `y` represent small variations around a given value of the parameters, or around some specific initial condition. Such constructions are exploited in the package [`TaylorIntegration.jl`](https://github.com/PerezHz/TaylorIntegration.jl). + +```@meta +CurrentModule = nothing +``` diff --git a/examples/User guide.ipynb b/examples/User guide.ipynb index a319377f..3151617e 100644 --- a/examples/User guide.ipynb +++ b/examples/User guide.ipynb @@ -1267,7 +1267,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The function `get_coeff(a,v)`\n", + "The function `getcoeff(a,v)`\n", "gives the coefficient of `a::TaylorN` that corresponds to the monomial\n", "specified by the vector of powers `v`:" ] @@ -1313,7 +1313,7 @@ } ], "source": [ - "rationalize(get_coeff(exy, [3,5]))" + "rationalize(getcoeff(exy, [3,5]))" ] }, { @@ -2057,7 +2057,7 @@ } ], "source": [ - "get_coeff(f2,[1,6,7,20]) # coef x^1 y^6 z^7 w^{20}" + "getcoeff(f2,[1,6,7,20]) # coef x^1 y^6 z^7 w^{20}" ] }, { diff --git a/src/TaylorSeries.jl b/src/TaylorSeries.jl index fa792c02..cc077f42 100644 --- a/src/TaylorSeries.jl +++ b/src/TaylorSeries.jl @@ -20,32 +20,41 @@ module TaylorSeries using Compat -if !isdefined(Base, :iszero) - import Compat: iszero - export iszero +if VERSION <= v"0.7.0-DEV.2004" + import Base: norm, gradient, A_mul_B! + const mul! = A_mul_B! + export mul! else - import Base: iszero + using InteractiveUtils: subtypes + using SparseArrays: SparseMatrixCSC + import LinearAlgebra: norm, gradient, mul! + import Base: lastindex + using Markdown end import Base: ==, +, -, *, /, ^ -import Base: zero, one, zeros, ones, isinf, isnan, +import Base: zero, one, zeros, ones, isinf, isnan, iszero, convert, promote_rule, promote, eltype, length, show, real, imag, conj, ctranspose, - rem, mod, mod2pi, abs, gradient, + rem, mod, mod2pi, abs, abs2, sqrt, exp, log, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, - A_mul_B!, power_by_squaring, - getindex, setindex!, endof + power_by_squaring, + getindex, setindex!, endof, + rtoldefault, isfinite, isapprox, rad2deg, deg2rad export Taylor1, TaylorN, HomogeneousPolynomial, AbstractSeries -export get_coeff, derivative, integrate, +export getcoeff, derivative, integrate, evaluate, evaluate!, inverse, - show_params_TaylorN, + show_params_TaylorN, show_monomials, displayBigO, get_order, get_numvars, set_variables, get_variables, - ∇, jacobian, jacobian!, hessian, hessian! + get_variable_names, get_variable_symbols, + ∇, jacobian, jacobian!, hessian, hessian!, + taylor_expand, update!, constant_term, + get_nested_numvars, set_nested_variables include("parameters.jl") include("hash_tables.jl") @@ -58,6 +67,8 @@ include("functions.jl") include("other_functions.jl") include("evaluate.jl") include("calculus.jl") +include("dictmutfunct.jl") include("printing.jl") +include("nested_taylor.jl") end # module diff --git a/src/arithmetic.jl b/src/arithmetic.jl index d0943409..10d36737 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -12,9 +12,9 @@ for T in (:Taylor1, :TaylorN) @eval begin - =={T<:Number,S<:Number}(a::$T{T}, b::$T{S}) = ==(promote(a,b)...) + ==(a::$T{T}, b::$T{S}) where {T<:Number, S<:Number} = ==(promote(a,b)...) - function =={T<:Number}(a::$T{T}, b::$T{T}) + function ==(a::$T{T}, b::$T{T}) where {T<:Number} if a.order != b.order a, b = fixorder(a, b) end @@ -23,8 +23,10 @@ for T in (:Taylor1, :TaylorN) end end -==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) = a.coeffs == b.coeffs - +function ==(a::HomogeneousPolynomial, b::HomogeneousPolynomial) + a.order == b.order && return a.coeffs == b.coeffs + return iszero(a.coeffs) && iszero(b.coeffs) +end for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) @eval iszero(a::$T) = iszero(a.coeffs) @@ -32,34 +34,34 @@ end ## zero and one ## for T in (:Taylor1, :TaylorN), f in (:zero, :one) - @eval ($f)(a::$T) = $T(($f)(a[1]), a.order) + @eval ($f)(a::$T) = $T(($f)(a[0]), a.order) end -function zero{T<:Number}(a::HomogeneousPolynomial{T}) - v = zeros(a.coeffs) +function zero(a::HomogeneousPolynomial{T}) where {T<:Number} + v = zero.(a.coeffs) return HomogeneousPolynomial(v, a.order) end -function zeros{T<:Number}(::HomogeneousPolynomial{T}, order::Int) +function zeros(::HomogeneousPolynomial{T}, order::Int) where {T<:Number} order == 0 && return [HomogeneousPolynomial([zero(T)], 0)] - v = Array{HomogeneousPolynomial{T}}(order+1) + @compat v = Array{HomogeneousPolynomial{T}}(undef, order+1) @simd for ord in eachindex(v) @inbounds v[ord] = HomogeneousPolynomial([zero(T)], ord-1) end return v end -zeros{T<:Number}(::Type{HomogeneousPolynomial{T}}, order::Int) = +zeros(::Type{HomogeneousPolynomial{T}}, order::Int) where {T<:Number} = zeros( HomogeneousPolynomial([zero(T)], 0), order) -function one{T<:Number}(a::HomogeneousPolynomial{T}) - v = ones(a.coeffs) - return HomogeneousPolynomial{T}(v, a.order) +function one(a::HomogeneousPolynomial{T}) where {T<:Number} + v = one.(a.coeffs) + return HomogeneousPolynomial(v, a.order) end -function ones{T<:Number}(a::HomogeneousPolynomial{T}, order::Int) +function ones(a::HomogeneousPolynomial{T}, order::Int) where {T<:Number} order == 0 && return [HomogeneousPolynomial([one(a[1])], 0)] - v = Array{HomogeneousPolynomial{T}}(order+1) + @compat v = Array{HomogeneousPolynomial{T}}(undef, order+1) @simd for ord in eachindex(v) @inbounds num_coeffs = size_table[ord] @inbounds v[ord] = HomogeneousPolynomial(ones(T, num_coeffs), ord-1) @@ -67,7 +69,7 @@ function ones{T<:Number}(a::HomogeneousPolynomial{T}, order::Int) return v end -ones{T<:Number}(::Type{HomogeneousPolynomial{T}}, order::Int) = +ones(::Type{HomogeneousPolynomial{T}}, order::Int) where {T<:Number} = ones( HomogeneousPolynomial([one(T)], 0), order) @@ -77,9 +79,10 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!))) for T in (:Taylor1, :TaylorN) @eval begin - ($f){T<:Number,S<:Number}(a::$T{T}, b::$T{S}) = $f(promote(a,b)...) + ($f)(a::$T{T}, b::$T{S}) where {T<:Number, S<:Number} = + $f(promote(a,b)...) - function $f{T<:Number}(a::$T{T}, b::$T{T}) + function $f(a::$T{T}, b::$T{T}) where {T<:Number} if a.order != b.order a, b = fixorder(a, b) end @@ -94,37 +97,55 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!))) return $T(v, a.order) end - ($f){T<:Number,S<:Number}(a::$T{T}, b::S) = $f(promote(a,b)...) + ($f)(a::$T{T}, b::S) where {T<:Number, S<:Number} = + $f(promote(a,b)...) - function $f{T<:Number}(a::$T{T}, b::T) + function $f(a::$T{T}, b::T) where {T<:Number} coeffs = copy(a.coeffs) - @inbounds coeffs[1] = $f(a[1], b) + @inbounds coeffs[1] = $f(a[0], b) return $T(coeffs, a.order) end - ($f){T<:Number,S<:Number}(b::S, a::$T{T}) = $f(promote(b,a)...) + ($f)(b::S, a::$T{T}) where {T<:Number, S<:Number} = + $f(promote(b,a)...) - function $f{T<:Number}(b::T, a::$T{T}) + function $f(b::T, a::$T{T}) where {T<:Number} coeffs = similar(a.coeffs) @__dot__ coeffs = ($f)(a.coeffs) - @inbounds coeffs[1] = $f(b, a[1]) + @inbounds coeffs[1] = $f(b, a[0]) return $T(coeffs, a.order) end ## add! and subst! ## + function ($fc)(v::$T{T}, a::$T{T}, k::Int) where {T} + @inbounds v[k] = ($f)(a[k]) + return nothing + end + function ($fc)(v::$T, a::NumberNotSeries, k::Int) + @inbounds v[k] = k==0 ? ($f)(zero(v[0]),a) : zero(v[k]) + return nothing + end function ($fc)(v::$T, a::$T, b::$T, k::Int) - @inbounds v[k+1] = ($f)(a[k+1], b[k+1]) + @inbounds v[k] = ($f)(a[k], b[k]) + return nothing + end + function ($fc)(v::$T, a::$T, b::NumberNotSeries, k::Int) + @inbounds v[k] = k==0 ? ($f)(a[0], b) : ($f)(a[k], zero(b)) + return nothing + end + function ($fc)(v::$T, a::NumberNotSeries, b::$T, k::Int) + @inbounds v[k] = k==0 ? ($f)(a, b[0]) : ($f)(zero(a), b[k]) return nothing end end end @eval begin - ($f){T<:NumberNotSeriesN,S<:NumberNotSeriesN}(a::HomogeneousPolynomial{T}, - b::HomogeneousPolynomial{S}) = $f(promote(a,b)...) + ($f)(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{S}) where + {T<:NumberNotSeriesN,S<:NumberNotSeriesN} = $f(promote(a,b)...) - function $f{T<:NumberNotSeriesN}(a::HomogeneousPolynomial{T}, - b::HomogeneousPolynomial{T}) + function $f(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where + T<:NumberNotSeriesN @assert a.order == b.order v = similar(a.coeffs) @__dot__ v = $f(a.coeffs, b.coeffs) @@ -137,41 +158,41 @@ for (f, fc) in ((:+, :(add!)), (:-, :(subst!))) return HomogeneousPolynomial(v, a.order) end - function ($f){T<:NumberNotSeries,S<:NumberNotSeries}( - a::TaylorN{Taylor1{T}}, b::Taylor1{S}) - @inbounds aux = $f(a[1][1], b) + function ($f)(a::TaylorN{Taylor1{T}}, b::Taylor1{S}) where + {T<:NumberNotSeries, S<:NumberNotSeries} + @inbounds aux = $f(a[0][1], b) R = eltype(aux) - coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(a.order+1) + @compat coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1) coeffs .= a.coeffs @inbounds coeffs[1] = aux return TaylorN(coeffs, a.order) end - function ($f){T<:NumberNotSeries,S<:NumberNotSeries}( - b::Taylor1{S}, a::TaylorN{Taylor1{T}}) - @inbounds aux = $f(b, a[1][1]) + function ($f)(b::Taylor1{S}, a::TaylorN{Taylor1{T}}) where + {T<:NumberNotSeries,S<:NumberNotSeries} + @inbounds aux = $f(b, a[0][1]) R = eltype(aux) - coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(a.order+1) + @compat coeffs = Array{HomogeneousPolynomial{Taylor1{R}}}(undef, a.order+1) @__dot__ coeffs = $f(a.coeffs) @inbounds coeffs[1] = aux return TaylorN(coeffs, a.order) end - function ($f){T<:NumberNotSeries,S<:NumberNotSeries}( - a::Taylor1{TaylorN{T}}, b::TaylorN{S}) - @inbounds aux = $f(a[1], b) + function ($f)(a::Taylor1{TaylorN{T}}, b::TaylorN{S}) where + {T<:NumberNotSeries,S<:NumberNotSeries} + @inbounds aux = $f(a[0], b) R = eltype(aux) - coeffs = Array{TaylorN{R}}(a.order+1) + @compat coeffs = Array{TaylorN{R}}(undef, a.order+1) coeffs .= a.coeffs @inbounds coeffs[1] = aux return Taylor1(coeffs, a.order) end - function ($f){T<:NumberNotSeries,S<:NumberNotSeries}( - b::TaylorN{S}, a::Taylor1{TaylorN{T}}) - @inbounds aux = $f(b, a[1]) + function ($f)(b::TaylorN{S}, a::Taylor1{TaylorN{T}}) where + {T<:NumberNotSeries,S<:NumberNotSeries} + @inbounds aux = $f(b, a[0]) R = eltype(aux) - coeffs = Array{TaylorN{R}}(a.order+1) + @compat coeffs = Array{TaylorN{R}}(undef, a.order+1) @__dot__ coeffs = $f(a.coeffs) @inbounds coeffs[1] = aux return Taylor1(coeffs, a.order) @@ -189,62 +210,67 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) *(a::$T, b::Bool) = b * a - function *{T<:NumberNotSeries}(a::T, b::$T) - @inbounds aux = a * b[1] - v = Array{typeof(aux)}(length(b.coeffs)) + function *(a::T, b::$T) where {T<:NumberNotSeries} + @inbounds aux = a * b.coeffs[1] + @compat v = Array{typeof(aux)}(undef, length(b.coeffs)) @__dot__ v = a * b.coeffs return $T(v, b.order) end - *{T<:NumberNotSeries}(b::$T, a::T) = a * b + *(b::$T, a::T) where {T<:NumberNotSeries} = a * b end end for T in (:HomogeneousPolynomial, :TaylorN) @eval begin - function *{T<:NumberNotSeries,S<:NumberNotSeries}(a::Taylor1{T}, b::$T{Taylor1{S}}) - @inbounds aux = a * b[1] + function *(a::Taylor1{T}, b::$T{Taylor1{S}}) where + {T<:NumberNotSeries, S<:NumberNotSeries} + + @inbounds aux = a * b.coeffs[1] R = typeof(aux) - coeffs = Array{R}(length(b.coeffs)) + @compat coeffs = Array{R}(undef, length(b.coeffs)) @__dot__ coeffs = a * b.coeffs return $T(coeffs, b.order) end - *{T<:NumberNotSeries,R<:NumberNotSeries}(b::$T{Taylor1{R}}, a::Taylor1{T}) = a * b + *(b::$T{Taylor1{R}}, a::Taylor1{T}) where + {T<:NumberNotSeries,R<:NumberNotSeries} = a * b - function *{T<:NumberNotSeries,S<:NumberNotSeries}(a::$T{T}, b::Taylor1{$T{S}}) - @inbounds aux = a * b[1] + function *(a::$T{T}, b::Taylor1{$T{S}}) where {T<:NumberNotSeries, S<:NumberNotSeries} + @inbounds aux = a * b[0] R = typeof(aux) - coeffs = Array{R}(length(b.coeffs)) + @compat coeffs = Array{R}(undef, length(b.coeffs)) @__dot__ coeffs = a * b.coeffs return Taylor1(coeffs, b.order) end - *{T<:NumberNotSeries,S<:NumberNotSeries}(b::Taylor1{$T{S}}, a::$T{T}) = a * b + *(b::Taylor1{$T{S}}, a::$T{T}) where + {T<:NumberNotSeries,S<:NumberNotSeries} = a * b end end for (T, W) in ((:Taylor1, :Number), (:TaylorN, :NumberNotSeriesN)) - @eval *{T<:$W, S<:$W}(a::$T{T}, b::$T{S}) = *(promote(a,b)...) + @eval *(a::$T{T}, b::$T{S}) where {T<:$W, S<:$W} = *(promote(a,b)...) - @eval function *{T<:$W}(a::$T{T}, b::$T{T}) + @eval function *(a::$T{T}, b::$T{T}) where {T<:$W} if a.order != b.order a, b = fixorder(a, b) end - c = $T(zero(a[1]), a.order) + c = $T(zero(a[0]), a.order) for ord = 0:c.order - mul!(c, a, b, ord) # updates c[ord+1] + mul!(c, a, b, ord) # updates c[ord] end return c end end -*{T<:NumberNotSeriesN,S<:NumberNotSeriesN}(a::HomogeneousPolynomial{T}, - b::HomogeneousPolynomial{S}) = *(promote(a,b)...) +*(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{S}) where + {T<:NumberNotSeriesN,S<:NumberNotSeriesN} = *(promote(a,b)...) -function *{T<:NumberNotSeriesN}(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) +function *(a::HomogeneousPolynomial{T}, b::HomogeneousPolynomial{T}) where + {T<:NumberNotSeriesN} order = a.order + b.order @@ -258,25 +284,34 @@ end # Internal multiplication functions for T in (:Taylor1, :TaylorN) - @eval @inline function mul!(c::$T, a::$T, b::$T, k::Int) + @eval @inline function mul!(c::$T{T}, a::$T{T}, b::$T{T}, k::Int) where {T} - # c[k+1] = zero( a[k+1] ) + # c[k] = zero( a[k] ) @inbounds for i = 0:k if $T == Taylor1 - c[k+1] += a[i+1] * b[k-i+1] + c[k] += a[i] * b[k-i] else - mul!(c[k+1], a[i+1], b[k-i+1]) + mul!(c[k], a[i], b[k-i]) end end return nothing end + + @eval @inline function mul!(v::$T, a::$T, b::NumberNotSeries, k::Int) + @inbounds v[k] = a[k] * b + return nothing + end + @eval @inline function mul!(v::$T, a::NumberNotSeries, b::$T, k::Int) + @inbounds v[k] = a * b[k] + return nothing + end end -doc""" +@doc doc""" mul!(c, a, b, k::Int) --> nothing -Update the `k`-th expansion coefficient `c[k+1]` of `c = a * b`, +Update the `k`-th expansion coefficient `c[k]` of `c = a * b`, where all `c`, `a`, and `b` are either `Taylor1` or `TaylorN`. The coefficients are given by @@ -310,12 +345,12 @@ Return `c = a*b` with no allocation; all arguments are `HomogeneousPolynomial`. @inbounds for na = 1:num_coeffs_a ca = a[na] - ca == zero(T) && continue + iszero(ca) && continue inda = indTa[na] @inbounds for nb = 1:num_coeffs_b cb = b[nb] - cb == zero(T) && continue + iszero(cb) && continue indb = indTb[nb] pos = posTb[inda + indb] @@ -329,38 +364,38 @@ end ## Division ## -function /{T<:Integer, S<:NumberNotSeries}(a::Taylor1{Rational{T}}, b::S) - R = typeof( a[1] // b) - v = Array{R}(a.order+1) +function /(a::Taylor1{Rational{T}}, b::S) where {T<:Integer, S<:NumberNotSeries} + R = typeof( a[0] // b) + @compat v = Array{R}(undef, a.order+1) @__dot__ v = a.coeffs // b return Taylor1(v, a.order) end for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) - @eval function /{T<:NumberNotSeries,S<:NumberNotSeries}(a::$T{T}, b::S) + @eval function /(a::$T{T}, b::S) where {T<:NumberNotSeries, S<:NumberNotSeries} R = promote_type(T,S) return convert($T{R}, a) * inv(convert(R, b)) end - @eval /{T<:NumberNotSeries}(a::$T, b::T) = a * inv(b) + @eval /(a::$T, b::T) where {T<:NumberNotSeries} = a * inv(b) end for T in (:HomogeneousPolynomial, :TaylorN) - @eval function /{T<:NumberNotSeries,S<:NumberNotSeries}( - b::$T{Taylor1{S}}, a::Taylor1{T}) - @inbounds aux = b[1] / a + @eval function /(b::$T{Taylor1{S}}, a::Taylor1{T}) where + {T<:NumberNotSeries,S<:NumberNotSeries} + @inbounds aux = b.coeffs[1] / a R = typeof(aux) - coeffs = Array{R}(length(b.coeffs)) + @compat coeffs = Array{R}(undef, length(b.coeffs)) @__dot__ coeffs = b.coeffs / a return $T(coeffs, b.order) end - @eval function /{T<:NumberNotSeries,S<:NumberNotSeries}( - b::Taylor1{$T{S}}, a::$T{T}) - @inbounds aux = b[1] / a + @eval function /(b::Taylor1{$T{S}}, a::$T{T}) where + {T<:NumberNotSeries,S<:NumberNotSeries} + @inbounds aux = b[0] / a R = typeof(aux) - coeffs = Array{R}(length(b.coeffs)) + @compat coeffs = Array{R}(undef, length(b.coeffs)) @__dot__ coeffs = b.coeffs / a return Taylor1(coeffs, b.order) end @@ -368,9 +403,9 @@ end -/{T<:Number,S<:Number}(a::Taylor1{T}, b::Taylor1{S}) = /(promote(a,b)...) +/(a::Taylor1{T}, b::Taylor1{S}) where {T<:Number, S<:Number} = /(promote(a,b)...) -function /{T<:Number}(a::Taylor1{T}, b::Taylor1{T}) +function /(a::Taylor1{T}, b::Taylor1{T}) where {T<:Number} if a.order != b.order a, b = fixorder(a, b) end @@ -379,17 +414,17 @@ function /{T<:Number}(a::Taylor1{T}, b::Taylor1{T}) ordfact, cdivfact = divfactorization(a, b) c = Taylor1(cdivfact, a.order) - for ord = 1:a.order-ordfact - div!(c, a, b, ord, ordfact) # updates c[ord+1] + for ord = 1:a.order + div!(c, a, b, ord, ordfact) # updates c[ord] end return c end -/{T<:NumberNotSeriesN,S<:NumberNotSeriesN}(a::TaylorN{T}, b::TaylorN{S}) = - /(promote(a,b)...) +/(a::TaylorN{T}, b::TaylorN{S}) where + {T<:NumberNotSeriesN, S<:NumberNotSeriesN} = /(promote(a,b)...) -function /{T<:NumberNotSeriesN}(a::TaylorN{T}, b::TaylorN{T}) +function /(a::TaylorN{T}, b::TaylorN{T}) where {T<:NumberNotSeriesN} @assert !iszero(constant_term(b)) if a.order != b.order @@ -397,10 +432,10 @@ function /{T<:NumberNotSeriesN}(a::TaylorN{T}, b::TaylorN{T}) end # first coefficient - @inbounds cdivfact = a[1] / constant_term(b) + @inbounds cdivfact = a[0] / constant_term(b) c = TaylorN(cdivfact, a.order) for ord in 1:a.order - div!(c, a, b, ord) # updates c[ord+1] + div!(c, a, b, ord) # updates c[ord] end return c @@ -414,12 +449,12 @@ function divfactorization(a1::Taylor1, b1::Taylor1) a1nz = a1nz ≥ 0 ? a1nz : a1.order b1nz = b1nz ≥ 0 ? b1nz : a1.order ordfact = min(a1nz, b1nz) - cdivfact = a1[ordfact+1] / b1[ordfact+1] + cdivfact = a1[ordfact] / b1[ordfact] # Is the polynomial factorizable? - iszero(b1[ordfact+1]) && throw( ArgumentError( + iszero(b1[ordfact]) && throw( ArgumentError( """Division does not define a Taylor1 polynomial; - order k=$(ordfact) => coeff[$(ordfact+1)]=$(cdivfact).""") ) + order k=$(ordfact) => coeff[$(ordfact)]=$(cdivfact).""") ) return ordfact, cdivfact end @@ -429,10 +464,10 @@ end # Homogeneous coefficient for the division -doc""" +@doc doc""" div!(c, a, b, k::Int, ordfact::Int=0) -Compute the `k-th` expansion coefficient `c[k+1]` of `c = a / b`, +Compute the `k-th` expansion coefficient `c[k]` of `c = a / b`, where all `c`, `a` and `b` are either `Taylor1` or `TaylorN`. The coefficients are given by @@ -443,30 +478,44 @@ c_k = \frac{1}{b_0} \big(a_k - \sum_{j=0}^{k-1} c_j b_{k-j}\big). For `Taylor1` polynomials, `ordfact` is the order of the factorized term of the denominator. -""" +""" div! + @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] + @inbounds c[0] = a[ordfact] / b[ordfact] return nothing end @inbounds for i = 0:k-1 - c[k+1] += c[i+1] * b[k+ordfact-i+1] + k+ordfact-i > b.order && continue + c[k] += c[i] * b[k+ordfact-i] + end + if k+ordfact ≤ b.order + @inbounds c[k] = (a[k+ordfact]-c[k]) / b[ordfact] + else + @inbounds c[k] = - c[k] / b[ordfact] end - @inbounds c[k+1] = (a[k+ordfact+1]-c[k+1]) / b[ordfact+1] return nothing end +@inline function div!(v::Taylor1, a::Taylor1, b::NumberNotSeries, k::Int) + @inbounds v[k] = a[k] / b + return nothing +end + +div!(v::Taylor1, b::NumberNotSeries, a::Taylor1, k::Int) = + div!(v::Taylor1, Taylor1(b, a.order), a, k) + @inline function div!(c::TaylorN, a::TaylorN, b::TaylorN, k::Int) if k==0 - @inbounds c[1] = a[1] / constant_term(b) + @inbounds c[0] = a[0] / constant_term(b) return nothing end @inbounds for i = 0:k-1 - mul!(c[k+1], c[i+1], b[k-i+1]) + mul!(c[k], c[i], b[k-i]) end - @inbounds c[k+1] = (a[k+1] - c[k+1]) / constant_term(b) + @inbounds c[k] = (a[k] - c[k]) / constant_term(b) return nothing end @@ -475,13 +524,13 @@ end """ - A_mul_B!(Y, A, B) + mul!(Y, A, B) Multiply A*B and save the result in Y. """ -function A_mul_B!{T<:Number}(y::Vector{Taylor1{T}}, +function mul!(y::Vector{Taylor1{T}}, a::Union{Matrix{T},SparseMatrixCSC{T}}, - b::Vector{Taylor1{T}}) + b::Vector{Taylor1{T}}) where {T<:Number} n, k = size(a) @assert (length(y)== n && length(b)== k) @@ -489,16 +538,18 @@ function A_mul_B!{T<:Number}(y::Vector{Taylor1{T}}, # determine the maximal order of b order = maximum([b1.order for b1 in b]) - # Use matrices of coefficients (of proper size) and A_mul_B! - B = zeros(T, k, order+1) + # Use matrices of coefficients (of proper size) and mul! + # B = zeros(T, k, order+1) + @compat B = Array{T}(undef, k, order+1) + B = zero.(B) for i = 1:k @inbounds ord = b[i].order @inbounds for j = 1:ord+1 - B[i,j] = b[i][j] + B[i,j] = b[i][j-1] end end - Y = Array{T}(n, order+1) - A_mul_B!(Y, a, B) + @compat Y = Array{T}(undef, n, order+1) + @compat mul!(Y, a, B) @inbounds for i = 1:n y[i] = Taylor1( collect(Y[i,:]), order) end diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 781a0fe1..7d6cbb68 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -14,7 +14,7 @@ If the length of `coeffs` is smaller than `order+1`, it resizes `coeffs` appropriately filling it with zeros. """ -function resize_coeffs1!{T<:Number}(coeffs::Array{T,1}, order::Int) +function resize_coeffs1!(coeffs::Array{T,1}, order::Int) where {T<:Number} lencoef = length(coeffs) resize!(coeffs, order+1) if order > lencoef-1 @@ -31,7 +31,7 @@ If the length of `coeffs` is smaller than the number of coefficients correspondinf to `order` (given by `size_table[order+1]`), it resizes `coeffs` appropriately filling it with zeros. """ -function resize_coeffsHP!{T<:Number}(coeffs::Array{T,1}, order::Int) +function resize_coeffsHP!(coeffs::Array{T,1}, order::Int) where {T<:Number} lencoef = length( coeffs ) @inbounds num_coeffs = size_table[order+1] @assert order ≤ get_order() && lencoef ≤ num_coeffs @@ -43,7 +43,7 @@ function resize_coeffsHP!{T<:Number}(coeffs::Array{T,1}, order::Int) end ## Minimum order of an HomogeneousPolynomial compatible with the vector's length -function orderH{T}(coeffs::Array{T,1}) +function orderH(coeffs::Array{T,1}) where {T<:Number} ord = 0 ll = length(coeffs) for i = 1:get_order()+1 @@ -55,7 +55,7 @@ function orderH{T}(coeffs::Array{T,1}) end ## Maximum order of a HomogeneousPolynomial vector; used by TaylorN constructor -function maxorderH{T<:Number}(v::Array{HomogeneousPolynomial{T},1}) +function maxorderH(v::Array{HomogeneousPolynomial{T},1}) where {T<:Number} m = 0 @inbounds for i in eachindex(v) m = max(m, v[i].order) @@ -65,111 +65,151 @@ end -## get_coeff ## +## getcoeff ## """ - get_coeff(a, n) + getcoeff(a, n) 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]) +getcoeff(a::Taylor1, n::Int) = (@assert 0 ≤ n ≤ a.order; return a[n]) -getindex(a::Taylor1, n::Int) = a.coeffs[n] -getindex(a::Taylor1, n::UnitRange) = view(a.coeffs, n) +getindex(a::Taylor1, n::Int) = a.coeffs[n+1] +getindex(a::Taylor1, u::UnitRange) = view(a.coeffs, u .+ 1 ) getindex(a::Taylor1, c::Colon) = view(a.coeffs, c) -setindex!{T<:Number}(a::Taylor1{T}, x::T, n::Int) = a.coeffs[n] = x -setindex!{T<:Number}(a::Taylor1{T}, x::T, n::UnitRange) = a.coeffs[n] = x -setindex!{T<:Number}(a::Taylor1{T}, x::Array{T,1}, n::UnitRange) = a.coeffs[n] = x -setindex!{T<:Number}(a::Taylor1{T}, x::T, c::Colon) = a.coeffs[c] = x -setindex!{T<:Number}(a::Taylor1{T}, x::Array{T,1}, c::Colon) = a.coeffs[c] = x + +setindex!(a::Taylor1{T}, x::T, n::Int) where {T<:Number} = a.coeffs[n+1] = x +setindex!(a::Taylor1{T}, x::T, u::UnitRange) where {T<:Number} = + a.coeffs[u .+ 1] = x +function setindex!(a::Taylor1{T}, x::Array{T,1}, u::UnitRange) where {T<:Number} + # a.coeffs[u .+ 1] .= x + @assert length(u) == length(x) + for ind in eachindex(x) + a.coeffs[u[ind]+1] = x[ind] + end +end +setindex!(a::Taylor1{T}, x::T, c::Colon) where {T<:Number} = a.coeffs[c] = x +setindex!(a::Taylor1{T}, x::Array{T,1}, c::Colon) where {T<:Number} = a.coeffs[c] = x """ - get_coeff(a, v) + getcoeff(a, v) -Return the coefficient of `a::HomogeneousPolynomial`, specified by -`v::Array{Int,1}` which has the indices of the specific monomial. +Return the coefficient of `a::HomogeneousPolynomial`, specified by `v`, +which is a tuple (or vector) with the indices of the specific +monomial. """ -function get_coeff(a::HomogeneousPolynomial, v::Array{Int,1}) - @assert length(v) == get_numvars() +function getcoeff(a::HomogeneousPolynomial, v::NTuple{N,Int}) where {N} + @assert N == get_numvars() && all(v .>= 0) kdic = in_base(get_order(),v) @inbounds n = pos_table[a.order+1][kdic] a[n] end +getcoeff(a::HomogeneousPolynomial, v::Array{Int,1}) = getcoeff(a, (v...,)) getindex(a::HomogeneousPolynomial, n::Int) = a.coeffs[n] getindex(a::HomogeneousPolynomial, n::UnitRange) = view(a.coeffs, n) getindex(a::HomogeneousPolynomial, c::Colon) = view(a.coeffs, c) -setindex!{T<:Number}(a::HomogeneousPolynomial{T}, x::T, n::Int) = + +setindex!(a::HomogeneousPolynomial{T}, x::T, n::Int) where {T<:Number} = a.coeffs[n] = x -setindex!{T<:Number}(a::HomogeneousPolynomial{T}, x::T, n::UnitRange) = +setindex!(a::HomogeneousPolynomial{T}, x::T, n::UnitRange) where {T<:Number} = a.coeffs[n] = x -setindex!{T<:Number}(a::HomogeneousPolynomial{T}, x::Array{T,1}, n::UnitRange) = +setindex!(a::HomogeneousPolynomial{T}, x::Array{T,1}, n::UnitRange) where {T<:Number} = a.coeffs[n] = x -setindex!{T<:Number}(a::HomogeneousPolynomial{T}, x::T, c::Colon) = +setindex!(a::HomogeneousPolynomial{T}, x::T, c::Colon) where {T<:Number} = a.coeffs[c] = x -setindex!{T<:Number}(a::HomogeneousPolynomial{T}, x::Array{T,1}, c::Colon) = +setindex!(a::HomogeneousPolynomial{T}, x::Array{T,1}, c::Colon) where {T<:Number} = a.coeffs[c] = x """ - get_coeff(a, v) + getcoeff(a, v) -Return the coefficient of `a::TaylorN`, specified by -`v::Array{Int,1}` which has the indices of the specific monomial. +Return the coefficient of `a::TaylorN`, specified by `v`, +which is a tuple (or vector) with the indices of the specific +monomial. """ -function get_coeff(a::TaylorN, v::Array{Int,1}) +function getcoeff(a::TaylorN, v::NTuple{N, Int}) where {N} order = sum(v) @assert order ≤ a.order - get_coeff(a[order+1], v) + getcoeff(a[order], v) end +getcoeff(a::TaylorN, v::Array{Int,1}) = getcoeff(a, (v...,)) -getindex(a::TaylorN, n::Int) = a.coeffs[n] -getindex(a::TaylorN, n::UnitRange) = view(a.coeffs, n) +getindex(a::TaylorN, n::Int) = a.coeffs[n+1] +getindex(a::TaylorN, u::UnitRange) = view(a.coeffs, u .+ 1) getindex(a::TaylorN, c::Colon) = view(a.coeffs, c) -setindex!{T<:Number}(a::TaylorN{T}, x::HomogeneousPolynomial{T}, n::Int) = - a.coeffs[n] = x -setindex!{T<:Number}(a::TaylorN{T}, x::T, n::Int) = - a.coeffs[n] = HomogeneousPolynomial(x, n-1) -setindex!{T<:Number}(a::TaylorN{T}, x::HomogeneousPolynomial{T}, n::UnitRange) = - a.coeffs[n] = x -setindex!{T<:Number}(a::TaylorN{T}, x::T, n::UnitRange) = - a.coeffs[n] = x -setindex!{T<:Number}(a::TaylorN{T}, x::Array{HomogeneousPolynomial{T},1}, n::UnitRange) = - a.coeffs[n] = x -setindex!{T<:Number}(a::TaylorN{T}, x::Array{T,1}, n::UnitRange) = - a.coeffs[n] = x -setindex!{T<:Number}(a::TaylorN{T}, x::HomogeneousPolynomial{T}, c::Colon) = - a.coeffs[c] = x -setindex!{T<:Number}(a::TaylorN{T}, x::T, c::Colon) = - a.coeffs[c] = x -setindex!{T<:Number}(a::TaylorN{T}, x::Array{HomogeneousPolynomial{T},1}, c::Colon) = - a.coeffs[c] = x -setindex!{T<:Number}(a::TaylorN{T}, x::Array{T,1}, c::Colon) = - a.coeffs[c] = x + +function setindex!(a::TaylorN{T}, x::HomogeneousPolynomial{T}, n::Int) where + {T<:Number} + @assert x.order == n + a.coeffs[n+1] = x +end +setindex!(a::TaylorN{T}, x::T, n::Int) where {T<:Number} = + a.coeffs[n+1] = HomogeneousPolynomial(x, n) +function setindex!(a::TaylorN{T}, x::T, u::UnitRange) where {T<:Number} + for ind in u + a[ind] = x + end + a[u] +end +function setindex!(a::TaylorN{T}, x::Array{HomogeneousPolynomial{T},1}, u::UnitRange) where {T<:Number} + # a[u[:]] .= x[:] + @assert length(u) == length(x) + for ind in eachindex(x) + a[u[ind]] = x[ind] + end +end +function setindex!(a::TaylorN{T}, x::Array{T,1}, u::UnitRange) where {T<:Number} + # a[u] .= x + @assert length(u) == length(x) + for ind in eachindex(x) + a[u[ind]] = x[ind] + end +end +setindex!(a::TaylorN{T}, x::T, ::Colon) where {T<:Number} = + (a[0:end] = x; a[:]) +setindex!(a::TaylorN{T}, x::Array{HomogeneousPolynomial{T},1}, ::Colon) where + {T<:Number} = (a[0:end] = x; a[:]) +setindex!(a::TaylorN{T}, x::Array{T,1}, ::Colon) where {T<:Number} = + (a[0:end] = x; a[:]) ## eltype, length, endof, get_order ## -for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) +for T in (:Taylor1, :TaylorN) @eval begin - eltype{T<:Number}(::$T{T}) = T - + eltype(::$T{S}) where {S<:Number} = S length(a::$T) = length(a.coeffs) - - endof(a::$T) = length(a.coeffs) - + endof(a::$T) = a.order + @compat lastindex(a::$T) = a.order get_order(a::$T) = a.order + + # Use `a[i0:i1]` or `a[:]` for iterations; see discussion in #140 + # start(a::$T) = start(a.coeffs)-1 + # next(a::$T, ord) = ($T(a[ord], ord), ord+1) + # done(a::$T, ord) = ord == a.order+1 end end +eltype(::HomogeneousPolynomial{S}) where {S<:Number} = S +length(a::HomogeneousPolynomial) = length(a.coeffs) +endof(a::HomogeneousPolynomial) = length(a.coeffs) +@compat lastindex(a::HomogeneousPolynomial) = length(a.coeffs) +get_order(a::HomogeneousPolynomial) = a.order + +# Use `a[i0:i1]` or `a[:]` for iterations; see discussion in #140 +# start(a::HomogeneousPolynomial) = start(a.coeffs) +# next(a::HomogeneousPolynomial, ord) = (HomogeneousPolynomial(a[ord]), ord+1) +# done(a::HomogeneousPolynomial, ord) = ord > length(a.coeffs) + ## 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) + 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 @@ -181,39 +221,17 @@ function fixorder(a::HomogeneousPolynomial, b::HomogeneousPolynomial) end -""" - zero_korder(a) - -For `a::Taylor1` returns `zero(a[1])` while for `a::TaylorN` returns -a zero of a k-th order `HomogeneousPolynomial` of proper type. -""" -zero_korder(a::Taylor1, ::Int) = zero(a[1]) - -zero_korder(a::TaylorN, k::Int) = HomogeneousPolynomial(zero(a[1][1]), k) - - # Finds the first non zero entry; extended to Taylor1 -Base.findfirst{T<:Number}(a::Taylor1{T}) = findfirst(a.coeffs)-1 - - -""" - order_posTb(order, nv, ord) - -Return a vector with the positions, in a `HomogeneousPolynomial` of -order `order`, where the variable `nv` has order `ord`. -""" -function order_posTb(order::Int, nv::Int, ord::Int) - @assert order ≤ get_order() - @inbounds indTb = coeff_table[order+1] - @inbounds num_coeffs = size_table[order+1] - posV = Int[] - for pos = 1:num_coeffs - @inbounds indTb[pos][nv] != ord && continue - push!(posV, pos) - end - posV +function Base.findfirst(a::Taylor1{T}) where {T<:Number} + first = findfirst(a.coeffs) + @compat isa(first, Nothing) && return -1 + return first-1 +end +function Base.findlast(a::Taylor1{T}) where {T<:Number} + last = findlast(a.coeffs) + @compat isa(last, Nothing) && return -1 + return last-1 end - """ constant_term(a) @@ -221,6 +239,10 @@ end Return the constant value (zero order coefficient) for `Taylor1` and `TaylorN`. """ -constant_term(a::Taylor1) = a[1] +constant_term(a::Taylor1) = a[0] + +constant_term(a::TaylorN) = a[0][1] + +constant_term(a::Vector{T}) where {T<:Number}= a[1] -constant_term(a::TaylorN) = a[1][1] +constant_term(a::Number) = a diff --git a/src/calculus.jl b/src/calculus.jl index 1e9699d2..29bc61b9 100644 --- a/src/calculus.jl +++ b/src/calculus.jl @@ -14,12 +14,60 @@ Return the `Taylor1` polynomial of the differential of `a::Taylor1`. The last coefficient is set to zero. """ function derivative(a::Taylor1) - coeffs = zero(a.coeffs) - @inbounds coeffs[1] = a[2] - @inbounds for i = 1:a.order - coeffs[i] = i*a[i+1] + res = zero(a) + @inbounds for ord = 0:a.order-1 + derivative!(res, a, ord) + end + return res +end + +""" + derivative!(res, a) --> nothing + +In-place version of `derivative`. Compute the `Taylor1` polynomial of the +differential of `a::Taylor1` and save it into `res`. The last coefficient is +set to zero. +""" +function derivative!(res::Taylor1, a::Taylor1) + @inbounds for ord = 0:a.order-1 + derivative!(res, a, ord) + end + res[a.order] = zero(a[0]) + nothing +end + +""" + derivative!(p, a, k) --> nothing + +Update in-place the `k-th` expansion coefficient `p[k]` of `p = derivative(a)` +for both `p` and `a` `Taylor1`. + +The coefficients are given by + +```math +p_k = (k+1)a_{k+1}. +``` + +""" +derivative!(p::Taylor1, a::Taylor1, k::Int) = (p[k] = (k+1)*a[k+1]) + +""" + derivative(a, n) + +Compute recursively the `Taylor1` polynomial of the n-th derivative of +`a::Taylor1`. +""" +function derivative(a::Taylor1{T}, n::Int) where {T <: Number} + @assert a.order ≥ n ≥ 0 + if n==0 + return a + else + res = deepcopy(a) + for i in 1:n + derivative!(res, res) + end + return res end - return Taylor1(coeffs, a.order) end """ @@ -27,9 +75,9 @@ end Return the value of the `n`-th derivative of the polynomial `a`. """ -function derivative{T<:Number}(n::Int, a::Taylor1{T}) +function derivative(n::Int, a::Taylor1{T}) where {T<:Number} @assert a.order ≥ n ≥ 0 - factorial( widen(n) ) * a[n+1] :: T + factorial( widen(n) ) * a[n] :: T end ## Integrating ## @@ -39,16 +87,16 @@ end Return the integral of `a::Taylor1`. The constant of integration (0-th order coefficient) is set to `x`, which is zero if ommitted. """ -function integrate{T<:Number, S<:Number}(a::Taylor1{T}, x::S) - R = promote_type(T, typeof(a[1] / 1), S) +function integrate(a::Taylor1{T}, x::S) where {T<:Number, S<:Number} + R = promote_type(T, typeof(a[0] / 1), S) coeffs = zeros(R, a.order+1) @inbounds for i = 1:a.order - coeffs[i+1] = a[i] / i + coeffs[i+1] = a[i-1] / i end @inbounds coeffs[1] = convert(R, x) return Taylor1(coeffs, a.order) end -integrate{T<:Number}(a::Taylor1{T}) = integrate(a, zero(T)) +integrate(a::Taylor1{T}) where {T<:Number} = integrate(a, zero(T)) @@ -81,23 +129,69 @@ function derivative(a::HomogeneousPolynomial, r::Int) return HomogeneousPolynomial{T}(coeffs, a.order-1) end +derivative(a::HomogeneousPolynomial, s::Symbol) = derivative(a, lookupvar(s)) """ - derivative(a, [r=1]) + derivative(a, r) Partial differentiation of `a::TaylorN` series with respect -to the `r`-th variable. +to the `r`-th variable. The `r`-th variable may be also +specified through its symbol. """ function derivative(a::TaylorN, r=1::Int) T = eltype(a) - coeffs = Array{HomogeneousPolynomial{T}}(a.order) + @compat coeffs = Array{HomogeneousPolynomial{T}}(undef, a.order) - @inbounds for ord in eachindex(coeffs) - ord == a.order+1 && continue - coeffs[ord] = derivative( a[ord+1], r) + @inbounds for ord = 1:a.order + coeffs[ord] = derivative( a[ord], r) end return TaylorN{T}( coeffs, a.order ) end +derivative(a::TaylorN, s::Symbol) = derivative(a, lookupvar(s)) + +""" + derivative(a::TaylorN{T}, ntup::NTuple{N,Int}) + +Return a `TaylorN` with the partial derivative of `a` defined +by `ntup::NTuple{N,Int}`, where the first entry is the number +of derivatives with respect to the first variable, the second is +the number of derivatives with respect to the second, and so on. +""" +function derivative(a::TaylorN, ntup::NTuple{N,Int}) where {N} + + @assert N == get_numvars() && all(ntup .>= 0) + + sum(ntup) > a.order && return zero(a) + sum(ntup) == 0 && return copy(a) + + aa = copy(a) + for nvar in 1:get_numvars() + for numder in 1:ntup[nvar] + aa = derivative(aa, nvar) + end + end + + return aa +end + +""" + derivative(ntup::NTuple{N,Int}, a::TaylorN{T}) + +Returns the value of the coefficient of `a` specified by +`ntup::NTuple{N,Int}`, multiplied by the corresponding +factorials. +""" +function derivative(ntup::NTuple{N,Int}, a::TaylorN) where {N} + + @assert N == get_numvars() && all(ntup .>= 0) + + c = getcoeff(a, [ntup...]) + for ind = 1:get_numvars() + c *= factorial(ntup[ind]) + end + + return c +end ## Gradient, jacobian and hessian @@ -112,7 +206,7 @@ Compute the gradient of the polynomial `f::TaylorN`. function gradient(f::TaylorN) T = eltype(f) numVars = get_numvars() - grad = Array{TaylorN{T}}(numVars) + @compat grad = Array{TaylorN{T}}(undef, numVars) @inbounds for nv = 1:numVars grad[nv] = derivative(f, nv) end @@ -129,22 +223,24 @@ const ∇ = gradient Compute the jacobian matrix of `vf`, a vector of `TaylorN` polynomials, evaluated at the vector `vals`. If `vals` is ommited, it is evaluated at zero. """ -function jacobian{T<:Number}(vf::Array{TaylorN{T},1}) +function jacobian(vf::Array{TaylorN{T},1}) where {T<:Number} numVars = get_numvars() @assert length(vf) == numVars - jac = Array{T}(numVars,numVars) + @compat jac = Array{T}(undef, numVars, numVars) @inbounds for comp = 1:numVars - jac[:,comp] = vf[comp][2][1:end] + jac[:,comp] = vf[comp][1][1:end] end return transpose(jac) end -function jacobian{T<:Number,S<:Number}(vf::Array{TaylorN{T},1},vals::Array{S,1}) +function jacobian(vf::Array{TaylorN{T},1}, vals::Array{S,1}) where + {T<:Number, S<:Number} + R = promote_type(T,S) numVars = get_numvars() @assert length(vf) == numVars == length(vals) - jac = Array{R}(numVars,numVars) + @compat jac = Array{R}(undef, numVars, numVars) for comp = 1:numVars @inbounds grad = gradient( vf[comp] ) @@ -156,7 +252,7 @@ function jacobian{T<:Number,S<:Number}(vf::Array{TaylorN{T},1},vals::Array{S,1}) return transpose(jac) end -function jacobian{T<:Number}(vf::Array{Taylor1{TaylorN{T}},1}) +function jacobian(vf::Array{Taylor1{TaylorN{T}},1}) where {T<:Number} vv = convert(Array{TaylorN{Taylor1{T}},1}, vf) jacobian(vv) end @@ -168,20 +264,23 @@ end ``` Compute the jacobian matrix of `vf`, a vector of `TaylorN` polynomials -evaluated at the vector `vals`, and write results to `jac`. If `vals` is ommited, it is evaluated at zero. +evaluated at the vector `vals`, and write results to `jac`. If `vals` is ommited, +it is evaluated at zero. """ -function jacobian!{T<:Number}(jac::Array{T,2}, vf::Array{TaylorN{T},1}) +function jacobian!(jac::Array{T,2}, vf::Array{TaylorN{T},1}) where {T<:Number} numVars = get_numvars() @assert length(vf) == numVars @assert (numVars, numVars) == size(jac) for comp2 = 1:numVars for comp1 = 1:numVars - @inbounds jac[comp1,comp2] = vf[comp1][2][comp2] + @inbounds jac[comp1,comp2] = vf[comp1][1][comp2] end end nothing end -function jacobian!{T<:Number}(jac::Array{T,2}, vf::Array{TaylorN{T},1},vals::Array{T,1}) +function jacobian!(jac::Array{T,2}, vf::Array{TaylorN{T},1}, + vals::Array{T,1}) where {T<:Number} + numVars = get_numvars() @assert length(vf) == numVars == length(vals) @assert (numVars, numVars) == size(jac) @@ -204,9 +303,10 @@ Return the hessian matrix (jacobian of the gradient) of `f::TaylorN`, evaluated at the vector `vals`. If `vals` is ommited, it is evaluated at zero. """ -hessian{T<:Number,S<:Number}(f::TaylorN{T}, vals::Array{S,1}) = +hessian(f::TaylorN{T}, vals::Array{S,1}) where {T<:Number, S<:Number} = (R = promote_type(T,S); jacobian( gradient(f), vals::Array{R,1}) ) -hessian{T<:Number}(f::TaylorN{T}) = hessian( f, zeros(T, get_numvars()) ) + +hessian(f::TaylorN{T}) where {T<:Number} = hessian( f, zeros(T, get_numvars()) ) """ ``` @@ -218,7 +318,81 @@ Return the hessian matrix (jacobian of the gradient) of `f::TaylorN`, evaluated at the vector `vals`, and write results to `hes`. If `vals` is ommited, it is evaluated at zero. """ -hessian!{T<:Number}(hes::Array{T,2}, f::TaylorN{T}, vals::Array{T,1}) = jacobian!(hes, gradient(f), vals) -hessian!{T<:Number}(hes::Array{T,2}, f::TaylorN{T}) = jacobian!(hes, gradient(f)) +hessian!(hes::Array{T,2}, f::TaylorN{T}, vals::Array{T,1}) where {T<:Number} = + jacobian!(hes, gradient(f), vals) + +hessian!(hes::Array{T,2}, f::TaylorN{T}) where {T<:Number} = + jacobian!(hes, gradient(f)) + + +##Integration +""" + integrate(a, r) + +Integrate the `a::HomogeneousPolynomial` with respect to the `r`-th +variable. The returned `HomogeneousPolynomial` has no added constant of +integration. If the order of a corresponds to `get_order()`, a zero +`HomogeneousPolynomial` of 0-th order is returned. + +""" +function integrate(a::HomogeneousPolynomial, r::Int) + @assert 1 ≤ r ≤ get_numvars() + + order_max = get_order() + T = promote_type(eltype(a), eltype(a[1]/1)) + a.order == order_max && return HomogeneousPolynomial(zero(T), 0) + + @inbounds posTb = pos_table[a.order+2] + @inbounds num_coeffs = size_table[a.order+1] + + coeffs = zeros(T, size_table[a.order+2]) + + @inbounds for i = 1:num_coeffs + iind = coeff_table[a.order+1][i] + n = iind[r] + n == order_max && continue + iind[r] += 1 + kdic = in_base(get_order(), iind) + pos = posTb[kdic] + coeffs[pos] = a[i] / (n+1) + iind[r] -= 1 + end + + return HomogeneousPolynomial{T}(coeffs, a.order+1) +end +integrate(a::HomogeneousPolynomial, s::Symbol) = integrate(a, lookupvar(s)) + + +""" + integrate(a, r, [x0]) + +Integrate the `a::TaylorN` series with respect to the `r`-th variable, +where `x0` the integration constant and must be independent +of the `r`-th variable; if `x0` is ommitted, it is taken as zero. +""" +function integrate(a::TaylorN, r::Int) + T = promote_type(eltype(a), eltype(a[0]/1)) + order_max = min(get_order(), a.order+1) + coeffs = zeros(HomogeneousPolynomial{T}, order_max) + + @inbounds for ord = 0:order_max-1 + coeffs[ord+1] = integrate( a[ord], r) + end + + return TaylorN(coeffs) +end +function integrate(a::TaylorN, r::Int, x0::TaylorN) + # Check constant of integration is independent of re + @assert derivative(x0, r) == 0.0 """ + The integration constant ($x0) must be independent of the + $(_params_TaylorN_.variable_names[r]) variable""" + + res = integrate(a, r) + return x0+res +end +integrate(a::TaylorN, r::Int, x0::NumberNotSeries) = + integrate(a,r,TaylorN(HomogeneousPolynomial([convert(eltype(a),x0)], 0))) -## TODO: Integration... +integrate(a::TaylorN, s::Symbol) = integrate(a, lookupvar(s)) +integrate(a::TaylorN, s::Symbol, x0::TaylorN) = integrate(a, lookupvar(s), x0) +integrate(a::TaylorN, s::Symbol, x0::NumberNotSeries) = integrate(a, lookupvar(s), x0) diff --git a/src/constructors.jl b/src/constructors.jl index f0724815..8d610d28 100644 --- a/src/constructors.jl +++ b/src/constructors.jl @@ -12,46 +12,49 @@ Parameterized abstract type for [`Taylor1`](@ref), [`HomogeneousPolynomial`](@ref) and [`TaylorN`](@ref). """ -@compat abstract type AbstractSeries{T<:Number} <: Number end +abstract type AbstractSeries{T<:Number} <: Number end ## Constructors ## ######################### Taylor1 -doc""" +""" Taylor1{T<:Number} <: AbstractSeries{T} DataType for polynomial expansions in one independent variable. **Fields:** -- `coeffs :: Array{T,1}` Expansion coefficients; the $i$-th - component is the coefficient of degree $i-1$ of the expansion. +- `coeffs :: Array{T,1}` Expansion coefficients; the ``i``-th + component is the coefficient of degree ``i-1`` of the expansion. - `order :: Int64` Maximum order (degree) of the polynomial. + +Note that `Taylor1` variables are callable. For more information, see +[`evaluate`](@ref). """ -immutable Taylor1{T<:Number} <: AbstractSeries{T} +struct Taylor1{T<:Number} <: AbstractSeries{T} coeffs :: Array{T,1} order :: Int ## Inner constructor ## - function (::Type{Taylor1{T}}){T}(coeffs::Array{T,1}, order::Int) - (order != length(coeffs)-1) && resize_coeffs1!(coeffs, order) + function Taylor1{T}(coeffs::Array{T,1}, order::Int) where T<:Number + resize_coeffs1!(coeffs, order) return new{T}(coeffs, order) end end ## Outer constructors ## -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(coeffs, length(coeffs)-1) -function Taylor1{T<:Number}(x::T, order::Int) +Taylor1(x::Taylor1{T}) where {T<:Number} = x +Taylor1(coeffs::Array{T,1}, order::Int) where {T<:Number} = Taylor1{T}(coeffs, order) +Taylor1(coeffs::Array{T,1}) where {T<:Number} = Taylor1(coeffs, length(coeffs)-1) +function Taylor1(x::T, order::Int) where {T<:Number} v = zeros(T, order+1) v[1] = x return Taylor1(v, order) end # Shortcut to define Taylor1 independent variables -doc""" +""" Taylor1([T::Type=Float64], [order::Int=1]) Shortcut to define the independent variable of a `Taylor1{T}` polynomial of @@ -65,12 +68,12 @@ julia> Taylor1(Rational{Int}, 4) 1//1 t + 𝒪(t⁵) ``` """ -Taylor1{T<:Number}(::Type{T}, order::Int64=1) = Taylor1{T}( [zero(T), one(T)], order) +Taylor1(::Type{T}, order::Int64=1) where {T<:Number} = Taylor1( [zero(T), one(T)], order) Taylor1(order::Int=1) = Taylor1(Float64, order) ######################### HomogeneousPolynomial -doc""" +""" HomogeneousPolynomial{T<:Number} <: AbstractSeries{T} DataType for homogenous polynomials in many (>1) independent variables. @@ -78,27 +81,30 @@ DataType for homogenous polynomials in many (>1) independent variables. **Fields:** - `coeffs :: Array{T,1}` Expansion coefficients of the homogeneous -polynomial; the $i$-th component is related to a monomial, where the degrees +polynomial; the ``i``-th component is related to a monomial, where the degrees of the independent variables are specified by `coeff_table[order+1][i]`. - `order :: Int` order (degree) of the homogenous polynomial. + +Note that `HomogeneousPolynomial` variables are callable. For more information, +see [`evaluate`](@ref). """ -immutable HomogeneousPolynomial{T<:Number} <: AbstractSeries{T} +struct HomogeneousPolynomial{T<:Number} <: AbstractSeries{T} coeffs :: Array{T,1} order :: Int - function (::Type{HomogeneousPolynomial{T}}){T}(coeffs::Array{T,1}, order::Int) + function HomogeneousPolynomial{T}(coeffs::Array{T,1}, order::Int) where T<:Number resize_coeffsHP!(coeffs, order) return new{T}(coeffs, order) end end -HomogeneousPolynomial{T<:Number}(x::HomogeneousPolynomial{T}) = x -HomogeneousPolynomial{T<:Number}(coeffs::Array{T,1}, order::Int) = +HomogeneousPolynomial(x::HomogeneousPolynomial{T}) where {T<:Number} = x +HomogeneousPolynomial(coeffs::Array{T,1}, order::Int) where {T<:Number} = HomogeneousPolynomial{T}(coeffs, order) -HomogeneousPolynomial{T<:Number}(coeffs::Array{T,1}) = - HomogeneousPolynomial{T}(coeffs, orderH(coeffs)) -HomogeneousPolynomial{T<:Number}(x::T, order::Int) = - HomogeneousPolynomial{T}([x], order) +HomogeneousPolynomial(coeffs::Array{T,1}) where {T<:Number} = + HomogeneousPolynomial(coeffs, orderH(coeffs)) +HomogeneousPolynomial(x::T, order::Int) where {T<:Number} = + HomogeneousPolynomial([x], order) # Shortcut to define HomogeneousPolynomial independent variable """ @@ -115,7 +121,7 @@ julia> HomogeneousPolynomial(Rational{Int}, 2) 1//1 x₂ ``` """ -function HomogeneousPolynomial{T<:Number}(::Type{T}, nv::Int) +function HomogeneousPolynomial(::Type{T}, nv::Int) where {T<:Number} @assert 0 < nv ≤ get_numvars() v = zeros(T, get_numvars()) @inbounds v[nv] = one(T) @@ -126,7 +132,7 @@ HomogeneousPolynomial(nv::Int) = HomogeneousPolynomial(Float64, nv) ######################### TaylorN -doc""" +""" TaylorN{T<:Number} <: AbstractSeries{T} DataType for polynomial expansions in many (>1) independent variables. @@ -134,15 +140,18 @@ DataType for polynomial expansions in many (>1) independent variables. **Fields:** - `coeffs :: Array{HomogeneousPolynomial{T},1}` Vector containing the -`HomogeneousPolynomial` entries. The $i$-th component corresponds to the -homogeneous polynomial of degree $i-1$. +`HomogeneousPolynomial` entries. The ``i``-th component corresponds to the +homogeneous polynomial of degree ``i-1``. - `order :: Int` maximum order of the polynomial expansion. + +Note that `TaylorN` variables are callable. For more information, see +[`evaluate`](@ref). """ -immutable TaylorN{T<:Number} <: AbstractSeries{T} +struct TaylorN{T<:Number} <: AbstractSeries{T} coeffs :: Array{HomogeneousPolynomial{T},1} order :: Int - function (::Type{TaylorN{T}}){T}(v::Array{HomogeneousPolynomial{T},1}, order::Int) + function TaylorN{T}(v::Array{HomogeneousPolynomial{T},1}, order::Int) where T<:Number m = maxorderH(v) order = max( m, order ) coeffs = zeros(HomogeneousPolynomial{T}, order) @@ -154,15 +163,15 @@ immutable TaylorN{T<:Number} <: AbstractSeries{T} end end -TaylorN{T<:Number}(x::TaylorN{T}) = x -TaylorN{T<:Number}(x::Array{HomogeneousPolynomial{T},1}, order::Int) = +TaylorN(x::TaylorN{T}) where T<:Number = x +TaylorN(x::Array{HomogeneousPolynomial{T},1}, order::Int) where {T<:Number} = TaylorN{T}(x, order) -TaylorN{T<:Number}(x::Array{HomogeneousPolynomial{T},1}) = TaylorN{T}(x,0) -TaylorN{T<:Number}(x::HomogeneousPolynomial{T},order::Int) = - TaylorN{T}([x], order) -TaylorN{T<:Number}(x::HomogeneousPolynomial{T}) = TaylorN{T}([x], x.order) -TaylorN{T<:Number}(x::T, order::Int) = - TaylorN{T}([HomogeneousPolynomial([x], 0)], order) +TaylorN(x::Array{HomogeneousPolynomial{T},1}) where {T<:Number} = TaylorN(x,0) +TaylorN(x::HomogeneousPolynomial{T},order::Int) where {T<:Number} = + TaylorN([x], order) +TaylorN(x::HomogeneousPolynomial{T}) where {T<:Number} = TaylorN([x], x.order) +TaylorN(x::T, order::Int) where {T<:Number} = + TaylorN([HomogeneousPolynomial([x], 0)], order) # Shortcut to define TaylorN independent variables """ @@ -181,7 +190,7 @@ julia> TaylorN(Rational{Int},2) 1//1 x₂ + 𝒪(‖x‖⁷) ``` """ -TaylorN{T<:Number}(::Type{T}, nv::Int; order::Int=get_order()) = +TaylorN(::Type{T}, nv::Int; order::Int=get_order()) where {T<:Number} = return TaylorN( [HomogeneousPolynomial(T, nv)], order ) TaylorN(nv::Int; order::Int=get_order()) = TaylorN(Float64, nv, order=order) @@ -189,6 +198,6 @@ TaylorN(nv::Int; order::Int=get_order()) = TaylorN(Float64, nv, order=order) # A `Number` which is not an `AbstractSeries` const NumberNotSeries = Union{setdiff(subtypes(Number), [AbstractSeries])...} -# A `Number` which is not a TaylorN nor a HomogeneousPolynomial +# A `Number` which is not `TaylorN` nor a `HomogeneousPolynomial` const NumberNotSeriesN = Union{setdiff(subtypes(Number), [AbstractSeries])..., Taylor1} diff --git a/src/conversion.jl b/src/conversion.jl index 3b44d4b2..9692b2df 100644 --- a/src/conversion.jl +++ b/src/conversion.jl @@ -7,112 +7,118 @@ # ## Conversion -convert{T<:Number}(::Type{Taylor1{T}}, a::Taylor1) = +convert(::Type{Taylor1{T}}, a::Taylor1) where {T<:Number} = Taylor1(convert(Array{T,1}, a.coeffs), a.order) -function convert{T<:Integer, S<:AbstractFloat}(::Type{Taylor1{Rational{T}}}, - a::Taylor1{S}) +function convert(::Type{Taylor1{Rational{T}}}, a::Taylor1{S}) where + {T<:Integer, S<:AbstractFloat} + la = length(a.coeffs) - v = Array{Rational{T}}(la) - v .= rationalize.(a[1:la]) + @compat v = Array{Rational{T}}(undef, la) + v .= rationalize.(a[0:la-1]) return Taylor1(v) end -convert{T<:Number}(::Type{Taylor1{T}}, b::Array{T,1}) = Taylor1(b, length(b)-1) +convert(::Type{Taylor1{T}}, b::Array{T,1}) where {T<:Number} = + Taylor1(b, length(b)-1) -convert{T<:Number, S<:Number}(::Type{Taylor1{T}}, b::Array{S,1}) = +convert(::Type{Taylor1{T}}, b::Array{S,1}) where {T<:Number, S<:Number} = Taylor1(convert(Array{T,1},b)) -convert{T<:Number, S<:Number}(::Type{Taylor1{T}}, b::S) = Taylor1([convert(T,b)], 0) +convert(::Type{Taylor1{T}}, b::S) where {T<:Number, S<:Number} = + Taylor1([convert(T,b)], 0) + +convert(::Type{Taylor1{T}}, b::T) where {T<:Number} = Taylor1([b], 0) -convert{T<:Number}(::Type{Taylor1{T}}, b::T) = Taylor1([b], 0) +convert(::Type{Taylor1}, a::T) where {T<:Number} = Taylor1(a, 0) -convert(::Type{Taylor1}, a::Number) = Taylor1(a, 0) +convert(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) where {T<:Number} = + HomogeneousPolynomial(convert(Array{T,1}, a.coeffs), a.order) -convert{T<:Number}(::Type{HomogeneousPolynomial{T}}, a::HomogeneousPolynomial) = - HomogeneousPolynomial{T}(convert(Array{T,1}, a.coeffs), a.order) +function convert(::Type{HomogeneousPolynomial{Rational{T}}}, + a::HomogeneousPolynomial{S}) where {T<:Integer, S<:AbstractFloat} -function convert{T<:Integer, S<:AbstractFloat}( - ::Type{HomogeneousPolynomial{Rational{T}}}, a::HomogeneousPolynomial{S}) la = length(a.coeffs) - v = Array{Rational{T}}(la) + @compat v = Array{Rational{T}}(undef, la) v .= rationalize.(a[1:la], tol=eps(one(S))) return HomogeneousPolynomial(v, a.order) end -convert{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}}, b::Array{S,1}) = - HomogeneousPolynomial{T}(convert(Array{T,1}, b), orderH(b)) -convert{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}}, b::S) = - HomogeneousPolynomial{T}([convert(T,b)], 0) +convert(::Type{HomogeneousPolynomial{T}}, b::Array{S,1}) where {T<:Number, S<:Number} = + HomogeneousPolynomial(convert(Array{T,1}, b), orderH(b)) + +convert(::Type{HomogeneousPolynomial{T}}, b::S) where {T<:Number, S<:Number}= + HomogeneousPolynomial([convert(T,b)], 0) -convert{T<:Number}(::Type{HomogeneousPolynomial{T}}, b::Array{T,1}) = - HomogeneousPolynomial{T}(b, orderH(b)) +convert(::Type{HomogeneousPolynomial{T}}, b::Array{T,1}) where {T<:Number} = + HomogeneousPolynomial(b, orderH(b)) -convert{T<:Number}(::Type{HomogeneousPolynomial{T}}, b::T) = - HomogeneousPolynomial{T}([b], 0) +convert(::Type{HomogeneousPolynomial{T}}, b::T) where {T<:Number} = + HomogeneousPolynomial([b], 0) convert(::Type{HomogeneousPolynomial}, a::Number) = HomogeneousPolynomial([a],0) -convert{T<:Number}(::Type{TaylorN{T}}, a::TaylorN) = - TaylorN{T}( convert(Array{HomogeneousPolynomial{T},1}, a.coeffs), a.order) +convert(::Type{TaylorN{T}}, a::TaylorN) where {T<:Number} = + TaylorN( convert(Array{HomogeneousPolynomial{T},1}, a.coeffs), a.order) -convert{T<:Number, S<:Number}(::Type{TaylorN{T}}, b::HomogeneousPolynomial{S}) = - TaylorN{T}( [convert(HomogeneousPolynomial{T}, b)], b.order) +convert(::Type{TaylorN{T}}, b::HomogeneousPolynomial{S}) where {T<:Number, S<:Number} = + TaylorN( [convert(HomogeneousPolynomial{T}, b)], b.order) -convert{T<:Number, S<:Number}(::Type{TaylorN{T}}, - b::Array{HomogeneousPolynomial{S},1}) = - TaylorN{T}( convert(Array{HomogeneousPolynomial{T},1}, b), length(b)-1) +convert(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{S},1}) where {T<:Number, S<:Number} = + TaylorN( convert(Array{HomogeneousPolynomial{T},1}, b), length(b)-1) -convert{T<:Number, S<:Number}(::Type{TaylorN{T}}, b::S) = +convert(::Type{TaylorN{T}}, b::S) where {T<:Number, S<:Number} = TaylorN( [HomogeneousPolynomial([convert(T, b)], 0)], 0) -convert{T<:Number}(::Type{TaylorN{T}}, b::HomogeneousPolynomial{T}) = - TaylorN{T}( [b], b.order) +convert(::Type{TaylorN{T}}, b::HomogeneousPolynomial{T}) where {T<:Number} = + TaylorN( [b], b.order) -convert{T<:Number}(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{T},1}) = - TaylorN{T}( b, length(b)-1) +convert(::Type{TaylorN{T}}, b::Array{HomogeneousPolynomial{T},1}) where {T<:Number} = + TaylorN( b, length(b)-1) -convert{T<:Number}(::Type{TaylorN{T}}, b::T) = +convert(::Type{TaylorN{T}}, b::T) where {T<:Number} = TaylorN( [HomogeneousPolynomial([b], 0)], 0) convert(::Type{TaylorN}, b::Number) = TaylorN( [HomogeneousPolynomial([b], 0)], 0) -function convert{T<:Number}(::Type{TaylorN{Taylor1{T}}}, s::Taylor1{TaylorN{T}}) +function convert(::Type{TaylorN{Taylor1{T}}}, s::Taylor1{TaylorN{T}}) where {T<:Number} + orderN = get_order() r = zeros(HomogeneousPolynomial{Taylor1{T}}, orderN) v = zeros(T, s.order+1) - @inbounds for ordT in 1:s.order+1 - v[ordT] = one(T) - @inbounds for ordHP in eachindex(s[ordT].coeffs) + @inbounds for ordT in 0:s.order + v[ordT+1] = one(T) + @inbounds for ordHP in 0:s[ordT].order @inbounds for ic in eachindex(s[ordT][ordHP].coeffs) coef = s[ordT][ordHP][ic] - r[ordHP][ic] += Taylor1(coef*v) + r[ordHP+1][ic] += Taylor1( coef.*v ) end end - v[ordT] = zero(T) + v[ordT+1] = zero(T) end return TaylorN(r) end -function convert{T<:Number}(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) +function convert(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) where {T<:Number} + ordert = 0 for ordHP in eachindex(s.coeffs) - ordert = max(ordert, s[ordHP][1].order) + ordert = max(ordert, s[ordHP-1][1].order) end - vT = Array{TaylorN{T}}(ordert+1) + @compat vT = Array{TaylorN{T}}(undef, ordert+1) @inbounds for ordT in eachindex(vT) vT[ordT] = TaylorN(zero(T), s.order) end @inbounds for ordN in eachindex(s.coeffs) - vHP = HomogeneousPolynomial(zeros(T,length(s[ordN]))) - @inbounds for ihp in eachindex(s[ordN].coeffs) - @inbounds for ind in eachindex(s[ordN][ihp].coeffs) - c = s[ordN][ihp][ind] + vHP = HomogeneousPolynomial(zeros(T,length(s[ordN-1]))) + @inbounds for ihp in eachindex(s[ordN-1].coeffs) + @inbounds for ind in eachindex(s[ordN-1][ihp].coeffs) + c = s[ordN-1][ihp][ind-1] vHP[ihp] = c vT[ind] += TaylorN(vHP) vHP[ihp] = zero(T) @@ -122,18 +128,20 @@ function convert{T<:Number}(::Type{Taylor1{TaylorN{T}}}, s::TaylorN{Taylor1{T}}) return Taylor1(vT) end -function convert{T<:Number,N}(::Type{Array{TaylorN{Taylor1{T}},N}}, - s::Array{Taylor1{TaylorN{T}},N}) - v = Array{TaylorN{Taylor1{T}}}(size(s)) +function convert(::Type{Array{TaylorN{Taylor1{T}},N}}, + s::Array{Taylor1{TaylorN{T}},N}) where {T<:Number, N} + + @compat v = Array{TaylorN{Taylor1{T}}}(undef, size(s)) for ind in eachindex(s) v[ind] = convert(TaylorN{Taylor1{T}}, s[ind]) end return v end -function convert{T<:Number,N}(::Type{Array{Taylor1{TaylorN{T}},N}}, - s::Array{TaylorN{Taylor1{T}},N}) - v = Array{Taylor1{TaylorN{T}}}(size(s)) +function convert(::Type{Array{Taylor1{TaylorN{T}},N}}, + s::Array{TaylorN{Taylor1{T}},N}) where {T<:Number, N} + + @compat v = Array{Taylor1{TaylorN{T}}}(undef, size(s)) for ind in eachindex(s) v[ind] = convert(Taylor1{TaylorN{T}}, s[ind]) end @@ -143,40 +151,54 @@ end # Promotion -promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{Taylor1{T}}) = Taylor1{T} +promote_rule(::Type{Taylor1{T}}, ::Type{Taylor1{T}}) where {T<:Number} = Taylor1{T} + +promote_rule(::Type{Taylor1{T}}, ::Type{Taylor1{S}}) where {T<:Number, S<:Number} = + Taylor1{promote_type(T,S)} + +promote_rule(::Type{Taylor1{T}}, ::Type{Array{T,1}}) where {T<:Number} = Taylor1{T} -promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{Taylor1{S}}) = - Taylor1{promote_type(T, S)} +promote_rule(::Type{Taylor1{T}}, ::Type{Array{S,1}}) where {T<:Number, S<:Number} = + Taylor1{promote_type(T,S)} -promote_rule{T<:Number}(::Type{Taylor1{T}}, ::Type{Array{T,1}}) = Taylor1{T} +promote_rule(::Type{Taylor1{T}}, ::Type{T}) where {T<:Number} = Taylor1{T} -promote_rule{T<:Number, S<:Number}(::Type{Taylor1{T}}, ::Type{Array{S,1}}) = - Taylor1{promote_type(T, S)} +promote_rule(::Type{Taylor1{T}}, ::Type{S}) where {T<:Number, S<:Number} = + 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)} +promote_rule(::Type{HomogeneousPolynomial{T}}, + ::Type{HomogeneousPolynomial{S}}) where {T<:Number, S<:Number} = + HomogeneousPolynomial{promote_type(T,S)} +promote_rule(::Type{HomogeneousPolynomial{T}}, + ::Type{Array{S,1}}) where {T<:Number, S<:Number} = HomogeneousPolynomial{promote_type(T,S)} -promote_rule{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}}, - ::Type{HomogeneousPolynomial{S}}) = HomogeneousPolynomial{promote_type(T,S)} +promote_rule(::Type{HomogeneousPolynomial{T}}, ::Type{S}) where + {T<:Number, S<:NumberNotSeries} = HomogeneousPolynomial{promote_type(T,S)} -promote_rule{T<:Number, S<:Number}(::Type{HomogeneousPolynomial{T}}, - ::Type{Array{S,1}}) = HomogeneousPolynomial{promote_type(T, S)} -promote_rule{T<:Number,S<:Union{Real,Complex}}(::Type{HomogeneousPolynomial{T}}, - ::Type{S}) = HomogeneousPolynomial{promote_type(T,S)} +promote_rule(::Type{TaylorN{T}}, ::Type{TaylorN{S}}) where {T<:Number, S<:Number}= + TaylorN{promote_type(T,S)} +promote_rule(::Type{TaylorN{T}}, ::Type{HomogeneousPolynomial{S}}) where + {T<:Number, S<:Number} = TaylorN{promote_type(T,S)} -promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, ::Type{TaylorN{S}}) = - TaylorN{promote_type(T, S)} +promote_rule(::Type{TaylorN{T}}, ::Type{Array{HomogeneousPolynomial{S},1}}) where + {T<:Number, S<:Number} = TaylorN{promote_type(T,S)} -promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, - ::Type{HomogeneousPolynomial{S}}) = TaylorN{promote_type(T, S)} +promote_rule(::Type{TaylorN{T}}, ::Type{S}) where {T<:Number, S<:Number} = + TaylorN{promote_type(T,S)} -promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, - ::Type{Array{HomogeneousPolynomial{S},1}}) = TaylorN{promote_type(T, S)} -promote_rule{T<:Number, S<:Number}(::Type{TaylorN{T}}, ::Type{S}) = - TaylorN{promote_type(T, S)} +# Order may matter +promote_rule(::Type{S}, ::Type{T}) where {S<:NumberNotSeries, T<:AbstractSeries} = + promote_rule(T,S) + +if VERSION <= v"0.7.0-DEV" + promote_rule(::Type{S}, ::Type{T}) where + {S<:Irrational, T<:AbstractSeries} = promote_rule(T,S) +else + promote_rule(::Type{S}, ::Type{T}) where + {S<:AbstractIrrational, T<:AbstractSeries} = promote_rule(T,S) +end diff --git a/src/dictmutfunct.jl b/src/dictmutfunct.jl new file mode 100644 index 00000000..72a1affe --- /dev/null +++ b/src/dictmutfunct.jl @@ -0,0 +1,175 @@ +# This file is part of the TaylorSeries.jl Julia package, MIT license +# +# Luis Benet & David P. Sanders +# UNAM +# +# MIT Expat license +# + +#= +This file contains some dictionary and functions to build +the dictionaries `_dict_unary_calls` and `_dict_binary_calls` +which allow to call the internal mutating functions. This +may be used to improve memory usage, e.g., to construct the +jet-coefficients used to integrate ODEs. +=# + +""" +`_InternalMutFuncs` + +Contains parameters and expressions that allow a simple +programatic construction for calling the internal mutating +functions. +""" +struct _InternalMutFuncs + namef :: Symbol # internal name of the function + argsf :: NTuple # arguments + defexpr :: Expr # defining expr + auxexpr :: Expr # auxiliary expr +end + +# Constructor +function _InternalMutFuncs( namef::Vector ) + if length(namef) == 3 + return _InternalMutFuncs( namef[1], namef[2], namef[3], Expr(:nothing) ) + else + return _InternalMutFuncs( namef...) + end +end + + +""" +`_dict_binary_ops` + +`Dict{Symbol, Array{Any,1}}` with the information to +construct the `_InternalMutFuncs` related to binary +operations. + +The keys correspond to the function symbols. + +The arguments of the array are the function name (e.g. `add!`), a tuple +with the function arguments, and an `Expr` with the calling pattern. The +convention for the arguments of the functions and the calling pattern +is to use `:_res` for the (mutated) result, `:_arg1` and `_arg2` +for the required arguments, and `:_k` for the computed order +of `:_res`. + +""" +const _dict_binary_ops = Dict( + :+ => [:add!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 + _arg2)], + :- => [:subst!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 - _arg2)], + :* => [:mul!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 * _arg2)], + :/ => [:div!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 / _arg2)], + :^ => [:pow!, (:_res, :_arg1, :_arg2, :_k), :(_res = _arg1 ^ _arg2)], +); + +""" +`_dict_binary_ops` + +`Dict{Symbol, Array{Any,1}}` with the information to +construct the `_InternalMutFuncs` related to unary +operations. + +The keys correspond to the function symbols. + +The arguments of the array are the function name (e.g. `add!`), a tuple +with the function arguments, and an `Expr` with the calling pattern. The +convention for the arguments of the functions and the calling pattern +is to use `:_res` for the (mutated) result, `:_arg1`, for the required +argument, possibly `:_aux` when there is an auxiliary expression +needed, and `:_k` for the computed order of `:_res`. When an auxiliary +expression is required, and `Expr` defining its calling pattern is +added as the last entry of the vector. + +""" +const _dict_unary_ops = Dict( + :+ => [:add!, (:_res, :_arg1, :_k), :(_res = + _arg1)], + :- => [:subst!, (:_res, :_arg1, :_k), :(_res = - _arg1)], + :sqr => [:sqr!, (:_res, :_arg1, :_k), :(_res = sqr(_arg1))], + :sqrt => [:sqrt!, (:_res, :_arg1, :_k), :(_res = sqrt(_arg1))], + :exp => [:exp!, (:_res, :_arg1, :_k), :(_res = exp(_arg1))], + :log => [:log!, (:_res, :_arg1, :_k), :(_res = log(_arg1))], + :identity => [:identity!, (:_res, :_arg1, :_k), :(_res = identity(_arg1))], + :zero => [:zero!, (:_res, :_arg1, :_k), :(_res = zero(_arg1))], + :one => [:one!, (:_res, :_arg1, :_k), :(_res = one(_arg1))], + :abs => [:abs!, (:_res, :_arg1, :_k), :(_res = abs(_arg1))], + :abs2 => [:abs2!, (:_res, :_arg1, :_k), :(_res = abs2(_arg1))], + :deg2rad => [:deg2rad!, (:_res, :_arg1, :_k), :(_res = deg2rad(_arg1))], + :rad2deg => [:rad2deg!, (:_res, :_arg1, :_k), :(_res = rad2deg(_arg1))], + # + :sin => [:sincos!, (:_res, :_aux, :_arg1, :_k), :(_res = sin(_arg1)), + :(_aux = cos(_arg1))], + :cos => [:sincos!, (:_aux, :_res, :_arg1, :_k), :(_res = cos(_arg1)), + :(_aux = sin(_arg1))], + :tan => [:tan!, (:_res, :_arg1, :_aux, :_k), :(_res = tan(_arg1)), + :(_aux = tan(_arg1)^2)], + :asin => [:asin!, (:_res, :_arg1, :_aux, :_k), :(_res = asin(_arg1)), + :(_aux = sqrt(1 - _arg1^2))], + :acos => [:acos!, (:_res, :_arg1, :_aux, :_k), :(_res = acos(_arg1)), + :(_aux = sqrt(1 - _arg1^2))], + :atan => [:atan!, (:_res, :_arg1, :_aux, :_k), :(_res = atan(_arg1)), + :(_aux = 1 + _arg1^2)], + :sinh => [:sinhcosh!, (:_res, :_aux, :_arg1, :_k), :(_res = sinh(_arg1)), + :(_aux = cosh(_arg1))], + :cosh => [:sinhcosh!, (:_aux, :_res, :_arg1, :_k), :(_res = cosh(_arg1)), + :(_aux = sinh(_arg1))], + :tanh => [:tanh!, (:_res, :_arg1, :_aux, :_k), :(_res = tanh(_arg1)), + :(_aux = tanh(_arg1)^2)], +); + + + +""" +``` +_internalmutfunc_call( fn :: _InternalMutFuncs ) +``` + +Creates the appropriate call to the internal mutating +function defined by the `_InternalMutFuncs` object. +This is used to construct [`_dict_unary_calls`](@ref) +and [`_dict_binary_calls`](@ref). +The call contains the prefix `TaylorSeries.`. +""" +@compat _internalmutfunc_call( fn :: _InternalMutFuncs ) = ( + Expr( :call, Meta.parse("TaylorSeries.$(fn.namef)"), fn.argsf... ), fn.defexpr, fn.auxexpr ) + + + +""" +`_dict_unary_calls::Dict{Symbol, NTuple{2,Expr}}` + +Dictionary with the expressions that define the +internal unary functions and the auxiliary functions, +whenever they exist. The keys correspond to those +functions, passed as symbols, with the defined +internal mutating functions. + +Evaluating the entries generates expressions that represent +the actual calls to the internal mutating functions. +""" +const _dict_unary_calls = Dict{Symbol, NTuple{3,Expr}}() + +#Populates the constant vector `_dict_unary_calls`. +for kk in keys(_dict_unary_ops) + res = _internalmutfunc_call( _InternalMutFuncs(_dict_unary_ops[kk]) ) + push!(_dict_unary_calls, kk => res ) +end + +""" +`_dict_binary_calls::Dict{Symbol, NTuple{2,Expr}}` + +Dictionary with the expressions that define the +internal binary functions and the auxiliary functions, +whenever they exist. The keys correspond to those +functions, passed as symbols, with the defined +internal mutating functions. + +Evaluating the entries generates symbols that represent +the actual calls to the internal mutating functions. +""" +const _dict_binary_calls = Dict{Symbol, NTuple{3,Expr}}() +#Populates the constant vector `_dict_binary_calls`. +for kk in keys(_dict_binary_ops) + res = _internalmutfunc_call( _InternalMutFuncs(_dict_binary_ops[kk]) ) + push!(_dict_binary_calls, kk => res ) +end diff --git a/src/evaluate.jl b/src/evaluate.jl index adf26561..96ae4402 100644 --- a/src/evaluate.jl +++ b/src/evaluate.jl @@ -11,43 +11,67 @@ evaluate(a, [dx]) Evaluate a `Taylor1` polynomial using Horner's rule (hand coded). If `dx` is -ommitted, its value is considered as zero. +ommitted, its value is considered as zero. Note that the syntax `a(dx)` is +equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`. """ -function evaluate{T<:Number}(a::Taylor1{T}, dx::T) +function evaluate(a::Taylor1{T}, dx::T) where {T<:Number} @inbounds suma = a[end] - @inbounds for k = a.order:-1:1 + @inbounds for k = a.order-1:-1:0 suma = suma*dx + a[k] end suma end -function evaluate{T<:Number,S<:Number}(a::Taylor1{T}, dx::S) +function evaluate(a::Taylor1{T}, dx::S) where {T<:Number, S<:Number} R = promote_type(T,S) @inbounds suma = convert(R, a[end]) - @inbounds for k = a.order:-1:1 + @inbounds for k = a.order-1:-1:0 suma = suma*dx + a[k] end suma end -evaluate{T<:Number}(a::Taylor1{T}) = a[1] +evaluate(a::Taylor1{T}) where {T<:Number} = a[0] -doc""" +""" evaluate(x, δt) -Evaluates each element of `x::Array{Taylor1{T},1}`, representing -the dependent variables of an ODE, at *time* δt. +Evaluates each element of `x::Union{ Vector{Taylor1{T}}, Matrix{Taylor1{T}} }`, +representing the dependent variables of an ODE, at *time* δt. Note that the +syntax `x(δt)` is equivalent to `evaluate(x, δt)`, and `x()` +is equivalent to `evaluate(x)`. """ -function evaluate{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S) +function evaluate(x::Union{Array{Taylor1{T},1}, SubArray{Taylor1{T},1}}, δt::S) where {T<:Number, S<:Number} R = promote_type(T,S) return evaluate(convert(Array{Taylor1{R},1},x), convert(R,δt)) end -function evaluate{T<:Number}(x::Array{Taylor1{T},1}, δt::T) - xnew = Array{T}( length(x) ) +function evaluate(x::Array{Taylor1{T},1}, δt::T) where {T<:Number} + @compat xnew = Array{T}(undef, length(x) ) evaluate!(x, δt, xnew) return xnew end -evaluate{T<:Number}(a::Array{Taylor1{T},1}) = evaluate(a, zero(T)) -doc""" +evaluate(a::Array{Taylor1{T},1}) where {T<:Number} = evaluate(a, zero(T)) +evaluate(a::SubArray{Taylor1{T},1}) where {T<:Number} = evaluate(a, zero(T)) + +function evaluate(A::Union{Array{Taylor1{T},2}, SubArray{Taylor1{T},2}}, δt::S) where {T<:Number, S<:Number} + R = promote_type(T,S) + return evaluate(convert(Array{Taylor1{R},2},A), convert(R,δt)) +end +function evaluate(A::Array{Taylor1{T},2}, δt::T) where {T<:Number} + n,m = size(A) + @compat Anew = Array{T}(undef, n, m ) + @compat xnew = Array{T}(undef, n ) + + for i in 1:m + evaluate!(A[:,i], δt, xnew) + Anew[:,i] = xnew + end + + return Anew +end +evaluate(A::Array{Taylor1{T},2}) where {T<:Number} = evaluate.(A) +evaluate(A::SubArray{Taylor1{T},2}) where {T<:Number} = evaluate.(A) + +""" evaluate!(x, δt, x0) Evaluates each element of `x::Array{Taylor1{T},1}`, @@ -55,14 +79,18 @@ representing the Taylor expansion for the dependent variables of an ODE at *time* `δt`. It updates the vector `x0` with the computed values. """ -function evaluate!{T<:Number}(x::Array{Taylor1{T},1}, δt::T, x0::Union{Array{T,1},SubArray{T,1}}) +function evaluate!(x::Array{Taylor1{T},1}, δt::T, + x0::Union{Array{T,1},SubArray{T,1}}) where {T<:Number} + @assert length(x) == length(x0) @inbounds for i in eachindex(x) x0[i] = evaluate( x[i], δt ) end nothing end -function evaluate!{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S, x0::Union{Array{T,1},SubArray{T,1}}) +function evaluate!(x::Array{Taylor1{T},1}, δt::S, + x0::Union{Array{T,1},SubArray{T,1}}) where {T<:Number, S<:Number} + @assert length(x) == length(x0) @inbounds for i in eachindex(x) x0[i] = evaluate( x[i], δt ) @@ -70,38 +98,88 @@ function evaluate!{T<:Number, S<:Number}(x::Array{Taylor1{T},1}, δt::S, x0::Uni nothing end - """ evaluate(a, x) Substitute `x::Taylor1` as independent variable in a `a::Taylor1` polynomial. +Note that the syntax `a(x)` is equivalent to `evaluate(a, x)`. """ -evaluate{T<:Number,S<:Number}(a::Taylor1{T}, x::Taylor1{S}) = +evaluate(a::Taylor1{T}, x::Taylor1{S}) where {T<:Number, S<:Number} = evaluate(promote(a,x)...) -function evaluate{T<:Number}(a::Taylor1{T}, x::Taylor1{T}) + +function evaluate(a::Taylor1{T}, x::Taylor1{T}) where {T<:Number} if a.order != x.order a, x = fixorder(a, x) end @inbounds suma = a[end] - @inbounds for k = a.order:-1:1 + @inbounds for k = a.order-1:-1:0 suma = suma*x + a[k] end suma end +function evaluate(a::Taylor1{Taylor1{T}}, x::Taylor1{T}) where {T<:Number} + @inbounds suma = a[end] + @inbounds for k = a.order-1:-1:0 + suma = suma*x + a[k] + end + suma +end + +evaluate(p::Taylor1{T}, x::Array{S}) where {T<:Number, S<:Number} = + evaluate.(p, x) + +#function-like behavior for Taylor1 +(p::Taylor1)(x) = evaluate(p, x) + +(p::Taylor1)() = evaluate(p) +#function-like behavior for Vector{Taylor1} +(p::Array{Taylor1{T},1})(x) where {T<:Number} = evaluate(p, x) +(p::SubArray{Taylor1{T},1})(x) where {T<:Number} = evaluate(p, x) +(p::Array{Taylor1{T},1})() where {T<:Number} = evaluate.(p) +(p::SubArray{Taylor1{T},1})() where {T<:Number} = evaluate.(p) + +#function-like behavior for Matrix{Taylor1} +(p::Array{Taylor1{T},2})(x) where {T<:Number} = evaluate(p, x) +(p::SubArray{Taylor1{T},2})(x) where {T<:Number} = evaluate(p, x) +(p::Array{Taylor1{T},2})() where {T<:Number} = evaluate.(p) +(p::SubArray{Taylor1{T},2})() where {T<:Number} = evaluate.(p) ## Evaluation of multivariable -function evaluate!{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1}, - x0::Array{T,1}) +function evaluate!(x::Array{TaylorN{T},1}, δx::Array{T,1}, + x0::Array{T,1}) where {T<:Number} + @assert length(x) == length(x0) @inbounds for i in eachindex(x) x0[i] = evaluate( x[i], δx ) end nothing end -function evaluate!{T<:Number}(x::Array{Taylor1{TaylorN{T}},1}, δt::T, - x0::Array{TaylorN{T},1}) + +function evaluate!(x::Array{TaylorN{T},1}, δx::Array{Taylor1{T},1}, + x0::Array{Taylor1{T},1}) where {T<:NumberNotSeriesN} + + @assert length(x) == length(x0) + @inbounds for i in eachindex(x) + x0[i] = evaluate( x[i], δx ) + end + nothing +end + +function evaluate!(x::Array{TaylorN{T},1}, δx::Array{TaylorN{T},1}, + x0::Array{TaylorN{T},1}) where {T<:NumberNotSeriesN} + + @assert length(x) == length(x0) + @inbounds for i in eachindex(x) + x0[i] = evaluate( x[i], δx ) + end + nothing +end + +function evaluate!(x::Array{Taylor1{TaylorN{T}},1}, δt::T, + x0::Array{TaylorN{T},1}) where {T<:Number} + @assert length(x) == length(x0) @inbounds for i in eachindex(x) x0[i] = evaluate( x[i], δt ) @@ -110,53 +188,73 @@ function evaluate!{T<:Number}(x::Array{Taylor1{TaylorN{T}},1}, δt::T, end """ - evaluate(a, vals) + evaluate(a, [vals]) -Evaluate a `HomogeneousPolynomial` polynomial at `vals`. +Evaluate a `HomogeneousPolynomial` polynomial at `vals`. If `vals` is ommitted, +it's evaluated at zero. Note that the syntax `a(vals)` is equivalent to +`evaluate(a, vals)`; and `a()` is equivalent to `evaluate(a)`. """ -function evaluate{T<:Number,S<:NumberNotSeriesN}(a::HomogeneousPolynomial{T}, - vals::Array{S,1} ) - @assert length(vals) == get_numvars() +function evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,S} ) where + {T<:Number, S<:NumberNotSeriesN, N} + + @assert N == get_numvars() - num_vars = get_numvars() ct = coeff_table[a.order+1] R = promote_type(T,S) suma = zero(R) - for (i,a_coeffs) in enumerate(a.coeffs) + for (i,a_coeff) in enumerate(a.coeffs) + iszero(a_coeff) && continue tmp = vals[1]^(ct[i][1]) - for n in 2:num_vars + for n in 2:N tmp *= vals[n]^(ct[i][n]) end - suma += a_coeffs * tmp + suma += a_coeff * tmp end return suma end -evaluate(a::HomogeneousPolynomial) = zero(a[1]) +evaluate(a::HomogeneousPolynomial{T}, vals::Array{S,1} ) where + {T<:Number, S<:NumberNotSeriesN} = evaluate(a, (vals...,)) + +evaluate(a::HomogeneousPolynomial, v, vals...) = evaluate(a, (v, vals...,)) + +function evaluate(a::HomogeneousPolynomial) + a.order == 0 && return a[1] + zero(a[1]) +end + +#function-like behavior for HomogeneousPolynomial +(p::HomogeneousPolynomial)(x) = evaluate(p, x) + +(p::HomogeneousPolynomial)(x, v...) = evaluate(p, (x, v...,)) + +(p::HomogeneousPolynomial)() = evaluate(p) """ evaluate(a, [vals]) Evaluate the `TaylorN` polynomial `a` at `vals`. If `vals` is ommitted, it's evaluated at zero. +Note that the syntax `a(vals)` is equivalent to `evaluate(a, vals)`; and `a()` +is equivalent to `evaluate(a)`. """ -function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T}, - vals::Array{S,1}) - @assert length(vals) == get_numvars() +function evaluate(a::TaylorN{T}, vals::NTuple{N,S}) where + {T<:Number,S<:NumberNotSeries, N} + + @assert N == get_numvars() - num_vars = get_numvars() - ct = coeff_table R = promote_type(T,S) a_length = length(a) - suma = zeros(R,a_length) + suma = zeros(R, a_length) for homPol in 1:length(a) sun = zero(R) - for (i,a_coeff) in enumerate(a.coeffs[homPol].coeffs) - tmp = vals[1]^(ct[homPol][i][1]) - for n in 2:num_vars - tmp *= vals[n]^(ct[homPol][i][n]) + for (i, a_coeff) in enumerate(a.coeffs[homPol].coeffs) + iszero(a_coeff) && continue + tmp = vals[1]^(coeff_table[homPol][i][1]) + for n in 2:N + tmp *= vals[n]^(coeff_table[homPol][i][n]) end sun += a_coeff * tmp end @@ -165,23 +263,82 @@ function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T}, return sum( sort!(suma, by=abs2) ) end -function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T}, - vals::Array{Taylor1{S},1}) - @assert length(vals) == get_numvars() - num_vars = get_numvars() - ct = coeff_table +evaluate(a::TaylorN{T}, vals::Array{S,1}) where + {T<:Number, S<:NumberNotSeriesN} = evaluate(a, (vals...,)) + +evaluate(a::TaylorN, v, vals...) = evaluate(a, (v, vals...,)) + +function evaluate(a::TaylorN{T}, vals::NTuple{N,Taylor1{S}}) where + {T<:Number, S<:NumberNotSeries, N} + + @assert N == get_numvars() + R = promote_type(T,S) a_length = length(a) ord = maximum( get_order.(vals) ) suma = Taylor1(zeros(R, ord)) for homPol in 1:length(a) - sun = zero(R) - for (i,a_coeff) in enumerate(a.coeffs[homPol].coeffs) - tmp = vals[1]^(ct[homPol][i][1]) + for (i, a_coeff) in enumerate(a.coeffs[homPol].coeffs) + iszero(a_coeff) && continue + tmp = vals[1]^(coeff_table[homPol][i][1]) + for n in 2:N + tmp *= vals[n]^(coeff_table[homPol][i][n]) + end + suma += a_coeff * tmp + end + end + + return suma +end + +evaluate(a::TaylorN{T}, vals::Array{Taylor1{S},1}) where + {T<:Number, S<:NumberNotSeriesN} = evaluate(a, (vals...,)) + +function evaluate(a::TaylorN{Taylor1{T}}, vals::NTuple{N, Taylor1{T}}) where + {T<:NumberNotSeries, N} + + @assert N == get_numvars() + + a_length = length(a) + ord = maximum( get_order.(vals) ) + suma = Taylor1(zeros(T, ord)) + + for homPol in 1:length(a) + for (i, a_coeff) in enumerate(a.coeffs[homPol].coeffs) + iszero(a_coeff) && continue + tmp = vals[1]^(coeff_table[homPol][i][1]) + for n in 2:N + tmp *= vals[n]^(coeff_table[homPol][i][n]) + end + suma += a_coeff * tmp + end + end + + return suma +end + +evaluate(a::TaylorN{Taylor1{T}}, vals::Array{Taylor1{T},1}) where + {T<:NumberNotSeries} = evaluate(a, (vals...,)) + +function evaluate(a::TaylorN{T}, vals::NTuple{N, TaylorN{S}}) where + {T<:Number, S<:NumberNotSeries, N} + + @assert length(vals) == get_numvars() + + num_vars = get_numvars() + R = promote_type(T,eltype(S)) + a_length = length(a) + ord = maximum( get_order.(vals) ) + suma = zero(TaylorN{R}) + + for homPol in 1:length(a) + for (i, a_coeff) in enumerate(a.coeffs[homPol].coeffs) + iszero(a_coeff) && continue + tmp = vals[1]^(coeff_table[homPol][i][1]) for n in 2:num_vars - tmp *= vals[n]^(ct[homPol][i][n]) + tmp *= vals[n]^(coeff_table[homPol][i][n]) end suma += a_coeff * tmp end @@ -190,10 +347,72 @@ function evaluate{T<:Number,S<:NumberNotSeries}(a::TaylorN{T}, return suma end -evaluate{T<:Number}(a::TaylorN{T}) = a[1][1] +evaluate(a::TaylorN{T}, vals::Array{TaylorN{S},1}) where + {T<:Number, S<:NumberNotSeries} = evaluate(a, (vals...,)) + +function evaluate(a::TaylorN{T}, s::Symbol, val::S) where + {T<:Number, S<:NumberNotSeriesN} + vars = get_variables() + ind = lookupvar(s) + vars[ind] = val + evaluate(a, vars) +end + +evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T<:Number, S<:NumberNotSeriesN} = + evaluate(a, first(x), last(x)) + +evaluate(a::TaylorN{T}) where {T<:Number} = a[0][1] -function evaluate{T<:Number}(x::Array{TaylorN{T},1}, δx::Array{T,1}) - x0 = Array{T}( length(x) ) +#Vector evaluation +function evaluate(x::Union{Array{TaylorN{T},1},SubArray{TaylorN{T},1}}, δx::Vector{S}) where {T<:Number, S<:Number} + R = promote_type(T,S) + return evaluate(convert(Array{TaylorN{R},1},x), convert(Vector{R},δx)) +end + +function evaluate(x::Array{TaylorN{T},1}, δx::Array{T,1}) where {T<:Number} + @compat x0 = Array{T}(undef, length(x) ) evaluate!( x, δx, x0 ) return x0 end + +evaluate(x::Array{TaylorN{T},1}) where {T<:Number} = evaluate.(x) +evaluate(x::SubArray{TaylorN{T},1}) where {T<:Number} = evaluate.(x) + +#Matrix evaluation +function evaluate(A::Union{Array{TaylorN{T},2}, SubArray{TaylorN{T},2}}, δx::Vector{S}) where {T<:Number, S<:Number} + R = promote_type(T,S) + return evaluate(convert(Array{TaylorN{R},2},A), convert(Vector{R},δx)) +end +function evaluate(A::Array{TaylorN{T},2}, δx::Vector{T}) where {T<:Number} + n,m = size(A) + @compat Anew = Array{T}(undef, n, m ) + @compat xnew = Array{T}(undef, n ) + + for i in 1:m + evaluate!(A[:,i], δx, xnew) + Anew[:,i] = xnew + end + + return Anew +end +evaluate(A::Array{TaylorN{T},2}) where {T<:Number} = evaluate.(A) +evaluate(A::SubArray{TaylorN{T},2}) where {T<:Number} = evaluate.(A) + +#function-like behavior for TaylorN +(p::TaylorN)(x) = evaluate(p, x) +(p::TaylorN)() = evaluate(p) +(p::TaylorN)(s::Symbol, x) = evaluate(p, s, x) +(p::TaylorN)(x::Pair) = evaluate(p, first(x), last(x)) +(p::TaylorN)(x, v...) = evaluate(p, (x, v...,)) + +#function-like behavior for Vector{TaylorN} +(p::Array{TaylorN{T},1})(x) where {T<:Number} = evaluate(p, x) +(p::SubArray{TaylorN{T},1})(x) where {T<:Number} = evaluate(p, x) +(p::Array{TaylorN{T},1})() where {T<:Number} = evaluate(p) +(p::SubArray{TaylorN{T},1})() where {T<:Number} = evaluate(p) + +#function-like behavior for Matrix{TaylorN} +(p::Array{TaylorN{T},2})(x) where {T<:Number} = evaluate(p, x) +(p::SubArray{TaylorN{T},2})(x) where {T<:Number} = evaluate(p, x) +(p::Array{TaylorN{T},2})() where {T<:Number} = evaluate.(p) +(p::SubArray{TaylorN{T},2})() where {T<:Number} = evaluate.(p) diff --git a/src/functions.jl b/src/functions.jl index 60fd23dc..c0cda296 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -12,24 +12,27 @@ for T in (:Taylor1, :TaylorN) @eval begin function exp(a::$T) order = a.order - c = $T( exp( constant_term(a) ), order ) + aux = exp(constant_term(a)) + aa = one(aux) * a + c = $T( aux, order ) for k = 1:order - exp!(c, a, k) + exp!(c, aa, k) end return c end function log(a::$T) - constant_term(a) == zero(constant_term(a)) && - throw(ArgumentError(""" + iszero(constant_term(a)) && throw(ArgumentError(""" The 0-th order `TaylorN` coefficient must be non-zero in order to expand `log` around 0. """)) order = a.order - c = $T( log( constant_term(a) ), order ) + aux = log(constant_term(a)) + aa = one(aux) * a + c = $T( aux, order ) for k = 1:order - log!(c, a, k) + log!(c, aa, k) end return c end @@ -38,10 +41,12 @@ for T in (:Taylor1, :TaylorN) cos(a::$T) = sincos(a)[2] function sincos(a::$T) order = a.order - s = $T( sin(constant_term(a)), order ) + aux = sin(constant_term(a)) + aa = one(aux) * a + s = $T( aux, order ) c = $T( cos(constant_term(a)), order ) for k = 1:order - sincos!(s, c, a, k) + sincos!(s, c, aa, k) end return s, c end @@ -49,10 +54,11 @@ for T in (:Taylor1, :TaylorN) function tan(a::$T) order = a.order aux = tan(constant_term(a)) + aa = one(aux) * a c = $T(aux, order) c2 = $T(aux^2, order) for k = 1:order - tan!(c, a, c2, k) + tan!(c, aa, c2, k) end return c end @@ -66,10 +72,12 @@ for T in (:Taylor1, :TaylorN) """)) order = a.order - c = $T( asin(a0), order ) + aux = asin(a0) + aa = one(aux) * a + c = $T( aux, order ) r = $T( sqrt(1 - a0^2), order ) for k in 1:order - asin!(c, a, r, k) + asin!(c, aa, r, k) end return c end @@ -82,32 +90,31 @@ for T in (:Taylor1, :TaylorN) in the denominator. """)) - # The following exploits a trick that the series of - # `acos(a)` is generated as the series of `asin(a)` except - # for a sign and the constant term; see `acos!`. order = a.order - c = $T( asin(a0), order ) + aux = acos(a0) + aa = one(aux) * a + c = $T( aux, order ) r = $T( sqrt(1 - a0^2), order ) for k in 1:order - acos!(c, a, r, k) + acos!(c, aa, r, k) end - @inbounds c[1] = acos(a0) return c end function atan(a::$T) order = a.order a0 = constant_term(a) - c = $T( atan(a0), order) - r = $T(1 + a0^2, order) - constant_term(r) == zero(constant_term(a)) && - throw(ArgumentError( + aux = atan(a0) + aa = one(aux) * a + c = $T( aux, order) + r = $T(one(aux) + a0^2, order) + iszero(constant_term(r)) && throw(ArgumentError( """ Recursion formula has a pole. """)) for k in 1:order - atan!(c, a, r, k) + atan!(c, aa, r, k) end return c end @@ -116,10 +123,12 @@ for T in (:Taylor1, :TaylorN) cosh(a::$T) = sinhcosh(a)[2] function sinhcosh(a::$T) order = a.order - s = $T( sinh(constant_term(a)), order) + aux = sinh(constant_term(a)) + aa = one(aux) * a + s = $T( aux, order) c = $T( cosh(constant_term(a)), order) for k = 1:order - sinhcosh!(s, c, a, k) + sinhcosh!(s, c, aa, k) end return s, c end @@ -127,10 +136,11 @@ for T in (:Taylor1, :TaylorN) function tanh(a::$T) order = a.order aux = tanh( constant_term(a) ) + aa = one(aux) * a c = $T( aux, order) c2 = $T( aux^2, order) for k = 1:order - tanh!(c, a, c2, k) + tanh!(c, aa, c2, k) end return c end @@ -141,176 +151,219 @@ end # Recursive functions (homogeneous coefficients) for T in (:Taylor1, :TaylorN) @eval begin - @inline function exp!(c::$T, a::$T, k::Int) + @inline function identity!(c::$T{T}, a::$T{T}, k::Int) where {T} + @inbounds c[k] = identity(a[k]) + return nothing + end + + @inline function zero!(c::$T{T}, a::$T{T}, k::Int) where{T} + @inbounds c[k] = zero(a[k]) + return nothing + end + + @inline function one!(c::$T{T}, a::$T{T}, k::Int) where{T} + if k == 0 + @inbounds c[0] = one(a[0]) + else + @inbounds c[k] = zero(a[k]) + end + return nothing + end + + @inline function abs!(c::$T{T}, a::$T{T}, k::Int) where {T} + z = zero(constant_term(a)) + if constant_term(constant_term(a)) > constant_term(z) + return add!(c, a, k) + elseif constant_term(constant_term(a)) < constant_term(z) + return subst!(c, a, k) + else + throw(ArgumentError( + """The 0th order coefficient must be non-zero + (abs(x) is not differentiable at x=0).""")) + end + return nothing + end + + @inline abs2!(c::$T{T}, a::$T{T}, k::Int) where {T} = sqr!(c, a, k) + + @inline function exp!(c::$T{T}, a::$T{T}, k::Int) where {T} if k == 0 - @inbounds c[1] = exp(constant_term(a)) + @inbounds c[0] = exp(constant_term(a)) return nothing end @inbounds for i = 0:k-1 if $T == Taylor1 - c[k+1] += (k-i) * a[k-i+1] * c[i+1] + c[k] += (k-i) * a[k-i] * c[i] else - mul!(c[k+1], (k-i) * a[k-i+1], c[i+1]) + mul!(c[k], (k-i) * a[k-i], c[i]) end end - @inbounds c[k+1] = c[k+1] / k + @inbounds c[k] = c[k] / k return nothing end - @inline function log!(c::$T, a::$T, k::Int) + @inline function log!(c::$T{T}, a::$T{T}, k::Int) where {T} if k == 0 - @inbounds c[1] = log(constant_term(a)) + @inbounds c[0] = log(constant_term(a)) return nothing end @inbounds for i = 1:k-1 if $T == Taylor1 - c[k+1] += (k-i) * a[i+1] * c[k-i+1] + c[k] += (k-i) * a[i] * c[k-i] else - mul!(c[k+1], (k-i)*a[i+1], c[k-i+1]) + mul!(c[k], (k-i)*a[i], c[k-i]) end end - @inbounds c[k+1] = (a[k+1] - c[k+1]/k) / constant_term(a) + @inbounds c[k] = (a[k] - c[k]/k) / constant_term(a) return nothing end - @inline function sincos!(s::$T, c::$T, a::$T, k::Int) + @inline function sincos!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T} if k == 0 a0 = constant_term(a) - @inbounds s[1], c[1] = sin( a0 ), cos( a0 ) + @inbounds s[0], c[0] = sin( a0 ), cos( a0 ) return nothing end @inbounds for i = 1:k - x = i * a[i+1] + x = i * a[i] if $T == Taylor1 - s[k+1] += x * c[k-i+1] - c[k+1] -= x * s[k-i+1] + s[k] += x * c[k-i] + c[k] -= x * s[k-i] else - mul!(s[k+1], x, c[k-i+1]) - mul!(c[k+1], -x, s[k-i+1]) + mul!(s[k], x, c[k-i]) + mul!(c[k], -x, s[k-i]) end end - @inbounds s[k+1] = s[k+1] / k - @inbounds c[k+1] = c[k+1] / k + @inbounds s[k] = s[k] / k + @inbounds c[k] = c[k] / k return nothing end - @inline function tan!(c::$T, a::$T, c2::$T, k::Int) + @inline function tan!(c::$T{T}, a::$T{T}, c2::$T{T}, k::Int) where {T} if k == 0 @inbounds aux = tan( constant_term(a) ) - @inbounds c[1] = aux - @inbounds c2[1] = aux^2 + @inbounds c[0] = aux + @inbounds c2[0] = aux^2 return nothing end @inbounds for i = 0:k-1 if $T == Taylor1 - c[k+1] += (k-i) * a[k-i+1] * c2[i+1] + c[k] += (k-i) * a[k-i] * c2[i] else - mul!(c[k+1], (k-i) * a[k-i+1], c2[i+1]) + mul!(c[k], (k-i) * a[k-i], c2[i]) end end - @inbounds c[k+1] = a[k+1] + c[k+1]/k + @inbounds c[k] = a[k] + c[k]/k sqr!(c2, c, k) return nothing end - @inline function asin!(c::$T, a::$T, r::$T, k::Int) + @inline function asin!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T} if k == 0 a0 = constant_term(a) - @inbounds c[1] = asin( a0 ) - @inbounds r[1] = sqrt( 1 - a0^2) + @inbounds c[0] = asin( a0 ) + @inbounds r[0] = sqrt( 1 - a0^2) return nothing end @inbounds for i in 1:k-1 if $T == Taylor1 - c[k+1] += (k-i) * r[i+1] * c[k-i+1] + c[k] += (k-i) * r[i] * c[k-i] else - mul!(c[k+1], (k-i) * r[i+1], c[k-i+1]) + mul!(c[k], (k-i) * r[i], c[k-i]) end end sqrt!(r, 1-a^2, k) - @inbounds c[k+1] = (a[k+1] - c[k+1]/k) / constant_term(r) + @inbounds c[k] = (a[k] - c[k]/k) / constant_term(r) return nothing end - @inline function acos!(c::$T, a::$T, r::$T, k::Int) + @inline function acos!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T} if k == 0 a0 = constant_term(a) - @inbounds c[1] = acos( a0 ) - @inbounds r[1] = sqrt( 1 - a0^2) + @inbounds c[0] = acos( a0 ) + @inbounds r[0] = sqrt( 1 - a0^2) return nothing end - asin!(c, -a, r, k) + @inbounds for i in 1:k-1 + if $T == Taylor1 + c[k] += (k-i) * r[i] * c[k-i] + else + mul!(c[k], (k-i) * r[i], c[k-i]) + end + end + sqrt!(r, 1-a^2, k) + @inbounds c[k] = -(a[k] + c[k]/k) / constant_term(r) return nothing end - @inline function atan!(c::$T, a::$T, r::$T, k::Int) + @inline function atan!(c::$T{T}, a::$T{T}, r::$T{T}, k::Int) where {T} if k == 0 a0 = constant_term(a) - @inbounds c[1] = atan( a0 ) - @inbounds r[1] = 1 + a0^2 + @inbounds c[0] = atan( a0 ) + @inbounds r[0] = 1 + a0^2 return nothing end @inbounds for i in 1:k-1 if $T == Taylor1 - c[k+1] += (k-i) * r[i+1] * c[k-i+1] + c[k] += (k-i) * r[i] * c[k-i] else - mul!(c[k+1], (k-i) * r[i+1], c[k-i+1]) + mul!(c[k], (k-i) * r[i], c[k-i]) end end @inbounds sqr!(r, a, k) - @inbounds c[k+1] = (a[k+1] - c[k+1]/k) / constant_term(r) + @inbounds c[k] = (a[k] - c[k]/k) / constant_term(r) return nothing end - @inline function sinhcosh!(s::$T, c::$T, a::$T, k::Int) + @inline function sinhcosh!(s::$T{T}, c::$T{T}, a::$T{T}, k::Int) where {T} if k == 0 - @inbounds s[1] = sinh( constant_term(a) ) - @inbounds c[1] = cosh( constant_term(a) ) + @inbounds s[0] = sinh( constant_term(a) ) + @inbounds c[0] = cosh( constant_term(a) ) return nothing end @inbounds for i = 1:k - x = i * a[i+1] + x = i * a[i] if $T == Taylor1 - s[k+1] += x * c[k-i+1] - c[k+1] += x * s[k-i+1] + s[k] += x * c[k-i] + c[k] += x * s[k-i] else - mul!(s[k+1], x, c[k-i+1]) - mul!(c[k+1], x, s[k-i+1]) + mul!(s[k], x, c[k-i]) + mul!(c[k], x, s[k-i]) end end - s[k+1] = s[k+1] / k - c[k+1] = c[k+1] / k + s[k] = s[k] / k + c[k] = c[k] / k return nothing end - @inline function tanh!(c::$T, a::$T, c2::$T, k::Int) + @inline function tanh!(c::$T{T}, a::$T{T}, c2::$T{T}, k::Int) where {T} if k == 0 @inbounds aux = tanh( constant_term(a) ) - @inbounds c[1] = aux - @inbounds c2[1] = aux^2 + @inbounds c[0] = aux + @inbounds c2[0] = aux^2 return nothing end @inbounds for i = 0:k-1 if $T == Taylor1 - c[k+1] += (k-i) * a[k-i+1] * c2[i+1] + c[k] += (k-i) * a[k-i] * c2[i] else - mul!(c[k+1], (k-i) * a[k-i+1], c2[i+1]) + mul!(c[k], (k-i) * a[k-i], c2[i]) end end - @inbounds c[k+1] = a[k+1] - c[k+1]/k + @inbounds c[k] = a[k] - c[k]/k sqr!(c2, c, k) return nothing @@ -319,14 +372,14 @@ for T in (:Taylor1, :TaylorN) end -doc""" +@doc doc""" inverse(f) -Return the Taylor expansion of $f^{-1}(t)$, of order `N = f.order`, +Return the Taylor expansion of ``f^{-1}(t)``, of order `N = f.order`, for `f::Taylor1` polynomial if the first coefficient of `f` is zero. Otherwise, an `ArgumentError` is thrown. -The algorithm implements Lagrange inversion at $t=0$ if $f(0)=0$: +The algorithm implements Lagrange inversion at ``t=0`` if ``f(0)=0``: ```math \begin{equation*} f^{-1}(t) = \sum_{n=1}^{N} \frac{t^n}{n!} \left. @@ -335,9 +388,10 @@ f^{-1}(t) = \sum_{n=1}^{N} \frac{t^n}{n!} \left. \end{equation*} ``` -""" -function inverse{T<:Number}(f::Taylor1{T}) - if f[1] != zero(T) +""" inverse + +function inverse(f::Taylor1{T}) where {T<:Number} + if f[0] != zero(T) throw(ArgumentError( """ Evaluation of Taylor1 series at 0 is non-zero. For high accuracy, revert @@ -352,7 +406,7 @@ function inverse{T<:Number}(f::Taylor1{T}) coeffs[1] = zero(S) @inbounds for n = 1:f.order - coeffs[n+1] = zdivfpown[n]/n + coeffs[n+1] = zdivfpown[n-1]/n zdivfpown *= zdivf end Taylor1(coeffs, f.order) @@ -361,7 +415,7 @@ end # Documentation for the recursion relations -doc""" +@doc doc""" exp!(c, a, k) --> nothing Update the `k-th` expansion coefficient `c[k+1]` of `c = exp(a)` @@ -378,7 +432,7 @@ c_k = \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} c_j. """ exp! -doc""" +@doc doc""" log!(c, a, k) --> nothing Update the `k-th` expansion coefficient `c[k+1]` of `c = log(a)` @@ -395,7 +449,7 @@ c_k = \frac{1}{a_0} \big(a_k - \frac{1}{k} \sum_{j=0}^{k-1} j a_{k-j} c_j \big). """ log! -doc""" +@doc doc""" sincos!(s, c, a, k) --> nothing Update the `k-th` expansion coefficients `s[k+1]` and `c[k+1]` @@ -414,7 +468,7 @@ c_k &=& -\frac{1}{k}\sum_{j=0}^{k-1} (k-j) a_{k-j} s_j. """ sincos! -doc""" +@doc doc""" tan!(c, a, p, k::Int) --> nothing Update the `k-th` expansion coefficients `c[k+1]` of `c = tan(a)`, @@ -432,7 +486,7 @@ c_k = a_k + \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} p_j. """ tan! -doc""" +@doc doc""" asin!(c, a, r, k) Update the `k-th` expansion coefficients `c[k+1]` of `c = asin(a)`, @@ -449,7 +503,7 @@ c_k = \frac{1}{ \sqrt{r_0} } """ asin! -doc""" +@doc doc""" acos!(c, a, r, k) Update the `k-th` expansion coefficients `c[k+1]` of `c = acos(a)`, @@ -466,7 +520,7 @@ c_k = - \frac{1}{ r_0 } """ acos! -doc""" +@doc doc""" atan!(c, a, r, k) Update the `k-th` expansion coefficients `c[k+1]` of `c = atan(a)`, @@ -482,7 +536,7 @@ c_k = \frac{1}{r_0}\big(a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j\big). """ atan! -doc""" +@doc doc""" sinhcosh!(s, c, a, k) Update the `k-th` expansion coefficients `s[k+1]` and `c[k+1]` @@ -501,7 +555,7 @@ c_k &=& \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} s_j. """ sinhcosh! -doc""" +@doc doc""" tanh!(c, a, p, k) Update the `k-th` expansion coefficients `c[k+1]` of `c = tanh(a)`, diff --git a/src/hash_tables.jl b/src/hash_tables.jl index 689760ea..3f94f988 100644 --- a/src/hash_tables.jl +++ b/src/hash_tables.jl @@ -2,7 +2,7 @@ # Hash tables for HomogeneousPolynomial and TaylorN -doc""" +""" generate_tables(num_vars, order) Return the hash tables `coeff_table`, `index_table`, `size_table` @@ -12,22 +12,22 @@ and `pos_table`. Internally, these are treated as `const`. coeff_table :: Array{Array{Array{Int64,1},1},1} -The $i+1$-th component contains a vector with the vectors of all the possible -combinations of monomials of a `HomogeneousPolynomial` of order $i$. +The ``i+1``-th component contains a vector with the vectors of all the possible +combinations of monomials of a `HomogeneousPolynomial` of order ``i``. index_table :: Array{Array{Int64,1},1} -The $i+1$-th component contains a vector of (hashed) indices that represent -the distinct monomials of a `HomogeneousPolynomial` of order (degree) $i$. +The ``i+1``-th component contains a vector of (hashed) indices that represent +the distinct monomials of a `HomogeneousPolynomial` of order (degree) ``i``. size_table :: Array{Int64,1} -The $i+1$-th component contains the number of distinct monomials of the -`HomogeneousPolynomial` of order $i$, equivalent to `length(coeff_table[i])`. +The ``i+1``-th component contains the number of distinct monomials of the +`HomogeneousPolynomial` of order ``i``, equivalent to `length(coeff_table[i])`. pos_table :: Array{Dict{Int64,Int64},1} -The $i+1$-th component maps the hash index to the (lexicographic) position +The ``i+1``-th component maps the hash index to the (lexicographic) position of the corresponding monomial in `coeffs_table`. """ function generate_tables(num_vars, order) @@ -103,4 +103,22 @@ const coeff_table, index_table, size_table, pos_table = generate_tables(get_numvars(), get_order()) # Garbage-collect here to free memory -gc(); +@compat GC.gc(); + + +""" + show_monomials(ord::Int) --> nothing + +List the indices and corresponding of a `HomogeneousPolynomial` +of degree `ord`. +""" +function show_monomials(ord::Int) + z = zeros(Int, TaylorSeries.size_table[ord+1]) + for (index, value) in enumerate(TaylorSeries.coeff_table[ord+1]) + z[index] = 1 + pol = HomogeneousPolynomial(z) + println(" $index --> $(homogPol2str(pol)[4:end])") + z[index] = 0 + end + nothing +end diff --git a/src/nested_taylor.jl b/src/nested_taylor.jl new file mode 100644 index 00000000..92677c93 --- /dev/null +++ b/src/nested_taylor.jl @@ -0,0 +1,56 @@ +# This file is part of the TaylorSeries.jl Julia package, MIT license +# +# Luis Benet & David P. Sanders +# UNAM +# +# MIT Expat license +# + + +#number of variables in a nested polynomial +function get_nested_numvars{T<:Number}(a::Taylor1{T}) + a_aux = a.coeffs[1] + dim_a = 1 + while !(typeof(a_aux) <: TaylorSeries.NumberNotSeries) + dim_a += 1 + a_aux = a_aux.coeffs[1] + end + return dim_a +end + + +function nest{T<:Number}(a::Taylor1{T},num_vars::Integer,dim_a::Integer) + ord = a.order + + for nest in 1:(num_vars-dim_a) + a = Taylor1([a],ord) + end + + return a +end + +function nest{T<:Number}(a::Taylor1{T},num_vars::Int64) + dim_a = get_nested_numvars(a) + ord = a.order + + for nest in 1:(num_vars-dim_a) + a = Taylor1([a],ord) + end + + return a +end + +#Shortcut to define an array of independent variables of nested taylors +function set_nested_variables(numvars::Integer,order::Integer,T::Type=Float64) + nested = Taylor1([zero(T),one(T)],order) + nested_array = Vector{Number}(numvars) + numvars_nested = 1 + + for i in 1:numvars + nested_array[i] = nest(nested,numvars,numvars_nested) + nested = Taylor1([zero(nested),one(nested)],order) + numvars_nested += 1 + end + + return nested_array +end diff --git a/src/other_functions.jl b/src/other_functions.jl index 6a049f4e..3d0e9019 100644 --- a/src/other_functions.jl +++ b/src/other_functions.jl @@ -26,13 +26,13 @@ end for op in (:mod, :rem) for T in (:Taylor1, :TaylorN) @eval begin - function ($op){T<:Real}(a::$T{T}, x::T) + function ($op)(a::$T{T}, x::T) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = ($op)(constant_term(a), x) return $T(coeffs, a.order) end - function ($op){T<:Real, S<:Real}(a::$T{T}, x::S) + function ($op)(a::$T{T}, x::S) where {T<:Real, S<:Real} R = promote_type(T, S) a = convert($T{R}, a) return ($op)(a, convert(R,x)) @@ -41,25 +41,27 @@ for op in (:mod, :rem) end @eval begin - function ($op){T<:Real}(a::TaylorN{Taylor1{T}}, x::T) + function ($op)(a::TaylorN{Taylor1{T}}, x::T) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = ($op)(constant_term(a), x) return TaylorN( coeffs, a.order ) end - function ($op){T<:Real,S<:Real}(a::TaylorN{Taylor1{T}}, x::S) + function ($op)(a::TaylorN{Taylor1{T}}, x::S) where {T<:Real, S<:Real} R = promote_type(T,S) a = convert(TaylorN{Taylor1{R}}, a) return ($op)(a, convert(R,x)) end - function ($op){T<:Real}(a::Taylor1{TaylorN{T}}, x::T) + function ($op)(a::Taylor1{TaylorN{T}}, x::T) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = ($op)(constant_term(a), x) return Taylor1( coeffs, a.order ) end - @inbounds function ($op){T<:Real,S<:Real}(a::Taylor1{TaylorN{T}}, x::S) + @inbounds function ($op)(a::Taylor1{TaylorN{T}}, x::S) where + {T<:Real, S<:Real} + R = promote_type(T,S) a = convert(Taylor1{TaylorN{R}}, a) return ($op)(a, convert(R,x)) @@ -71,16 +73,16 @@ end ## mod2pi and abs ## for T in (:Taylor1, :TaylorN) @eval begin - function mod2pi{T<:Real}(a::$T{T}) + function mod2pi(a::$T{T}) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = mod2pi( constant_term(a) ) return $T( coeffs, a.order) end - function abs{T<:Real}(a::$T{T}) - if constant_term(a) > zero(T) + function abs(a::$T{T}) where {T<:Real} + if constant_term(a) > 0 return a - elseif constant_term(a) < zero(T) + elseif constant_term(a) < 0 return -a else throw(ArgumentError( @@ -88,25 +90,27 @@ for T in (:Taylor1, :TaylorN) (abs(x) is not differentiable at x=0).""")) end end + + abs2(a::$T) = a^2 end end -function mod2pi{T<:Real}(a::TaylorN{Taylor1{T}}) +function mod2pi(a::TaylorN{Taylor1{T}}) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = mod2pi( constant_term(a) ) return TaylorN( coeffs, a.order ) end -function mod2pi{T<:Real}(a::Taylor1{TaylorN{T}}) +function mod2pi(a::Taylor1{TaylorN{T}}) where {T<:Real} coeffs = copy(a.coeffs) @inbounds coeffs[1] = mod2pi( constant_term(a) ) return Taylor1( coeffs, a.order ) end -function abs{T<:Real}(a::TaylorN{Taylor1{T}}) - if constant_term(a)[1] > zero(T) +function abs(a::TaylorN{Taylor1{T}}) where {T<:Real} + if constant_term(a)[0] > 0 return a - elseif constant_term(a)[1] < zero(T) + elseif constant_term(a)[0] < 0 return -a else throw(ArgumentError( @@ -115,10 +119,10 @@ function abs{T<:Real}(a::TaylorN{Taylor1{T}}) end end -function abs{T<:Real}(a::Taylor1{TaylorN{T}}) - if constant_term(a[1]) > zero(T) +function abs(a::Taylor1{TaylorN{T}}) where {T<:Real} + if constant_term(a[0]) > 0 return a - elseif constant_term(a[1]) < zero(T) + elseif constant_term(a[0]) < 0 return -a else throw(ArgumentError( @@ -126,3 +130,171 @@ function abs{T<:Real}(a::Taylor1{TaylorN{T}}) (abs(x) is not differentiable at x=0).""")) end end + +@doc doc""" + abs(a) + +Returns `a` if `constant_term(a) > 0` and `-a` if `constant_term(a) < 0` for +`a <:Union{Taylor1,TaylorN}`. +Notice that `typeof(abs(a)) <: AbstractSeries`. + +""" abs + + +#norm +@doc doc""" + norm(x::AbstractSeries, p::Real) + +Returns the p-norm of an `x::AbstractSeries`, defined by + +```math +\begin{equation*} +\left\Vert x \right\Vert_p = \left( \sum_k | x_k |^p \right)^{\frac{1}{p}}, +\end{equation*} +``` +which returns a non-negative number. + +""" norm + +norm(x::AbstractSeries, p::Real=2) = norm( norm.(x.coeffs, p), p) +norm(x::Union{Taylor1{T},HomogeneousPolynomial{T}}, p::Real=2) where + {T<:NumberNotSeries} = norm(x.coeffs, p) +#norm for Taylor vectors +norm(v::Vector{T}, p::Real=2) where T<:AbstractSeries = norm( norm.(v, p), p) + +# rtoldefault +for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) + @eval rtoldefault(::Type{$T{T}}) where {T<:Number} = rtoldefault(T) + @eval rtoldefault(::$T{T}) where {T<:Number} = rtoldefault(T) +end + +# isfinite +""" + isfinite(x::AbstractSeries) -> Bool + +Test whether the coefficients of the polynomial `x` are finite. +""" +isfinite(x::AbstractSeries) = !isnan(x) && !isinf(x) + +# isapprox; modified from Julia's Base.isapprox +""" + isapprox(x::AbstractSeries, y::AbstractSeries; rtol::Real=sqrt(eps), atol::Real=0, nans::Bool=false) + +Inexact equality comparison between polynomials: returns `true` if +`norm(x-y,1) <= atol + rtol*max(norm(x,1), norm(y,1))`, where `x` and `y` are +polynomials. For more details, see [`Base.isapprox`](@ref). +""" +@compat function isapprox(x::T, y::S; rtol::Real=rtoldefault(x,y,0), atol::Real=0.0, + nans::Bool=false) where {T<:AbstractSeries, S<:AbstractSeries} + + x == y || (isfinite(x) && isfinite(y) && + norm(x-y,1) <= atol + rtol*max(norm(x,1), norm(y,1))) || + (nans && isnan(x) && isnan(y)) +end +#isapprox for vectors of Taylors +@compat function isapprox(x::Vector{T}, y::Vector{S}; rtol::Real=rtoldefault(T,S,0), atol::Real=0.0, + nans::Bool=false) where {T<:AbstractSeries, S<:AbstractSeries} + + x == y || norm(x-y,1) <= atol + rtol*max(norm(x,1), norm(y,1)) || + (nans && isnan(x) && isnan(y)) +end + +#taylor_expand function for Taylor1 +""" + taylor_expand(f, x0; order) + +Computes the Taylor expansion of the function `f` around the point `x0`. + +If `x0` is a scalar, a `Taylor1` expansion will be returned. If `x0` is a vector, +a `TaylorN` expansion will be computed. If the dimension of x0 (`length(x0)`) +is different from the variables set for `TaylorN` (`get_numvars()`), an +`AssertionError` will be thrown. +""" +function taylor_expand(f::Function; order::Int64=15) + a = Taylor1(order) + return f(a) +end + +function taylor_expand(f::Function, x0::T; order::Int64=15) where {T<:Number} + a = Taylor1([x0, one(T)], order) + return f(a) +end + +#taylor_expand function for TaylorN +function taylor_expand(f::Function, x0::Vector{T}; + order::Int64=get_order()) where {T<:Number} + + ll = length(x0) + @assert ll == get_numvars() && order <= get_order() + @compat X = Array{TaylorN{T}}(undef, ll) + + for i in eachindex(X) + X[i] = x0[i] + TaylorN(T, i, order=order) + end + + return f( X ) +end + +function taylor_expand(f::Function, x0...; order::Int64=get_order()) + x0 = promote(x0...) + T = eltype(x0[1]) + ll = length(x0) + @assert ll == get_numvars() && order <= get_order() + @compat X = Array{TaylorN{T}}(undef, ll) + + for i in eachindex(X) + X[i] = x0[i] + TaylorN(T, i, order=order) + end + + return f( X... ) +end + +#update! function for Taylor1 +""" + update!(a, x0) + +Takes `a <: Union{Taylo1,TaylorN}` and expands it around the coordinate `x0`. +""" +function update!(a::Taylor1, x0::T) where {T<:Number} + a.coeffs .= evaluate(a, Taylor1([x0, one(x0)], a.order) ).coeffs + nothing +end + +#update! function for TaylorN +function update!(a::TaylorN, vals::Vector{T}) where {T<:Number} + a.coeffs .= evaluate(a, get_variables() .+ vals).coeffs + nothing +end + +function update!(a::Union{Taylor1,TaylorN}) + #shifting around zero shouldn't change anything... + nothing +end + +for T in (:Taylor1, :TaylorN) + @eval deg2rad(z::$T{T}) where {T<:AbstractFloat} = z * (convert(T, pi) / 180) + @eval deg2rad(z::$T{T}) where {T<:Real} = z * (convert(float(T), pi) / 180) + + @eval rad2deg(z::$T{T}) where {T<:AbstractFloat} = z * (180 / convert(T, pi)) + @eval rad2deg(z::$T{T}) where {T<:Real} = z * (180 / convert(float(T), pi)) +end + +# 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} + @inbounds v[k] = a[k] * (convert(T, pi) / 180) + return nothing + end + @eval @inline function deg2rad!(v::$T{S}, a::$T{T}, k::Int) where {S<:AbstractFloat, T<:Real} + @inbounds v[k] = a[k] * (convert(float(T), pi) / 180) + return nothing + end + @eval @inline function rad2deg!(v::$T{T}, a::$T{T}, k::Int) where {T<:AbstractFloat} + @inbounds v[k] = a[k] * (180 / convert(T, pi)) + return nothing + end + @eval @inline function rad2deg!(v::$T{S}, a::$T{T}, k::Int) where {S<:AbstractFloat, T<:Real} + @inbounds v[k] = a[k] * (180 / convert(float(T), pi)) + return nothing + end +end diff --git a/src/parameters.jl b/src/parameters.jl index fa515f99..7d716c9a 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -3,7 +3,7 @@ # Parameters for HomogeneousPolynomial and TaylorN -doc""" +""" ParamsTaylorN DataType holding the current parameters for `TaylorN` and @@ -11,29 +11,52 @@ DataType holding the current parameters for `TaylorN` and **Fields:** -- `order :: Int` Order (degree) of the polynomials -- `num_vars :: Int` Number of variables -- `variable_names :: Array{String,1}` Name of the variables +- `order :: Int` Order (degree) of the polynomials +- `num_vars :: Int` Number of variables +- `variable_names :: Vector{String}` Names of the variables +- `variable_symbols :: Vector{Symbol}` Symbols of the variables -These parameters can be changed using `set_params_TaylorN(order, numVars)`. +These parameters can be changed using [`set_variables`](@ref) """ -type ParamsTaylorN - order :: Int - num_vars :: Int - variable_names :: Array{String,1} +mutable struct ParamsTaylorN + order :: Int + num_vars :: Int + variable_names :: Vector{String} + variable_symbols :: Vector{Symbol} end -const _params_TaylorN_ = ParamsTaylorN(6, 2, String["x₁", "x₂"]) + +ParamsTaylorN(order, num_vars, variable_names) = ParamsTaylorN(order, num_vars, variable_names, Symbol.(variable_names)) + +const _params_TaylorN_ = ParamsTaylorN(6, 2, ["x₁", "x₂"]) -## Utilities to get the maximum order and number of variables +## Utilities to get the maximum order, number of variables, their names and symbols get_order() = _params_TaylorN_.order get_numvars() = _params_TaylorN_.num_vars get_variable_names() = _params_TaylorN_.variable_names +get_variable_symbols() = _params_TaylorN_.variable_symbols +function lookupvar(s::Symbol) + @compat ind = findfirst(x -> x==s, _params_TaylorN_.variable_symbols) + @compat isa(ind, Nothing) && return 0 + return ind +end -set_variable_names{T<:AbstractString}(names::Vector{T}) = _params_TaylorN_.variable_names = names +function set_variable_names(varnames::Vector{T}) where {T<:AbstractString} + _params_TaylorN_.variable_names = varnames + _params_TaylorN_.variable_symbols = Symbol.(varnames) + nothing +end +""" + get_variables(;order=get_order()) -get_variables() = [TaylorN(i) for i in 1:get_numvars()] +Return a `TaylorN` vector with each entry representing an +independent variable. It takes the default `_params_TaylorN_` values +if `set_variables` hasn't been changed with the exception that `order` +can be explicitely established by the user without changing internal values +for `num_vars` or `variable_names`. +""" +get_variables(;order::Integer=get_order()) = [TaylorN(i,order=order) for i in 1:get_numvars()] """ set_variables([T::Type], names::String; [order=get_order(), numvars=-1]) @@ -66,13 +89,17 @@ julia> set_variables("x", order=6, numvars=2) 1.0 x₂ + 𝒪(‖x‖⁷) ``` """ -function set_variables{T<:AbstractString}(R::Type, names::Vector{T}; order=get_order()) +function set_variables(R::Type, names::Vector{T}; order=get_order()) where + {T<:AbstractString} + order ≥ 1 || error("Order must be at least 1") num_vars = length(names) num_vars ≥ 1 || error("Number of variables must be at least 1") _params_TaylorN_.variable_names = names + _params_TaylorN_.variable_symbols = Symbol.(names) + if !(order == get_order() && num_vars == get_numvars()) # if these are unchanged, no need to regenerate tables @@ -85,18 +112,25 @@ function set_variables{T<:AbstractString}(R::Type, names::Vector{T}; order=get_o resize!(size_table,order+1) resize!(pos_table,order+1) - coeff_table[:], index_table[:], size_table[:], pos_table[:] = generate_tables(num_vars, order) - gc(); + coeff_table[:], index_table[:], size_table[:], pos_table[:] = + generate_tables(num_vars, order) + @compat GC.gc(); end # return a list of the new variables TaylorN{R}[TaylorN(R,i) for i in 1:get_numvars()] end +set_variables(R::Type, symbs::Vector{T}; order=get_order()) where + {T<:Symbol} = set_variables(R, string.(symbs), order=order) -set_variables{T}(names::Vector{T}; order=get_order()) = +set_variables(names::Vector{T}; order=get_order()) where {T<:AbstractString} = set_variables(Float64, names, order=order) +set_variables(symbs::Vector{T}; order=get_order()) where {T<:Symbol} = + set_variables(Float64, symbs, order=order) + +function set_variables(R::Type, names::T; order=get_order(), numvars=-1) where + {T<:AbstractString} -function set_variables{T<:AbstractString}(R::Type, names::T; order=get_order(), numvars=-1) variable_names = split(names) if length(variable_names) == 1 && numvars ≥ 1 @@ -106,9 +140,13 @@ function set_variables{T<:AbstractString}(R::Type, names::T; order=get_order(), set_variables(R, variable_names, order=order) end +set_variables(R::Type, symbs::Symbol; order=get_order(), numvars=-1) = + set_variables(R, string(symbs), order=order, numvars=numvars) -set_variables{T<:AbstractString}(names::T; order=get_order(), numvars=-1) = +set_variables(names::T; order=get_order(), numvars=-1) where {T<:AbstractString} = set_variables(Float64, names, order=order, numvars=numvars) +set_variables(symbs::Symbol; order=get_order(), numvars=-1) = + set_variables(Float64, string(symbs), order=order, numvars=numvars) """ @@ -117,11 +155,25 @@ set_variables{T<:AbstractString}(names::T; order=get_order(), numvars=-1) = Display the current parameters for `TaylorN` and `HomogeneousPolynomial` types. """ function show_params_TaylorN() - info( """ + Compat.@info( """ Parameters for `TaylorN` and `HomogeneousPolynomial`: Maximum order = $(get_order()) Number of variables = $(get_numvars()) Variable names = $(get_variable_names()) + Variable symbols = $(Symbol.(get_variable_names())) """) nothing end + + +# Control the display of the big 𝒪 notation +const bigOnotation = Bool[true] + +""" + displayBigO(d::Bool) --> nothing + +Set/unset displaying of the big 𝒪 notation in the output +of `Taylor1` and `TaylorN` polynomials. The initial value is +`true`. +""" +displayBigO(d::Bool) = (bigOnotation[end] = d; d) diff --git a/src/power.jl b/src/power.jl index 98bc6bb2..cbfed99e 100644 --- a/src/power.jl +++ b/src/power.jl @@ -15,20 +15,58 @@ function ^(a::HomogeneousPolynomial, n::Integer) return power_by_squaring(a, n) end +#= The following three methods are coded like that, to use +preferentially a^float(n), but for cases like Taylor1{Interval{T}}^n +power_by_squaring is used. The latter is important when the +0-th order coefficient is/contains zero. +=# +function ^(a::Taylor1{T}, n::Integer) where {T<:Number} + n == 0 && return one(a) + n == 1 && return copy(a) + n == 2 && return square(a) + return a^float(n) +end + +# Method used for Taylor1{Interval{T}}^n +function ^(a::Taylor1{T}, n::Integer) where {T<:Real} + n == 0 && return one(a) + n == 1 && return copy(a) + n == 2 && return square(a) + n < 0 && return a^float(n) + return power_by_squaring(a, n) +end + +function ^(a::Taylor1{T}, n::Integer) where {T<:AbstractFloat} + n == 0 && return one(a) + n == 1 && return copy(a) + n == 2 && return square(a) + return a^float(n) +end + + +function ^(a::TaylorN{T}, n::Integer) where {T<:Number} + n == 0 && return one(a) + n == 1 && return copy(a) + n == 2 && return square(a) + n < 0 && return inv( a^(-n) ) + return power_by_squaring(a, n) +end + + for T in (:Taylor1, :TaylorN) - @eval function ^{T<:Number}(a::$T{T}, n::Integer) + @eval function ^(a::$T{T}, n::Integer) where {T<:Integer} n == 0 && return one(a) n == 1 && return copy(a) n == 2 && return square(a) - n < 0 && return inv( a^(-n) ) + n < 0 && throw(DomainError()) return power_by_squaring(a, n) end - @eval function ^{T<:Integer}(a::$T{T}, n::Integer) + @eval function ^(a::$T{Rational{T}}, n::Integer) where {T<:Integer} n == 0 && return one(a) n == 1 && return copy(a) n == 2 && return square(a) - n < 0 && throw(DomainError()) + n < 0 && return inv( a^(-n) ) return power_by_squaring(a, n) end @@ -36,7 +74,7 @@ for T in (:Taylor1, :TaylorN) @eval ^(a::$T, b::$T) = exp( b*log(a) ) - @eval ^{T<:Complex}(a::$T, x::T) = exp( x*log(a) ) + @eval ^(a::$T, x::T) where {T<:Complex} = exp( x*log(a) ) end @@ -69,44 +107,57 @@ for T in (:Taylor1, :HomogeneousPolynomial, :TaylorN) end ## Real power ## -function ^{S<:Real}(a::Taylor1, r::S) - r == zero(r) && return one(a) - r == one(r)/2 && return sqrt(a) - isinteger(r) && return a^round(Int,r) +function ^(a::Taylor1, r::S) where {S<:Real} + a0 = constant_term(a) + aux = one(a0^r) + + iszero(r) && return Taylor1(aux, a.order) + aa = aux*a + r == 1 && return aa + r == 2 && return square(aa) + r == 1/2 && return sqrt(aa) # First non-zero coefficient l0nz = findfirst(a) - l0nz < 0 && return zero(a) + l0nz < 0 && return Taylor1(zero(aux), a.order) # The first non-zero coefficient of the result; must be integer !isinteger(r*l0nz) && throw(ArgumentError( """The 0th order Taylor1 coefficient must be non-zero to raise the Taylor1 polynomial to a non-integer exponent.""")) lnull = trunc(Int, r*l0nz ) + lnull > a.order && return Taylor1(zero(aux), a.order) - @inbounds aux = ( a[l0nz+1] )^r k0 = lnull+l0nz c = Taylor1( zero(aux), a.order) - @inbounds c[lnull+1] = aux - for k = k0+1:a.order - pow!(c, a, r, k, l0nz) + @inbounds c[lnull] = ( a[l0nz] )^r + for k = k0+1:a.order+l0nz + pow!(c, aa, r, k, l0nz) end return c end ## Real power ## -function ^{S<:Real}(a::TaylorN, r::S) - r == zero(r) && return TaylorN( one(eltype(a)), 0 ) - r == one(r)/2 && return sqrt(a) - isinteger(r) && return a^round(Int,r) - +function ^(a::TaylorN, r::S) where {S<:Real} a0 = constant_term(a) - @assert a0 != zero(a0) + aux = one(a0^r) + + iszero(r) && return TaylorN(aux, a.order) + aa = aux*a + r == 1 && return aa + r == 2 && return square(aa) + r == one(r)/2 && return sqrt(aa) + isinteger(r) && return aa^round(Int,r) + + # @assert !iszero(a0) + iszero(a0) && throw(ArgumentError( + """The 0-th order TaylorN coefficient must be non-zero + in order to expand `^` around 0.""")) c = TaylorN( a0^r, a.order) for ord in 1:a.order - pow!(c, a, r, ord) + pow!(c, aa, r, ord) end return c @@ -114,10 +165,10 @@ end # Homogeneous coefficients for real power -doc""" +@doc doc""" pow!(c, a, r::Real, k::Int, k0::Int=0) -Update the `k-th` expansion coefficient `c[k+1]` of `c = a^r`, for +Update the `k`-th expansion coefficient `c[k]` of `c = a^r`, for both `c` and `a` either `Taylor1` or `TaylorN`. The coefficients are given by @@ -129,34 +180,67 @@ c_k = \frac{1}{k a_0} \sum_{j=0}^{k-1} \big(r(k-j) -j\big)a_{k-j} c_j. For `Taylor1` polynomials, `k0` is the order of the first non-zero coefficient of `a`. -""" -@inline function pow!{S<:Real}(c::Taylor1, a::Taylor1, r::S, k::Int, l0::Int=0) +""" pow! + +@inline function pow!(c::Taylor1{T}, a::Taylor1{T}, r::S, k::Int, l0::Int=0) where + {T<:Number, S<:Real} + + if r == 0 + return one!(c, a, k) + elseif r == 1 + return identity!(c, a, k) + elseif r == 2 + return sqr!(c, a, k) + elseif r == 0.5 + return sqrt!(c, a, k) + end + if k == l0 - @inbounds c[1] = ( a[l0+1] )^r + @inbounds c[0] = ( a[l0] )^r + return nothing + end + + # Relevant for positive integer r, to avoid round-off errors + if isinteger(r) && (k-l0 > r*findlast(a) > 0) + @inbounds c[k-l0] = zero(c[0]) return nothing end - for i = 0:k-l0-1 + imin = max(0,k-a.order) + for i = imin:k-l0-1 + # k-i > a.order && continue aux = r*(k-i) - i - @inbounds c[k-l0+1] += aux * a[k-i+1] * c[i+1] + @inbounds c[k-l0] += aux * a[k-i] * c[i] end aux = k - l0*(r+1) - @inbounds c[k-l0+1] = c[k-l0+1] / (aux * a[l0+1]) + @inbounds c[k-l0] = c[k-l0] / (aux * a[l0]) return nothing end -@inline function pow!{S<:Real}(c::TaylorN, a::TaylorN, r::S, k::Int) +@inline function pow!(c::TaylorN{T}, a::TaylorN{T}, r::S, k::Int) where + {T<:NumberNotSeriesN, S<:Real} + + if r == 0 + return one!(c, a, k) + elseif r == 1 + return identity!(c, a, k) + elseif r == 2 + return sqr!(c, a, k) + elseif r == 0.5 + return sqrt!(c, a, k) + end + if k == 0 - @inbounds c[1] = ( constant_term(a) )^r + @inbounds c[0] = ( constant_term(a) )^r return nothing end for i = 0:k-1 aux = r*(k-i) - i - mul!(c[k+1], aux*a[k-i+1], c[i+1]) + mul!(c[k], aux*a[k-i], c[i]) end - @inbounds c[k+1] = c[k+1] / (k * constant_term(a)) + @inbounds c[k] = c[k] / (k * constant_term(a)) return nothing end @@ -193,10 +277,10 @@ end # Homogeneous coefficients for square -doc""" +@doc doc""" sqr!(c, a, k::Int) --> nothing -Update the `k-th` expansion coefficient `c[k+1]` of `c = a^2`, for +Update the `k-th` expansion coefficient `c[k]` of `c = a^2`, for both `c` and `a` either `Taylor1` or `TaylorN`. The coefficients are given by @@ -214,9 +298,9 @@ c_k & = & 2 \sum_{j=0}^{(k-2)/2} a_{k-j} a_j + (a_{k/2})^2, for T = (:Taylor1, :TaylorN) @eval begin - @inline function sqr!(c::$T, a::$T, k::Int) + @inline function sqr!(c::$T{T}, a::$T{T}, k::Int) where {T} if k == 0 - @inbounds c[1] = constant_term(a)^2 + @inbounds c[0] = constant_term(a)^2 return nothing end @@ -224,18 +308,18 @@ for T = (:Taylor1, :TaylorN) kend = div(k - 2 + kodd, 2) @inbounds for i = 0:kend if $T == Taylor1 - c[k+1] += a[i+1] * a[k-i+1] + c[k] += a[i] * a[k-i] else - mul!(c[k+1], a[i+1], a[k-i+1]) + mul!(c[k], a[i], a[k-i]) end end - @inbounds c[k+1] = 2 * c[k+1] + @inbounds c[k] = 2 * c[k] kodd == 1 && return nothing if $T == Taylor1 - @inbounds c[k+1] += a[div(k,2)+1]^2 + @inbounds c[k] += a[div(k,2)]^2 else - sqr!(c[k+1], a[div(k,2)+1]) + sqr!(c[k], a[div(k,2)]) end return nothing @@ -250,10 +334,10 @@ end Return `c = a*a` with no allocation; all parameters are `HomogeneousPolynomial`. """ -@inline function sqr!(c::HomogeneousPolynomial, a::HomogeneousPolynomial) +@inline function sqr!(c::HomogeneousPolynomial{T}, a::HomogeneousPolynomial{T}) where + {T<:NumberNotSeriesN} iszero(a) && return nothing - T = eltype(c) @inbounds num_coeffs_a = size_table[a.order+1] @inbounds posTb = pos_table[c.order+1] @@ -261,13 +345,13 @@ Return `c = a*a` with no allocation; all parameters are `HomogeneousPolynomial`. @inbounds for na = 1:num_coeffs_a ca = a[na] - ca == zero(T) && continue + iszero(ca) && continue inda = idxTb[na] pos = posTb[2*inda] - c[pos] += ca * ca + c[pos] += ca^2 @inbounds for nb = na+1:num_coeffs_a cb = a[nb] - cb == zero(T) && continue + iszero(cb) && continue indb = idxTb[nb] pos = posTb[inda+indb] c[pos] += 2 * ca * cb @@ -284,8 +368,9 @@ function sqrt(a::Taylor1) # First non-zero coefficient l0nz = findfirst(a) + aux = zero(sqrt( constant_term(a) )) if l0nz < 0 - return zero(a) + return Taylor1(aux, a.order) elseif l0nz%2 == 1 # l0nz must be pair throw(ArgumentError( """First non-vanishing Taylor1 coefficient must correspond @@ -294,13 +379,11 @@ function sqrt(a::Taylor1) # The last l0nz coefficients are set to zero. lnull = div(l0nz, 2) - @inbounds aux = sqrt( a[l0nz+1] ) - T = typeof(aux) - - c = Taylor1( zeros(T, a.order+1) ) - @inbounds c[lnull+1] = aux - for k = lnull+1:a.order-l0nz - sqrt!(c, a, k, lnull) + c = Taylor1( aux, a.order ) + @inbounds c[lnull] = sqrt( a[l0nz] ) + aa = one(aux) * a + for k = lnull+1:a.order + sqrt!(c, aa, k, lnull) end return c @@ -308,80 +391,87 @@ end function sqrt(a::TaylorN) @inbounds p0 = sqrt( constant_term(a) ) - if p0 == zero(p0) + if iszero(p0) throw(ArgumentError( """The 0-th order TaylorN coefficient must be non-zero in order to expand `sqrt` around 0.""")) end c = TaylorN( p0, a.order) + aa = one(p0)*a for k = 1:a.order - sqrt!(c, a, k) + sqrt!(c, aa, k) end return c end # Homogeneous coefficients for the square-root -doc""" +@doc doc""" sqrt!(c, a, k::Int, k0::Int=0) -Compute the `k-th` expansion coefficient `c[k+1]` of `c = sqrt(a)` +Compute the `k-th` expansion coefficient `c[k]` of `c = sqrt(a)` for both`c` and `a` either `Taylor1` or `TaylorN`. The coefficients are given by ```math \begin{eqnarray*} -c_k &=& \frac{1}{2 c_0} \big( a_k - 2 \sum_{j=0}^{(k-1)/2} c_{k-j}c_j\big), - \text{ if $k$ is odd,} \\ -c_k &=& \frac{1}{2 c_0} \big( a_k - 2 \sum_{j=0}^{(k-2)/2} c_{k-j}c_j\big) - - (c_{k/2})^2, \text{ if $k$ is even.} +c_k &=& \frac{1}{2 c_0} \big( a_k - 2 \sum_{j=1}^{(k-1)/2} c_{k-j}c_j\big), + \text{ if k is odd,} \\ +c_k &=& \frac{1}{2 c_0} \big( a_k - 2 \sum_{j=1}^{(k-2)/2} c_{k-j}c_j + - (c_{k/2})^2\big), \text{ if k is even.} \end{eqnarray*} ``` For `Taylor1` polynomials, `k0` is the order of the first non-zero coefficient, which must be even. -""" -@inline function sqrt!(c::Taylor1, a::Taylor1, k::Int, k0::Int=0) +""" sqrt! + +@inline function sqrt!(c::Taylor1{T}, a::Taylor1{T}, k::Int, k0::Int=0) where {T} if k == k0 - @inbounds c[k+1] = sqrt(a[2*k0+1]) + @inbounds c[k] = sqrt(a[2*k0]) return nothing end kodd = (k - k0)%2 kend = div(k - k0 - 2 + kodd, 2) @inbounds for i = k0+1:k0+kend - c[k+1] += c[i+1] * c[k+k0-i+1] + (k+k0-i > a.order) || (i > a.order) && continue + c[k] += c[i] * c[k+k0-i] + end + if k+k0 ≤ a.order + @inbounds aux = a[k+k0] - 2*c[k] + else + @inbounds aux = - 2*c[k] end - @inbounds aux = a[k+k0+1] - 2*c[k+1] if kodd == 0 - @inbounds aux = aux - (c[kend+k0+2])^2 + @inbounds aux = aux - (c[kend+k0+1])^2 end - @inbounds c[k+1] = aux / (2*c[k0+1]) + @inbounds c[k] = aux / (2*c[k0]) return nothing end -@inline function sqrt!(c::TaylorN, a::TaylorN, k::Int) +@inline function sqrt!(c::TaylorN{T}, a::TaylorN{T}, k::Int) where {T<:NumberNotSeriesN} if k == 0 - @inbounds c[1] = sqrt( constant_term(a) ) + @inbounds c[0] = sqrt( constant_term(a) ) return nothing end kodd = k%2 kend = div(k - 2 + kodd, 2) @inbounds for i = 1:kend - mul!(c[k+1], c[i+1], c[k-i+1]) + mul!(c[k], c[i], c[k-i]) end - @inbounds aux = a[k+1] - 2*c[k+1] + @inbounds aux = a[k] - 2*c[k] if kodd == 0 - @inbounds aux = aux - (c[kend+2])^2 + @inbounds aux = aux - (c[kend+1])^2 end - @inbounds c[k+1] = aux / (2*constant_term(c)) + @inbounds c[k] = aux / (2*constant_term(c)) return nothing end diff --git a/src/printing.jl b/src/printing.jl index 60d5dbc1..8ea50ab9 100644 --- a/src/printing.jl +++ b/src/printing.jl @@ -19,17 +19,20 @@ function superscriptify(n::Int) end -function pretty_print{T<:Number}(a::Taylor1{T}) - z = zero(a[1]) +function pretty_print(a::Taylor1{T}) where {T<:Number} + z = zero(a[0]) space = string(" ") - bigO = string("+ 𝒪(t", superscriptify(a.order+1), ")") - a == zero(a) && return string(space, z, space, bigO) + # bigO = string("+ 𝒪(t", superscriptify(a.order+1), ")") + bigO = bigOnotation[end] ? + string("+ 𝒪(t", superscriptify(a.order+1), ")") : + string("") + iszero(a) && return string(space, z, space, bigO) strout::String = space ifirst = true for i in eachindex(a.coeffs) monom::String = i==1 ? string("") : i==2 ? string(" t") : string(" t", superscriptify(i-1)) - @inbounds c = a[i] + @inbounds c = a[i-1] c == z && continue cadena = numbr2str(c, ifirst) strout = string(strout, cadena, monom, space) @@ -38,17 +41,21 @@ function pretty_print{T<:Number}(a::Taylor1{T}) strout = strout * bigO strout end -function pretty_print{T<:Number}(a::Taylor1{HomogeneousPolynomial{T}}) - z = zero(a[1]) + +function pretty_print(a::Taylor1{HomogeneousPolynomial{T}}) where {T<:Number} + z = zero(a[0]) space = string(" ") - bigO = string("+ 𝒪(t", superscriptify(a.order+1), ")") - a == zero(a) && return string(space, z, space, bigO) + # bigO = string("+ 𝒪(t", superscriptify(a.order+1), ")") + bigO = bigOnotation[end] ? + string("+ 𝒪(t", superscriptify(a.order+1), ")") : + string("") + iszero(a) && return string(space, z, space, bigO) strout::String = space ifirst = true for i in eachindex(a.coeffs) monom::String = i==1 ? string("") : i==2 ? string(" t") : string(" t", superscriptify(i-1)) - @inbounds c = a[i] + @inbounds c = a[i-1] c == z && continue cadena = numbr2str(c, ifirst) ccad::String = i==1 ? cadena : ifirst ? string("(", cadena, ")") : @@ -59,17 +66,22 @@ function pretty_print{T<:Number}(a::Taylor1{HomogeneousPolynomial{T}}) strout = strout * bigO strout end -function pretty_print{T<:Number}(a::Taylor1{TaylorN{T}}) - z = zero(a[1]) + +function pretty_print(a::Taylor1{TaylorN{T}}) where {T<:Number} + z = zero(a[0]) space = string(" ") - bigO = string("+ 𝒪(t", superscriptify(a.order+1), ")") - a == zero(a) && return string(space, z, space, bigO) + # bigO = string("+ 𝒪(t", superscriptify(a.order+1), ")") + bigO = bigOnotation[end] ? + string("+ 𝒪(t", superscriptify(a.order+1), ")") : + string("") + iszero(a) && return string(space, z, space, bigO) + iszero(a) && return string(space, z, space, bigO) strout::String = space ifirst = true for i in eachindex(a.coeffs) monom::String = i==1 ? string("") : i==2 ? string(" t") : string(" t", superscriptify(i-1)) - @inbounds c = a[i] + @inbounds c = a[i-1] c == z && continue cadena = numbr2str(c, ifirst) ccad::String = i==1 ? cadena : ifirst ? string("(", cadena, ")") : @@ -80,23 +92,29 @@ function pretty_print{T<:Number}(a::Taylor1{TaylorN{T}}) strout = strout * bigO strout end -function pretty_print{T<:Number}(a::HomogeneousPolynomial{T}) + +function pretty_print(a::HomogeneousPolynomial{T}) where {T<:Number} z = zero(a[1]) space = string(" ") - a == zero(a) && return string(space, z) + iszero(a) && return string(space, z) strout::String = homogPol2str(a) strout end -function pretty_print{T<:Number}(a::TaylorN{T}) - z = zero(a[1]) + +function pretty_print(a::TaylorN{T}) where {T<:Number} + z = zero(a[0]) space = string("") - bigO::String = string(" + 𝒪(‖x‖", superscriptify(a.order+1), ")") - a == zero(a) && return string(space, z, bigO) + # bigO::String = string(" + 𝒪(‖x‖", superscriptify(a.order+1), ")") + bigO::String = bigOnotation[end] ? + string(" + 𝒪(‖x‖", superscriptify(a.order+1), ")") : + string("") + iszero(a) && return string(space, z, space, bigO) + iszero(a) && return string(space, z, bigO) strout::String = space#string("") ifirst = true for ord in eachindex(a.coeffs) - pol = a[ord] - pol == zero(a[ord]) && continue + pol = a[ord-1] + iszero(pol) && continue cadena::String = homogPol2str( pol ) strsgn = (ifirst || ord == 1 || cadena[2] == '-') ? string("") : string(" +") @@ -107,7 +125,7 @@ function pretty_print{T<:Number}(a::TaylorN{T}) strout end -function homogPol2str{T<:Number}(a::HomogeneousPolynomial{T}) +function homogPol2str(a::HomogeneousPolynomial{T}) where {T<:Number} numVars = get_numvars() order = a.order z = zero(T) @@ -132,9 +150,10 @@ function homogPol2str{T<:Number}(a::HomogeneousPolynomial{T}) strout = string(strout, cadena, monom, space) ifirst = false end - return strout[1:end-1] + return strout[1:prevind(strout, end)] end -function homogPol2str{T<:Number}(a::HomogeneousPolynomial{Taylor1{T}}) + +function homogPol2str(a::HomogeneousPolynomial{Taylor1{T}}) where {T<:Number} numVars = get_numvars() order = a.order z = zero(a[1]) @@ -162,23 +181,26 @@ function homogPol2str{T<:Number}(a::HomogeneousPolynomial{Taylor1{T}}) strout = string(strout, ccad, monom, space) ifirst = false end - return strout[1:end-1] + return strout[1:prevind(strout, end)] end function numbr2str(zz, ifirst::Bool=false) - zz == zero(zz) && return string( zz ) + iszero(zz) && return string( zz ) plusmin = ifelse( ifirst, string(""), string("+ ") ) return string(plusmin, zz) end -function numbr2str{T<:Union{AbstractFloat,Integer,Rational}}(zz::T, ifirst::Bool=false) - zz == zero(T) && return string( zz ) + +function numbr2str(zz::T, ifirst::Bool=false) where + {T<:Union{AbstractFloat,Integer,Rational}} + iszero(zz) && return string( zz ) plusmin = ifelse( zz < zero(T), string("- "), ifelse( ifirst, string(""), string("+ ")) ) return string(plusmin, abs(zz)) end -function numbr2str{T<:Complex}(zz::T, ifirst::Bool=false) + +function numbr2str(zz::T, ifirst::Bool=false) where {T<:Complex} zT = zero(zz.re) - zz == zero(zz) && return string(zT) + iszero(zz) && return string(zT) zre, zim = reim(zz) cadena = string("") if zre > zT @@ -219,10 +241,11 @@ end name_taylorNvar(i::Int) = string(" ", get_variable_names()[i]) - # summary -summary{T<:Number}(a::Taylor1{T}) = string(a.order, "-order ", typeof(a), ":") -function summary{T<:Number}(a::Union{HomogeneousPolynomial{T}, TaylorN{T}}) +summary(a::Taylor1{T}) where {T<:Number} = + string(a.order, "-order ", typeof(a), ":") + +function summary(a::Union{HomogeneousPolynomial{T}, TaylorN{T}}) where {T<:Number} string(a.order, "-order ", typeof(a), " in ", get_numvars(), " variables:") end diff --git a/test/REQUIRE b/test/REQUIRE new file mode 100644 index 00000000..ff005721 --- /dev/null +++ b/test/REQUIRE @@ -0,0 +1 @@ +IntervalArithmetic 0.13.0 diff --git a/test/fateman40.jl b/test/fateman40.jl index 1fe95083..addf2b7d 100644 --- a/test/fateman40.jl +++ b/test/fateman40.jl @@ -2,7 +2,11 @@ # using TaylorSeries -using Base.Test +if VERSION < v"0.7.0-DEV.2004" + using Base.Test +else + using Test +end @testset "Test inspired by Fateman (takes a few seconds)" begin x, y, z, w = set_variables(Int128, "x", numvars=4, order=40) @@ -27,7 +31,7 @@ using Base.Test f2 = fateman2(20) f3 = fateman3(20) - c = get_coeff(f2,[1,6,7,20]) + c = getcoeff(f2,[1,6,7,20]) @test c == 128358585324486316800 - @test get_coeff(f2,[1,6,7,20]) == c + @test getcoeff(f3,[1,6,7,20]) == c end diff --git a/test/identities_Euler.jl b/test/identities_Euler.jl index 74221402..cb127fe2 100644 --- a/test/identities_Euler.jl +++ b/test/identities_Euler.jl @@ -2,7 +2,11 @@ # using TaylorSeries -using Base.Test +if VERSION < v"0.7.0-DEV.2004" + using Base.Test +else + using Test +end @testset "Testing an identity proved by Euler (8 variables)" begin make_variable(name, index::Int) = string(name, TaylorSeries.subscriptify(index)) diff --git a/test/intervals.jl b/test/intervals.jl new file mode 100644 index 00000000..2b715799 --- /dev/null +++ b/test/intervals.jl @@ -0,0 +1,41 @@ +# This file is part of TaylorSeries.jl, MIT licensed +# + +using TaylorSeries, IntervalArithmetic + +if VERSION < v"0.7.0-DEV.2004" + using Base.Test + eeuler = Base.e +else + using Test + eeuler = Base.MathConstants.e +end + +@testset "Tests Taylor1 and TaylorN expansions over Intervals" begin + a = 1..2 + b = -1 .. 1 + p4(x, a) = x^4 + 4*a*x^3 + 6*a^2*x^2 + 4*a^3*x + a^4 + p5(x, a) = x^5 + 5*a*x^4 + 10*a^2*x^3 + 10*a^3*x^2 + 5*a^4*x + a^5 + + ti = Taylor1(Interval{Float64}, 10) + x, y = set_variables(Interval{Float64}, "x y") + + @test eltype(ti) == Interval{Float64} + @test eltype(x) == Interval{Float64} + + @test p4(ti,-a) == (ti-a)^4 + @test p5(ti,-a) == (ti-a)^5 + @test p4(ti,-b) == (ti-b)^4 + @test all((p5(ti,-b)).coeffs .⊆ ((ti-b)^5).coeffs) + + + @test p4(x,-y) == (x-y)^4 + @test p5(x,-y) == (x-y)^5 + @test p4(x,-a) == (x-a)^4 + @test p5(x,-a) == (x-a)^5 + @test p4(x,-b) == (x-b)^4 + for ind in eachindex((p5(x,-b)).coeffs) + @test all(((p5(x,-b)).coeffs[ind]).coeffs .⊆ (((x-b)^5).coeffs[ind]).coeffs) + end + +end diff --git a/test/manyvariables.jl b/test/manyvariables.jl index 3ad22790..716f1b1c 100644 --- a/test/manyvariables.jl +++ b/test/manyvariables.jl @@ -2,18 +2,38 @@ # using TaylorSeries -using Base.Test +using Compat + +if VERSION < v"0.7.0-DEV.2004" + using Base.Test + eeuler = Base.e +else + using Test + using LinearAlgebra + eeuler = Base.MathConstants.e +end @testset "Tests for HomogeneousPolynomial and TaylorN" begin @test HomogeneousPolynomial <: AbstractSeries @test HomogeneousPolynomial{Int} <: AbstractSeries{Int} @test TaylorN{Float64} <: AbstractSeries{Float64} + _taylorNparams = TaylorSeries.ParamsTaylorN(6, 2, String["x₁", "x₂"]) + @test _taylorNparams.order == get_order() + @test _taylorNparams.num_vars == get_numvars() + @test _taylorNparams.variable_names == get_variable_names() + @test _taylorNparams.variable_symbols == get_variable_symbols() + @test eltype(set_variables(Int, "x", numvars=2, order=6)) == TaylorN{Int} @test eltype(set_variables("x", numvars=2, order=6)) == TaylorN{Float64} @test eltype(set_variables(BigInt, "x y", order=6)) == TaylorN{BigInt} @test eltype(set_variables("x y", order=6)) == TaylorN{Float64} - @test typeof(show_params_TaylorN()) == Void + @test eltype(set_variables(Int, :x, numvars=2, order=6)) == TaylorN{Int} + @test eltype(set_variables(:x, numvars=2, order=6)) == TaylorN{Float64} + @test eltype(set_variables(BigInt, [:x,:y], order=6)) == TaylorN{BigInt} + @test eltype(set_variables([:x,:y], order=6)) == TaylorN{Float64} + @test @compat typeof(show_params_TaylorN()) == Nothing + @test @compat typeof(show_monomials(2)) == Nothing @test TaylorSeries.coeff_table[2][1] == [1,0] @test TaylorSeries.index_table[2][1] == 7 @@ -23,9 +43,18 @@ using Base.Test @test get_order() == 6 @test get_numvars() == 2 + @test get_variables()[1].order == get_order() + @test get_variables(order=2)[1].order == 2 + @test get_variables(order=3)[1] == TaylorN(1,order=3) + x, y = set_variables("x y", order=6) @test x.order == 6 - @test TaylorSeries.set_variable_names(["x","y"]) == ["x", "y"] + @test TaylorSeries.set_variable_names(["x","y"]) == nothing + @test TaylorSeries._params_TaylorN_.variable_names == ["x","y"] + @test TaylorSeries._params_TaylorN_.variable_symbols == [:x, :y] + @test get_variable_symbols() == [:x, :y] + @test TaylorSeries.lookupvar(:x) == 1 + @test TaylorSeries.lookupvar(:α) == 0 @test TaylorSeries.get_variable_names() == ["x", "y"] @test x == HomogeneousPolynomial(Float64, 1) @test x == HomogeneousPolynomial(1) @@ -35,7 +64,7 @@ using Base.Test set_variables("x", numvars=2, order=17) v = [1,2] - @test typeof(TaylorSeries.resize_coeffsHP!(v,2)) == Void + @test @compat typeof(TaylorSeries.resize_coeffsHP!(v,2)) == Nothing @test v == [1,2,0] @test_throws AssertionError TaylorSeries.resize_coeffsHP!(v,1) HomogeneousPolynomial(v)[3] = 3 @@ -54,8 +83,8 @@ using Base.Test @test one(HomogeneousPolynomial(1,1)) == HomogeneousPolynomial([1,1]) uT = one(convert(TaylorN{Float64},yT)) @test uT == one(HomogeneousPolynomial) - @test zeroT[1] == HomogeneousPolynomial(0, 0) - @test uT[1] == HomogeneousPolynomial(1, 0) + @test zeroT[0] == HomogeneousPolynomial(0, 0) + @test uT[0] == HomogeneousPolynomial(1, 0) @test ones(xH,1) == [1, xH+yH] @test typeof(ones(xH,2)) == Array{HomogeneousPolynomial{Int},1} @test length(ones(xH,2)) == 3 @@ -66,10 +95,10 @@ using Base.Test @test_throws AssertionError TaylorSeries.fixorder(zeros(xH,0)[1],yH) @test get_order(zeroT) == 1 - @test xT[2][1] == 1 + @test xT[1][1] == 1 @test yH[2] == 1 - @test get_coeff(xT,[1,0]) == 1 - @test get_coeff(yH,[1,0]) == 0 + @test getcoeff(xT,(1,0)) == getcoeff(xT,[1,0]) == 1 + @test getcoeff(yH,(1,0)) == getcoeff(yH,[1,0]) == 0 @test typeof(convert(HomogeneousPolynomial,1im)) == HomogeneousPolynomial{Complex{Int}} @test convert(HomogeneousPolynomial,1im) == @@ -102,28 +131,32 @@ using Base.Test @test get_order(xT) == 17 @test xT * true == xT @test false * yT == zero(yT) + @test HomogeneousPolynomial([1.0])*xH == xH @test xT == TaylorN([xH]) @test one(xT) == TaylorN(1,5) @test TaylorN(uT) == convert(TaylorN{Complex},1) @test get_numvars() == 2 @test length(uT) == get_order()+1 - @test eltype(convert(TaylorN{Complex128},1)) == Complex128 + @test eltype(convert(TaylorN{Complex{Float64}},1)) == Complex{Float64} @test 1+xT+yT == TaylorN(1,1) + TaylorN([xH,yH],1) @test xT-yT-1 == TaylorN([-1,xH-yH]) @test xT*yT == TaylorN([HomogeneousPolynomial([0,1,0],2)]) - @test (1/(1-xT))[4] == HomogeneousPolynomial([1.0],3) + @test (1/(1-xT))[3] == HomogeneousPolynomial([1.0],3) @test xH^20 == HomogeneousPolynomial([0], get_order()) - @test (yT/(1-xT))[5] == xH^3 * yH + @test (yT/(1-xT))[4] == xH^3 * yH @test mod(1+xT,1) == +xT - @test (rem(1+xT,1))[1] == 0 + @test (rem(1+xT,1))[0] == 0 @test mod(1+xT,1.0) == +xT - @test (rem(1+xT,1.0))[1] == 0 + @test (rem(1+xT,1.0))[0] == 0 @test abs(1-xT) == 1-xT @test abs(-1-xT) == 1+xT - @test derivative(mod2pi(2pi+yT^3),2) == derivative(yT^3,2) - @test derivative(yT) == zeroT + @test derivative(yH,1) == derivative(xH, :x₂) + @test derivative(mod2pi(2pi+yT^3),2) == derivative(yT^3, :x₂) + @test derivative(yT^3, :x₂) == derivative(yT^3, (0,1)) + @test derivative(yT) == zeroT == derivative(yT, (1,0)) + @test derivative((0,1), yT) == 1 @test -xT/3im == im*xT/3 @test (xH/3im)' == im*xH/3 @test xT/BigInt(3) == TaylorN(BigFloat,1)/3 @@ -131,14 +164,51 @@ using Base.Test @test (xH/complex(0,BigInt(3)))' == im*HomogeneousPolynomial([BigInt(1),0])/3 @test evaluate(xH) == zero(eltype(xH)) + @test xH() == zero(eltype(xH)) + @test xH([1,1]) == evaluate(xH, [1,1]) + @test xH((1,1)) == 1 + hp = -5.4xH+6.89yH + @test hp([1,1]) == evaluate(hp, [1,1]) + vr = rand(2) + @test hp(vr) == evaluate(hp, vr) + + @test integrate(yH,1) == integrate(xH, :x₂) + p = (xT-yT)^6 + @test integrate(derivative(p, 1), 1, yT^6) == p + @test integrate(derivative(p, :x₁), :x₁, yT^6) == p + @test derivative(integrate(p, 2), 2) == p + @test derivative(integrate(p, :x₂), :x₂) == p + @test derivative(TaylorN(1.0, get_order())) == TaylorN(0.0, get_order()) + @test integrate(TaylorN(6.0, get_order()), 1) == 6xT + @test integrate(TaylorN(0.0, get_order()), 2) == TaylorN(0.0, get_order()) + @test integrate(TaylorN(0.0, get_order()), 2, xT) == xT + @test integrate(TaylorN(0.0, get_order()), :x₂, xT) == xT + @test integrate(xT^17, 2) == TaylorN(0.0, get_order()) + @test integrate(xT^17, 1, yT) == yT + @test integrate(xT^17, 1, 2.0) == TaylorN(2.0, get_order()) + @test integrate(xT^17, :x₁, 2.0) == TaylorN(2.0, get_order()) + @test_throws AssertionError integrate(xT, 1, xT) + @test_throws AssertionError integrate(xT, :x₁, xT) + @test_throws AssertionError derivative(xT, (1,)) + @test_throws AssertionError derivative(xT, (1,2,3)) + @test_throws AssertionError derivative(xT, (-1,2)) + @test_throws AssertionError derivative((1,), xT) + @test_throws AssertionError derivative((1,2,3), xT) + @test_throws AssertionError derivative((-1,2), xT) - @test derivative(2xT*yT^2,1) == 2yT^2 + + @test derivative(2xT*yT^2, (8,8)) == 0 + @test derivative((8,8), 2xT*yT^2) == 0 + @test derivative(2xT*yT^2, 1) == 2yT^2 + @test derivative((1,0), 2xT*yT^2) == 0 + @test derivative(2xT*yT^2, (1,2)) == 4*one(yT) + @test derivative((1,2), 2xT*yT^2) == 4 @test xT*xT^3 == xT^4 txy = 1.0 + xT*yT - 0.5*xT^2*yT + (1/3)*xT^3*yT + 0.5*xT^2*yT^2 - @test getindex((1+TaylorN(1))^TaylorN(2),1:5) == txy.coeffs[1:5] + @test getindex((1+TaylorN(1))^TaylorN(2),0:4) == txy.coeffs[1:5] @test ( (1+TaylorN(1))^TaylorN(2) )[:] == ( (1+TaylorN(1))^TaylorN(2) ).coeffs[:] @test txy.coeffs[:] == txy[:] - @test txy.coeffs[:] == txy[1:end] + @test txy.coeffs[:] == txy[0:end] txy[:] .= ( -1.0 + 3xT*yT - xT^2*yT + (4/3)*xT^3*yT + (1/3)*xT*yT^3 + 0.5*xT^2*yT^2 + 0.5*xT*yT^3 )[:] @test txy[:] == ( -1.0 + 3xT*yT - xT^2*yT + (4/3)*xT^3*yT + (1/3)*xT*yT^3 + 0.5*xT^2*yT^2 + 0.5*xT*yT^3 )[:] txy[2:end-1] .= ( 1.0 - xT*yT + 0.5*xT^2*yT - (2/3)*xT*yT^3 - 0.5*xT^2*yT^2 + 7*xT^3*yT )[2:end-1] @@ -146,15 +216,17 @@ using Base.Test a = -5.0 + sin(xT+yT^2) b = deepcopy(a) - @test a[:] == a[1:end] + @test a[:] == a[0:end] @test a[:] == b[:] @test a[1:end] == b[1:end] @test a[end][:] == b[end][:] @test a[end][1:end] == b[end][1:end] - rv = a[end][:] .= rand.() + a[end][:] .= rand.() + rv = a[end][:] @test a[end][:] == rv @test a[end][:] != b[end][:] - rv = a[end][1:end] .= rand.() + a[end][1:end] .= rand.() + rv = a[end][1:end] @test a[end][1:end] == rv @test a[end][1:end] != b[end][1:end] @@ -177,21 +249,24 @@ using Base.Test pol = sin(xT+yT*xT)+yT^2-(1-xT)^3 q = deepcopy(pol) q[:] = 0.0 + @test get_order.(q[:]) == collect(0:q.order) @test q[:] == zero(q[:]) q[:] = pol.coeffs @test q == pol @test q[:] == pol[:] q[2:end-1] = 0.0 - @test q[2:end-1] == zero(q[2:end-1]) + @test q[2:end-1] == zero.(q[2:end-1]) @test q[1] == pol[1] @test q[end] == pol[end] + # q[:] = pol.coeffs + # zH0 = zero(HomogeneousPolynomial{Float64}) + q[:] = 1.0 + @test q[1] == HomogeneousPolynomial([1,0]) + @test q[2] == HomogeneousPolynomial([1,0,0]) q[:] = pol.coeffs - zH0 = zero(HomogeneousPolynomial{Float64}) - q[:] = zH0 - @test q[:] == zero(q[:]) - q[:] = pol.coeffs - q[2:end-1] = zH0 - @test q[2:end-1] == zero(q[2:end-1]) + q[2:end-1] = one.(q[2:end-1]) + @test q[2:end-1] == one.(q[2:end-1]) + @test q[2] == HomogeneousPolynomial([1,1,1]) @test q[1] == pol[1] @test q[end] == pol[end] q[:] = pol.coeffs @@ -199,52 +274,77 @@ using Base.Test q[:] = zHall @test q[:] == zHall q[:] = pol.coeffs - q[2:end-1] = zHall[2:end-1] - @test q[2:end-1] == zHall[2:end-1] + q[1:end-1] .= zHall[2:end-1] + @test q[1:end-1] == zHall[2:end-1] q[:] = pol.coeffs @test q[:] != zeros(q.order+1) q[:] = zeros(q.order+1) @test q[:] == zeros(q.order+1) q[:] = pol.coeffs - q[2:end-1] = zeros(q.order+1)[2:end-1] + q[1:end-1] = zeros(q.order+1)[2:end-1] @test q != pol @test q[:] != pol[:] - @test q[2:end-1] == zeros(q.order+1)[2:end-1] - @test q[1] == pol[1] + @test q[1:end-1] == zeros(q.order+1)[2:end-1] + @test q[0] == pol[0] @test q[end] == pol[end] q[:] = pol.coeffs pol2 = cos(sin(xT)-yT^3*xT)-3yT^2+sqrt(1-xT) - q[3:end-2] = pol2.coeffs[3:end-2] - @test q[1:2] == pol[1:2] - @test q[3:end-2] == pol2[3:end-2] + q[2:end-2] .= pol2.coeffs[3:end-2] + @test q[0:1] == pol[0:1] + @test q[2:end-2] == pol2[2:end-2] @test q[end-1:end] == pol[end-1:end] - @test_throws DomainError yT^(-2) - @test_throws DomainError yT^(-2.0) + if VERSION < v"0.7.0-DEV" + @test_throws DomainError yT^(-2) + else + @test_throws AssertionError yT^(-2) + end + @test_throws AssertionError yT^(-2.0) @test (1+xT)^(3//2) == ((1+xT)^0.5)^3 @test real(xH) == xH @test imag(xH) == zero(xH) @test conj(im*yH) == (im*yH)' @test conj(im*yT) == (im*yT)' @test real( exp(1im * xT)) == cos(xT) - @test get_coeff(convert(TaylorN{Rational{Int}},cos(xT)),[4,0]) == + @test getcoeff(convert(TaylorN{Rational{Int}},cos(xT)),(4,0)) == 1//factorial(4) cr = convert(TaylorN{Rational{Int}},cos(xT)) - @test get_coeff(cr,[4,0]) == 1//factorial(4) + @test getcoeff(cr,(4,0)) == 1//factorial(4) @test imag((exp(yT))^(-1im)') == sin(yT) exy = exp( xT+yT ) @test evaluate(exy) == 1 - @test evaluate(exy,[0.1im,0.01im]) == exp(0.11im) - @test isapprox(evaluate(exy, [1,1]), e^2) + @test evaluate(exy, 0.1im, 0.01im) == exp(0.11im) + @test exy(0.1im, 0.01im) == exp(0.11im) + @test evaluate(exy,(0.1im, 0.01im)) == exp(0.11im) + @test exy((0.1im, 0.01im)) == exp(0.11im) + @test evaluate(exy,[0.1im, 0.01im]) == exp(0.11im) + @test exy([0.1im, 0.01im]) == exp(0.11im) + @test isapprox(evaluate(exy, (1,1)), eeuler^2) + @test exy(:x₁, 0.0) == exp(yT) txy = tan(xT+yT) - @test get_coeff(txy,[8,7]) == 929569/99225 + @test getcoeff(txy,(8,7)) == 929569/99225 ptxy = xT + yT + (1/3)*( xT^3 + yT^3 ) + xT^2*yT + xT*yT^2 - @test getindex(tan(TaylorN(1)+TaylorN(2)),1:5) == ptxy.coeffs[1:5] - @test evaluate(xH*yH,[1.0,2.0]) == 2.0 + @test getindex(tan(TaylorN(1)+TaylorN(2)),0:4) == ptxy.coeffs[1:5] + @test evaluate(xH*yH, 1.0, 2.0) == (xH*yH)(1.0, 2.0) == 2.0 + @test evaluate(xH*yH, (1.0, 2.0)) == 2.0 + @test evaluate(xH*yH, [1.0, 2.0]) == 2.0 + @test ptxy(:x₁, -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2 + @test ptxy(:x₁ => -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2 + @test evaluate(ptxy, :x₁ => -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2 + @test evaluate(ptxy, :x₁, -1.0) == -1 + yT + (-1.0+yT^3)/3 + yT - yT^2 v = zeros(Int, 2) @test evaluate!([xT, yT], ones(Int, 2), v) == nothing @test v == ones(2) + A_TN = [xT 2xT 3xT; yT 2yT 3yT] + @test evaluate(A_TN, ones(2)) == [1.0 2.0 3.0; 1.0 2.0 3.0] + @test evaluate(A_TN) == [0.0 0.0 0.0; 0.0 0.0 0.0] + @test A_TN() == [0.0 0.0 0.0; 0.0 0.0 0.0] + @test (view(A_TN,:,:))() == [0.0 0.0 0.0; 0.0 0.0 0.0] + t = Taylor1(10) + @test A_TN([t,t^2]) == [t 2t 3t; t^2 2t^2 3t^2] + @test view(A_TN, :, :)(ones(2)) == A_TN(ones(2)) + @test view(A_TN, :, 1)(ones(2)) == A_TN[:,1](ones(2)) @test evaluate(sin(asin(xT+yT)), [1.0,0.5]) == 1.5 @test evaluate(asin(sin(xT+yT)), [1.0,0.5]) == 1.5 @@ -262,71 +362,149 @@ using Base.Test xx = 1.0*zeroT TaylorSeries.add!(xx, 1.0*xT, 2yT, 1) - @test xx[2] == HomogeneousPolynomial([1,2]) + @test xx[1] == HomogeneousPolynomial([1,2]) + TaylorSeries.add!(xx, 5.0, 0) + @test xx[0] == HomogeneousPolynomial([5.0]) + TaylorSeries.add!(xx, -5.0, 1) + @test xx[1] == zero(xx[1]) TaylorSeries.subst!(xx, 1.0*xT, yT, 1) - @test xx[2] == HomogeneousPolynomial([1,-1]) + @test xx[1] == HomogeneousPolynomial([1,-1]) + TaylorSeries.subst!(xx, 5.0, 0) + @test xx[0] == HomogeneousPolynomial([-5.0]) + TaylorSeries.subst!(xx, -5.0, 1) + @test xx[1] == zero(xx[end]) TaylorSeries.div!(xx, 1.0+xT, 1.0+xT, 0) - @test xx[1] == 1.0 + @test xx[0] == HomogeneousPolynomial([1.0]) + TaylorSeries.pow!(xx, 1.0+xT, 0.5, 1) + @test xx[1] == HomogeneousPolynomial([0.5,0.0]) + xx = 1.0*zeroT TaylorSeries.pow!(xx, 1.0+xT, 1.5, 0) - @test xx[1] == 1.0 + @test xx[0] == HomogeneousPolynomial([1.0]) + TaylorSeries.pow!(xx, 1.0+xT, 1.5, 1) + @test xx[1] == HomogeneousPolynomial([1.5,0.0]) + xx = 1.0*zeroT + TaylorSeries.pow!(xx, 1.0+xT, 0, 0) + @test xx[0] == HomogeneousPolynomial([1.0]) + TaylorSeries.pow!(xx, 1.0+xT, 1, 1) + @test xx[1] == HomogeneousPolynomial([1.0,0.0]) + xx = 1.0*zeroT + TaylorSeries.pow!(xx, 1.0+xT, 2, 0) + @test xx[0] == HomogeneousPolynomial([1.0]) + TaylorSeries.pow!(xx, 1.0+xT, 2, 1) + @test xx[1] == HomogeneousPolynomial([2.0,0.0]) + xx = 1.0*zeroT TaylorSeries.sqrt!(xx, 1.0+xT, 0) - @test xx[1] == 1.0 - TaylorSeries.exp!(xx, xT, 0) - @test xx[1] == 1.0 + TaylorSeries.sqrt!(xx, 1.0+xT, 1) + @test xx[0] == 1.0 + @test xx[1] == HomogeneousPolynomial([0.5,0.0]) + xx = 1.0*zeroT + TaylorSeries.exp!(xx, 1.0*xT, 0) + TaylorSeries.exp!(xx, 1.0*xT, 1) + @test xx[0] == 1.0 + @test xx[1] == HomogeneousPolynomial([1.0,0.0]) + xx = 1.0*zeroT TaylorSeries.log!(xx, 1.0+xT, 0) - @test xx[1] == 0.0 + TaylorSeries.log!(xx, 1.0+xT, 1) + @test xx[0] == 0.0 + @test xx[1] == HomogeneousPolynomial(1.0,1) + xx = 1.0*zeroT cxx = zero(xx) TaylorSeries.sincos!(xx, cxx, 1.0*xT, 0) - @test xx[1] == 0.0 - @test cxx[1] == 1.0 + TaylorSeries.sincos!(xx, cxx, 1.0*xT, 1) + @test xx[0] == 0.0 + @test xx[1] == HomogeneousPolynomial(1.0,1) + @test cxx[0] == 1.0 + @test cxx[1] == HomogeneousPolynomial(0.0,1) + xx = 1.0*zeroT + cxx = zero(xx) TaylorSeries.tan!(xx, 1.0*xT, cxx, 0) - @test xx[1] == 0.0 - @test cxx[1] == 0.0 + TaylorSeries.tan!(xx, 1.0*xT, cxx, 1) + @test xx[0] == 0.0 + @test xx[1] == HomogeneousPolynomial(1.0,1) + @test cxx[0] == 0.0 + @test cxx[1] == HomogeneousPolynomial(0.0,1) + xx = 1.0*zeroT + cxx = zero(xx) TaylorSeries.asin!(xx, 1.0*xT, cxx, 0) - @test xx[1] == 0.0 - @test cxx[1] == 1.0 + TaylorSeries.asin!(xx, 1.0*xT, cxx, 1) + @test xx[0] == 0.0 + @test xx[1] == HomogeneousPolynomial(1.0,1) + @test cxx[0] == 1.0 + @test cxx[1] == HomogeneousPolynomial(0.0,1) + xx = 1.0*zeroT + cxx = zero(xx) TaylorSeries.acos!(xx, 1.0*xT, cxx, 0) - @test xx[1] == acos(0.0) - @test cxx[1] == 1.0 + TaylorSeries.acos!(xx, 1.0*xT, cxx, 1) + @test xx[0] == acos(0.0) + @test xx[1] == HomogeneousPolynomial(-1.0,1) + @test cxx[0] == 1.0 + @test cxx[1] == HomogeneousPolynomial(0.0,1) + xx = 1.0*zeroT + cxx = zero(xx) TaylorSeries.atan!(xx, 1.0*xT, cxx, 0) - @test xx[1] == 0.0 - @test cxx[1] == 1.0 + TaylorSeries.atan!(xx, 1.0*xT, cxx, 1) + @test xx[0] == 0.0 + @test xx[1] == HomogeneousPolynomial(1.0,1) + @test cxx[0] == 1.0 + @test cxx[1] == HomogeneousPolynomial(0.0,1) + xx = 1.0*zeroT + cxx = zero(xx) TaylorSeries.sinhcosh!(xx, cxx, 1.0*xT, 0) - @test xx[1] == 0.0 - @test cxx[1] == 1.0 + TaylorSeries.sinhcosh!(xx, cxx, 1.0*xT, 1) + @test xx[0] == 0.0 + @test xx[1] == HomogeneousPolynomial(1.0,1) + @test cxx[0] == 1.0 + @test cxx[1] == HomogeneousPolynomial(0.0,1) + xx = 1.0*zeroT + cxx = zero(xx) TaylorSeries.tanh!(xx, 1.0*xT, cxx, 0) - @test xx[1] == 0.0 - @test cxx[1] == 0.0 + TaylorSeries.tanh!(xx, 1.0*xT, cxx, 1) + @test xx[0] == 0.0 + @test xx[1] == HomogeneousPolynomial(1.0,1) + @test cxx[0] == 0.0 + @test cxx[1] == HomogeneousPolynomial(0.0,1) - g1(xT,yT) = xT^3 + 3yT^2 - 2xT^2 * yT - 7xT + 2 - g2(xT,yT) = yT + xT^2 - xT^4 - f1 = g1(xT,yT) - f2 = g2(xT,yT) + g1(x, y) = x^3 + 3y^2 - 2x^2 * y - 7x + 2 + g2(x, y) = y + x^2 - x^4 + f1 = g1(xT, yT) + f2 = g2(xT, yT) @test gradient(f1) == [ 3*xT^2-4*xT*yT-TaylorN(7,0), 6*yT-2*xT^2 ] @test ∇(f2) == [2*xT - 4*xT^3, TaylorN(1,0)] @test jacobian([f1,f2], [2,1]) == jacobian( [g1(xT+2,yT+1), g2(xT+2,yT+1)] ) - jac = Array{Int64}(2, 2) + @compat jac = Array{Int64}(undef, 2, 2) jacobian!(jac, [g1(xT+2,yT+1), g2(xT+2,yT+1)]) @test jac == jacobian( [g1(xT+2,yT+1), g2(xT+2,yT+1)] ) jacobian!(jac, [f1,f2], [2,1]) @test jac == jacobian([f1,f2], [2,1]) - @test [xT yT]*hessian(f1*f2)*[xT, yT] == [ 2*TaylorN((f1*f2)[3]) ] + @test hessian( f1*f2 ) == + [derivative((2,0), f1*f2) derivative((1,1), (f1*f2)); + derivative((1,1), f1*f2) derivative((0,2), (f1*f2))] == [4 -7; -7 0] + @test hessian( f1*f2, [xT, yT] ) == + [derivative(f1*f2, (2,0)) derivative((f1*f2), (1,1)); + derivative(f1*f2, (1,1)) derivative((f1*f2), (0,2))] + @test [xT yT]*hessian(f1*f2)*[xT, yT] == [ 2*TaylorN((f1*f2)[2]) ] @test hessian(f1^2)/2 == [ [49,0] [0,12] ] @test hessian(f1-f2-2*f1*f2) == (hessian(f1-f2-2*f1*f2))' @test hessian(f1-f2,[1,-1]) == hessian(g1(xT+1,yT-1)-g2(xT+1,yT-1)) - hes = Array{Int64}(2, 2) + @compat hes = Array{Int64}(undef, 2, 2) hessian!(hes, f1*f2) @test hes == hessian(f1*f2) - @test [xT yT]*hes*[xT, yT] == [ 2*TaylorN((f1*f2)[3]) ] + @test [xT yT]*hes*[xT, yT] == [ 2*TaylorN((f1*f2)[2]) ] hessian!(hes, f1^2) @test hes/2 == [ [49,0] [0,12] ] hessian!(hes, f1-f2-2*f1*f2) @test hes == hes' - hes1 = hes2 = Array{Int64}(2, 2) - hessian!(hes1,f1-f2,[1,-1]) - hessian!(hes2,g1(xT+1,yT-1)-g2(xT+1,yT-1)) - @test hes1 == hes2 + @compat hes1 = Array{Int64}(undef, 2, 2) + hessian!(hes1, f1-f2,[1,-1]) + hessian!(hes, g1(xT+1,yT-1)-g2(xT+1,yT-1)) + @test hes1 == hes + displayBigO(false) + @test string(-xH) == " - 1 x₁" + @test string(xT^2) == " 1 x₁²" + @test string(1im*yT) == " ( 1 im ) x₂" + @test string(xT-im*yT) == " ( 1 ) x₁ - ( 1 im ) x₂" + displayBigO(true) @test string(-xH) == " - 1 x₁" @test string(xT^2) == " 1 x₁² + 𝒪(‖x‖¹⁸)" @test string(1im*yT) == " ( 1 im ) x₂ + 𝒪(‖x‖¹⁸)" @@ -341,4 +519,119 @@ using Base.Test @test_throws AssertionError cos(x)/sin(y) @test_throws BoundsError xH[20] @test_throws BoundsError xT[20] + + a = 3x + 4y +6x^2 + 8x*y + @test typeof( norm(x) ) == Float64 + @test norm(x) > 0 + @test norm(a) == norm([3,4,6,8.0]) + @test norm(a, 4) == sum([3,4,6,8.0].^4)^(1/4.) + @test norm(a, Inf) == 8.0 + @test norm((3.0 + 4im)*x) == abs(3.0 + 4im) + + @test TaylorSeries.rtoldefault(TaylorN{Int64}) == 0 + @test TaylorSeries.rtoldefault(TaylorN{Float64}) == sqrt(eps(Float64)) + @test TaylorSeries.rtoldefault(TaylorN{BigFloat}) == sqrt(eps(BigFloat)) + @test TaylorSeries.real(TaylorN{Float64}) == TaylorN{Float64} + @test TaylorSeries.real(TaylorN{Complex{Float64}}) == TaylorN{Float64} + @test isfinite(a) + @test a[0] ≈ a[0] + @test a[1] ≈ a[1] + @test a[2] ≈ a[2] + @test a[3] ≈ a[3] + @test a ≈ a + @test a .≈ a + b = deepcopy(a) + b[2][3] = Inf + @test !isfinite(b) + b[2][3] = NaN + @test !isfinite(b) + b[2][3] = a[2][3]+eps() + @test isapprox(a[2], b[2], rtol=eps()) + @test a ≈ b + b[2][2] = a[2][2]+sqrt(eps()) + @test a[2] ≈ b[2] + @test a ≈ b + + f11(a,b) = (a+b)^a - cos(a*b)*b + f22(a) = (a[1] + a[2])^a[1] - cos(a[1]*a[2])*a[2] + @test taylor_expand(f11, 1.0,2.0) == taylor_expand(f22, [1,2.0]) + @test evaluate(taylor_expand(x->x[1] + x[2], [1,2])) == 3.0 + f33(x,y) = 3x+y + @test eltype(taylor_expand(f33,1,1)) == eltype(1) + x,y = get_variables() + xysq = x^2 + y^2 + update!(xysq,[1.0,-2.0]) + @test xysq == (x+1.0)^2 + (y-2.0)^2 + update!(xysq,[-1.0,2.0]) + @test xysq == x^2 + y^2 + + #test function-like behavior for TaylorN + @test exy() == 1 + @test exy([0.1im,0.01im]) == exp(0.11im) + @test isapprox(exy([1,1]), eeuler^2) + @test sin(asin(xT+yT))([1.0,0.5]) == 1.5 + @test asin(sin(xT+yT))([1.0,0.5]) == 1.5 + @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(rand(2)) == 1 + @test ( -sinh(xT+yT)^2 + cosh(xT+yT)^2 )(zeros(2)) == 1 + #number of variables changed to 4... + dx = set_variables("x", numvars=4, order=10) + P = sin.(dx) + v = [1.0,2,3,4] + for i in 1:4 + @test P[i](v) == evaluate(P[i], v) + end + @test P.(fill(v, 4)) == fill(P(v), 4) + F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))] + Q = F(v+dx) + @test Q.( fill(v, 4) ) == fill(Q(v), 4) + vr = map(x->rand(4), 1:4) + @test Q.(vr) == map(x->Q(x), vr) + for i in 1:4 + @test P[i]() == evaluate(P[i]) + @test Q[i]() == evaluate(Q[i]) + end + @test P() == evaluate.(P) + @test P() == evaluate(P) + @test Q() == evaluate.(Q) + @test Q() == evaluate(Q) + @test Q[1:3]() == evaluate(Q[1:3]) + + dx = set_variables("x", numvars=4, order=10) + for i in 1:4 + @test deg2rad(180+dx[i]) == pi + deg2rad(1.0)dx[i] + @test rad2deg(pi+dx[i]) == 180.0+rad2deg(1.0)dx[i] + end + p = sin(exp(dx[1]*dx[2])+dx[3]*dx[2])/(1.0+dx[4]^2) + q = zero(p) + TaylorSeries.deg2rad!(q, p, 0) + @test q[0] == p[0]*(pi/180) + TaylorSeries.deg2rad!.(q, p, [1,3,5]) + for i in [0,1,3,5] + @test q[i] == p[i]*(pi/180) + end + TaylorSeries.rad2deg!(q, p, 0) + @test q[0] == p[0]*(180/pi) + TaylorSeries.rad2deg!.(q, p, [1,3,5]) + for i in [0,1,3,5] + @test q[i] == p[i]*(180/pi) + end + xT = 5+TaylorN(Int64, 1, order=10) + yT = TaylorN(2, order=10) + TaylorSeries.deg2rad!(yT, xT, 0) + @test yT[0] == xT[0]*(pi/180) + TaylorSeries.rad2deg!(yT, xT, 0) + @test yT[0] == xT[0]*(180/pi) +end + +@testset "Integrate for several variables" begin + + t, x, y = set_variables("t x y") + + @test integrate(t, 1) == 0.5*t^2 + @test integrate(t, 2) == t * x + @test integrate(t, 3) == t * y + @test integrate(x, 1) == t * x + @test integrate(x, 2) == 0.5*x^2 + @test integrate(y, 2) == x * y + end diff --git a/test/mixtures.jl b/test/mixtures.jl index 6ce96a9e..386cb493 100644 --- a/test/mixtures.jl +++ b/test/mixtures.jl @@ -2,7 +2,14 @@ # using TaylorSeries -using Base.Test +using Compat + +if VERSION < v"0.7.0-DEV.2004" + using Base.Test +else + using Test + using LinearAlgebra, SparseArrays +end @testset "Tests with mixtures of Taylor1 and TaylorN" begin @test TaylorSeries.NumberNotSeries == Union{Real,Complex} @@ -42,12 +49,16 @@ using Base.Test @test 3*xHt == HomogeneousPolynomial([3*one(t), zero(t)]) @test t*xHt == HomogeneousPolynomial([t, zero(t)]) @test complex(0,1)*xHt == HomogeneousPolynomial([1im*one(t), zero(1im*t)]) - @test eltype(complex(0,1)*xHt) == Taylor1{Complex128} + @test eltype(complex(0,1)*xHt) == Taylor1{Complex{Float64}} + @test (xHt+yHt)(1, 1) == 1+t + @test (xHt+yHt)([1, 1]) == (xHt+yHt)((1, 1)) tN1 = TaylorN([HomogeneousPolynomial([t]),xHt,yHt^2]) - @test tN1[1] == HomogeneousPolynomial([t]) + @test tN1[0] == HomogeneousPolynomial([t]) + @test tN1(t,one(t)) == 2t+t^2 + @test tN1([t,one(t)]) == tN1((t,one(t))) t1N = convert(Taylor1{TaylorN{Float64}}, tN1) - @test t1N[1] == HomogeneousPolynomial(1) + @test t1N[0] == HomogeneousPolynomial(1) ctN1 = convert(TaylorN{Taylor1{Float64}}, t1N) @test eltype(xHt) == Taylor1{Float64} @test eltype(tN1) == Taylor1{Float64} @@ -73,16 +84,16 @@ using Base.Test @test mod(tN1+1,1.0) == 0+tN1 @test mod(tN1-1.125,2) == 0.875+tN1 - @test (rem(tN1+1.125,1.0))[1][1] == 0.125 + t - @test (rem(tN1-1.125,2))[1][1] == -1.125 + t - @test mod2pi(-3pi+tN1)[1][1][1] ≈ pi - @test mod2pi(0.125+2pi+tN1)[1][1][1] ≈ 0.125 + @test (rem(tN1+1.125,1.0))[0][1] == 0.125 + t + @test (rem(tN1-1.125,2))[0][1] == -1.125 + t + @test mod2pi(-3pi+tN1)[0][1][0] ≈ pi + @test mod2pi(0.125+2pi+tN1)[0][1][0] ≈ 0.125 @test mod(t1N+1.125,1.0) == 0.125+t1N @test mod(t1N-1.125,2) == 0.875+t1N - @test (rem(t1N+1.125,1.0))[1] == 0.125 + t1N[1] - @test (rem(t1N-1.125,2))[1] == -1.125 + t1N[1] - @test mod2pi(-3pi+t1N)[1][1][1] ≈ pi - @test mod2pi(0.125+2pi+t1N)[1][1][1] ≈ 0.125 + @test (rem(t1N+1.125,1.0))[0] == 0.125 + t1N[0] + @test (rem(t1N-1.125,2))[0] == -1.125 + t1N[0] + @test mod2pi(-3pi+t1N)[0][0][1] ≈ pi + @test mod2pi(0.125+2pi+t1N)[0][0][1] ≈ 0.125 @test abs(tN1+1) == 1+tN1 @test abs(tN1-1) == 1-tN1 @@ -94,16 +105,22 @@ using Base.Test @test convert(Array{TaylorN{Taylor1{Float64}},1}, [t1N, t1N]) == [tN1, tN1] @test convert(Array{TaylorN{Taylor1{Float64}},2}, [t1N t1N]) == [tN1 tN1] + @test evaluate(t1N, 0.0) == TaylorN(xH, 2) + @test t1N() == TaylorN(xH, 2) @test string(evaluate(t1N, 0.0)) == " 1.0 x₁ + 𝒪(‖x‖³)" @test string(evaluate(t1N^2, 1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)" + @test string((t1N^2)(1.0)) == " 1.0 + 2.0 x₁ + 1.0 x₁² + 2.0 x₂² + 𝒪(‖x‖³)" v = zeros(TaylorN{Float64},2) @test evaluate!([t1N, t1N^2], 0.0, v) == nothing - @test v[1] == TaylorN([xHt]) - @test v[2] == TaylorN([xHt^2]) + @test v == [TaylorN(1), TaylorN(1)^2] + + vt = zeros(Taylor1{Float64},2) + @test evaluate!([tN1, tN1^2], [t, t], vt) == nothing + @test vt == [2t, 4t^2] tint = Taylor1(Int, 10) t = Taylor1(10) - x = TaylorN( [HomogeneousPolynomial(zero(t)), HomogeneousPolynomial([one(t),zero(t)])], 5) + x = TaylorN( [HomogeneousPolynomial(zero(t), 5), HomogeneousPolynomial([one(t),zero(t)])], 5) y = TaylorN(typeof(tint), 2, order=5) @test typeof(x) == TaylorN{Taylor1{Float64}} @test eltype(y) == Taylor1{Int} @@ -133,4 +150,127 @@ using Base.Test @test xx*δx + Taylor1(typeof(δx),5) == δx + δx^2 + Taylor1(typeof(δx),5) @test xx/(1+δx) == one(xx) @test typeof(xx+δx) == Taylor1{TaylorN{Float64}} + + #testing evaluate and function-like behavior of Taylor1, TaylorN for mixtures: + t = Taylor1(25) + p = cos(t) + q = sin(t) + a = [p,q] + dx = set_variables("x", numvars=4, order=10) + P = sin.(dx) + v = [1.0,2,3,4] + F(x) = [sin(sin(x[4]+x[3])), sin(cos(x[3]-x[2])), cos(sin(x[1]^2+x[2]^2)), cos(cos(x[2]*x[3]))] + Q = F(v+dx) + diff_evals = cos(sin(dx[1]))-p(P[1]) + @test norm( norm.(map(x->x.coeffs, diff_evals.coeffs),Inf) , Inf) < 1e-15 + #evaluate a Taylor1 at a TaylorN + @test p(P) == evaluate(p, P) + @test q(Q) == evaluate(q, Q) + #evaluate an array of Taylor1s at a TaylorN + aT1 = [p,q,p^2,log(1+q)] #an array of Taylor1s + @test aT1(Q[4]) == evaluate(aT1, Q[4]) + @test (aT1.^2)(Q[3]) == evaluate(aT1.^2, Q[3]) + #evaluate a TaylorN at an array of Taylor1s + @test P[1](aT1) == evaluate(P[1], aT1) + @test P[1](aT1) == evaluate(P[1], (aT1...,)) + @test Q[2](aT1) == evaluate(Q[2], [aT1...]) + #evaluate an array of TaylorN{Float64} at an array of Taylor1{Float64} + @test P(aT1) == evaluate(P, aT1) + @test Q(aT1) == evaluate(Q, aT1) + #test evaluation of an Array{TaylorN{Taylor1}} at an Array{Taylor1} + aH1 = [ + HomogeneousPolynomial([Taylor1(rand(2))]), + HomogeneousPolynomial([Taylor1(rand(2)),Taylor1(rand(2)), + Taylor1(rand(2)),Taylor1(rand(2))]) + ] + bH1 = [ + HomogeneousPolynomial([Taylor1(rand(2))]), + HomogeneousPolynomial([Taylor1(rand(2)),Taylor1(rand(2)), + Taylor1(rand(2)),Taylor1(rand(2))]) + ] + aTN1 = TaylorN(aH1); bTN1 = TaylorN(bH1) + x = [aTN1, bTN1] + δx = [Taylor1(rand(3)) for i in 1:4] + @test typeof(x) == Array{TaylorN{Taylor1{Float64}},1} + @test typeof(δx) == Array{Taylor1{Float64},1} + @compat x0 = Array{Taylor1{Float64}}(undef, length(x)) + eval_x_δx = evaluate(x,δx) + @test x(δx) == eval_x_δx + evaluate!(x,δx,x0) + @test x0 == eval_x_δx + @test typeof(evaluate(x[1],δx)) == Taylor1{Float64} + @test x() == map(y->y[0][1], x) + for i in eachindex(x) + @test evaluate(x[i],δx) == eval_x_δx[i] + @test x[i](δx) == eval_x_δx[i] + end + p11 = Taylor1([sin(t),cos(t)]) + @test evaluate(p11,t) == sin(t)+t*cos(t) + @test p11(t) == sin(t)+t*cos(t) + a11 = Taylor1([t,t^2,exp(-t),sin(t),cos(t)]) + b11 = t+t*(t^2)+(t^2)*(exp(-t))+(t^3)*sin(t)+(t^4)*cos(t) + diff_a11b11 = a11(t)-b11 + @test norm(diff_a11b11.coeffs,Inf) < 1E-19 + + X, Y = set_variables(Taylor1{Float64}, "x y") + @test typeof( norm(X) ) == Float64 + @test norm(X) > 0 + @test norm(X+Y) == sqrt(2) + @test norm(-10X+4Y,Inf) == 10. + + + X,Y = convert(Taylor1{TaylorN{Float64}},X), convert(Taylor1{TaylorN{Float64}},Y) + @test typeof( norm(X) ) == Float64 + @test norm(X) > 0 + @test norm(X+Y) == sqrt(2) + @test norm(-10X+4Y,Inf) == 10. + + + @test TaylorSeries.rtoldefault(TaylorN{Taylor1{Int64}}) == 0 + @test TaylorSeries.rtoldefault(Taylor1{TaylorN{Int64}}) == 0 + for T in (Float64, BigFloat) + @test TaylorSeries.rtoldefault(TaylorN{Taylor1{T}}) == sqrt(eps(T)) + @test TaylorSeries.rtoldefault(Taylor1{TaylorN{T}}) == sqrt(eps(T)) + @test TaylorSeries.real(TaylorN{Taylor1{T}}) == TaylorN{Taylor1{T}} + @test TaylorSeries.real(Taylor1{TaylorN{T}}) == Taylor1{TaylorN{T}} + @test TaylorSeries.real(TaylorN{Taylor1{Complex{T}}}) == TaylorN{Taylor1{T}} + @test TaylorSeries.real(Taylor1{TaylorN{Complex{T}}}) == Taylor1{TaylorN{T}} + end + + rndT1(ord1) = Taylor1(-1 .+ 2rand(ord1+1)) # generates a random Taylor1 with order `ord` + nmonod(s, d) = binomial(d+s-1, d) #number of monomials in s variables with exact degree d + #rndHP generates a random `ordHP`-th order homog. pol. of Taylor1s, each with order `ord1` + rndHP(ordHP, ord1) = HomogeneousPolynomial( [rndT1(ord1) for i in 1:nmonod(get_numvars(), ordHP)] ) + #rndTN generates a random `ordHP`-th order TaylorN of of Taylor1s, each with order `ord1` + rndTN(ordN, ord1) = TaylorN([rndHP(i, ord1) for i in 0:ordN]) + + P = rndTN(get_order(), 3) + @test P ≈ P + Q = deepcopy(P) + Q[2][2] = Taylor1([NaN, Inf]) + @test isnan(Q) + @test isinf(Q) + @test !isfinite(Q) + Q[2][2] = P[2][2]+sqrt(eps())/2 + @test isapprox(P, Q, rtol=1.0) + Q[2][2] = P[2][2]+10sqrt(eps()) + @test !isapprox(P, Q, atol=sqrt(eps()), rtol=0) + @test P ≉ Q^2 + Q[2][2] = P[2][2]+eps()/2 + @test isapprox(Q, Q, atol=eps(), rtol=0) + @test isapprox(Q, P, atol=eps(), rtol=0) + Q[2][1] = P[2][1]-10eps() + @test !isapprox(Q, P, atol=eps(), rtol=0) + @test P ≉ Q^2 + + X, Y = set_variables(BigFloat, "x y", numvars=2, order=6) + p1N = Taylor1([X^2,X*Y,Y+X,Y^2]) + q1N = Taylor1([X^2,(1.0+sqrt(eps(BigFloat)))*X*Y,Y+X,Y^2]) + @test p1N ≈ p1N + @test p1N ≈ q1N + + Pv = [rndTN(get_order(), 3), rndTN(get_order(), 3)] + Qv = convert.(Taylor1{TaylorN{Float64}}, Pv) + + @test jacobian(Pv) == jacobian(Qv) end diff --git a/test/mutatingfuncts.jl b/test/mutatingfuncts.jl new file mode 100644 index 00000000..ab2f7860 --- /dev/null +++ b/test/mutatingfuncts.jl @@ -0,0 +1,93 @@ +# This file is part of TaylorSeries.jl, MIT licensed +# + +using TaylorSeries +if VERSION < v"0.7.0-DEV.2004" + using Base.Test +else + using Test +end + +@testset "Mutating functions" begin + t1 = Taylor1(6) + + # Dictionaries with calls + @test length(TaylorSeries._dict_binary_ops) == 5 + @test length(TaylorSeries._dict_unary_ops) == 22 + + @test all([haskey(TaylorSeries._dict_binary_ops, op) + for op in [:+, :-, :*, :/, :^]]) + @test all([haskey(TaylorSeries._dict_binary_calls, op) + for op in [:+, :-, :*, :/, :^]]) + + for kk in keys(TaylorSeries._dict_binary_ops) + res = TaylorSeries._internalmutfunc_call( + TaylorSeries._InternalMutFuncs(TaylorSeries._dict_binary_ops[kk]) ) + @test TaylorSeries._dict_binary_calls[kk] == res + end + for kk in keys(TaylorSeries._dict_unary_ops) + res = TaylorSeries._internalmutfunc_call( + TaylorSeries._InternalMutFuncs(TaylorSeries._dict_unary_ops[kk]) ) + @test TaylorSeries._dict_unary_calls[kk] == res + end + + # Some examples + t1 = Taylor1(5) + t2 = zero(t1) + TaylorSeries.pow!(t2, t1, 2, 2) + @test t2[2] == 1.0 + # + res = zero(t1) + TaylorSeries.add!(res, t1, t2, 3) + @test res[3] == 0.0 + TaylorSeries.add!(res, 1, t2, 3) + @test res[3] == 0.0 + TaylorSeries.add!(res, t2, 3, 0) + @test res[0] == 3.0 + TaylorSeries.subst!(res, t1, t2, 2) + @test res[2] == -1.0 + TaylorSeries.subst!(res, t1, 1, 0) + @test res[0] == -1.0 + @test res[2] == -1.0 + TaylorSeries.subst!(res, 1, t2, 2) + @test res[2] == -1.0 + + res = zero(t1) + TaylorSeries.mul!(res, t1, t2, 3) + @test res[3] == 1.0 + TaylorSeries.mul!(res, res, 2, 3) + @test res[3] == 2.0 + TaylorSeries.mul!(res, 0.5, res, 3) + @test res[3] == 1.0 + + res = zero(t1) + TaylorSeries.div!(res, t2-1, 1+t1, 0) + TaylorSeries.div!(res, t2-1, 1+t1, 1) + @test res == t1-1 + TaylorSeries.div!(res, res, 2, 0) + @test res[0] == -0.5 + + res = zero(t1) + TaylorSeries.identity!(res, t1, 0) + @test res[0] == t1[0] + TaylorSeries.zero!(res, t1, 0) + TaylorSeries.zero!(res, t1, 1) + @test res[0] == zero(t1[0]) + @test res[1] == zero(t1[1]) + TaylorSeries.one!(res, t1, 0) + TaylorSeries.one!(res, t1, 0) + @test res[0] == one(t1[0]) + @test res[1] == zero(t1[1]) + + res = zero(t1) + TaylorSeries.abs!(res, -1-t2, 2) + @test res[2] == 1.0 + @test_throws ArgumentError TaylorSeries.abs!(res, t2, 2) + + res = zero(t1) + TaylorSeries.abs2!(res, 1-t1, 1) + @test res[1] == -2.0 + TaylorSeries.abs2!(res, t1, 2) + @test res[2] == 1.0 + +end diff --git a/test/onevariable.jl b/test/onevariable.jl index e102b8bf..37a26b16 100644 --- a/test/onevariable.jl +++ b/test/onevariable.jl @@ -2,7 +2,16 @@ # using TaylorSeries -using Base.Test +using Compat + +if VERSION < v"0.7.0-DEV.2004" + using Base.Test + eeuler = Base.e +else + using Test + using LinearAlgebra, SparseArrays + eeuler = Base.MathConstants.e +end @testset "Tests for Taylor1 expansions" begin ta(a) = Taylor1([a,one(a)],15) @@ -18,34 +27,40 @@ using Base.Test @test Taylor1([1,2,3,4,5], 2) == Taylor1([1,2,3]) @test get_order(Taylor1([1,2,3,4,5], 2)) == 2 + # ta = Taylor1([1,2,3,4,5]) + # vT = [h for h in ta] + # voT = [h.order for h in Taylor1([1,2,3,4,5])] + # @test vT == [Taylor1(1,0), Taylor1(2,1), Taylor1(3,2), Taylor1(4,3), Taylor1(5,4)] + # @test voT == [0, 1, 2, 3, 4] + v = [1,2] - @test typeof(TaylorSeries.resize_coeffs1!(v,3)) == Void + @test @compat typeof(TaylorSeries.resize_coeffs1!(v,3)) == Nothing @test v == [1,2,0,0] TaylorSeries.resize_coeffs1!(v,0) @test v == [1] TaylorSeries.resize_coeffs1!(v,3) - setindex!(Taylor1(v),3,3) + setindex!(Taylor1(v),3,2) @test v == [1,0,3,0] @test Taylor1(v)[:] == [1,0,3,0] @test Taylor1(v)[:] == Taylor1(v).coeffs[:] - setindex!(Taylor1(v),0,1:3) + setindex!(Taylor1(v),0,0:2) @test v == zero(v) setindex!(Taylor1(v),1,:) - @test v == [1,1,1,1] + @test v == ones(Int, 4) Taylor1(v)[:] = 0 @test v == zero(v) rv = [rand(0:3) for i in 1:4] - @test Taylor1(rv)[:] == Taylor1(rv)[1:end] + @test Taylor1(rv)[:] == rv y = sin(Taylor1(16)) - @test y[:] == sin(Taylor1(16))[:] + @test y[:] == y.coeffs y[:] .= cos(Taylor1(16))[:] @test y == cos(Taylor1(16)) @test y[:] == cos(Taylor1(16))[:] y = sin(Taylor1(16)) rv = rand(5) - y[1:5] = rv - @test y[1:5] == rv - @test y[6:end] == sin(Taylor1(16))[6:end] + y[0:4] = rv + @test y[0:4] == rv + @test y[5:end] == y.coeffs[6:end] rv = rand( length(y.coeffs) ) y[:] = rv @test y[:] == rv @@ -55,19 +70,21 @@ using Base.Test y[:] = 0.0 @test y[:] == zero(y[:]) y = sin(Taylor1(16)) - rv = y[1:5] .= rand.() - @test y[1:5] == rv + rv = rand.(length(0:4)) + y[0:4] .= rv + @test y[0:4] == rv @test y[6:end] == sin(Taylor1(16))[6:end] - rv = y[:] .= rand.() + rv = rand.(length(y)) + y[:] = rv @test y[:] == rv @test Taylor1([0,1,0,0]) == Taylor1(3) - @test get_coeff(Taylor1(Complex128,3),1) == complex(1.0,0.0) - @test Taylor1(Complex128,3)[2] == complex(1.0,0.0) - @test getindex(Taylor1(3),2) == 1.0 - @inferred convert(Taylor1{Complex128},ot) == Taylor1{Complex128} - @test eltype(convert(Taylor1{Complex128},ot)) == Complex128 - @test eltype(convert(Taylor1{Complex128},1)) == Complex128 + @test getcoeff(Taylor1(Complex{Float64},3),1) == complex(1.0,0.0) + @test Taylor1(Complex{Float64},3)[1] == complex(1.0,0.0) + @test getindex(Taylor1(3),1) == 1.0 + @inferred convert(Taylor1{Complex{Float64}},ot) == Taylor1{Complex{Float64}} + @test eltype(convert(Taylor1{Complex{Float64}},ot)) == Complex{Float64} + @test eltype(convert(Taylor1{Complex{Float64}},1)) == Complex{Float64} @test eltype(convert(Taylor1, 1im)) == Complex{Int} @test convert(Taylor1, 1im) == Taylor1(1im, 0) @test convert(Taylor1{Int},[0,2]) == 2*t @@ -84,7 +101,13 @@ using Base.Test @test length(TaylorSeries.fixorder(zt,Taylor1([1]))[2]) == 16 @test eltype(TaylorSeries.fixorder(zt,Taylor1([1]))[1]) == Int @test TaylorSeries.findfirst(t) == 1 + @test TaylorSeries.findfirst(t^2) == 2 + @test TaylorSeries.findfirst(ot) == 0 @test TaylorSeries.findfirst(zt) == -1 + @test TaylorSeries.findlast(t) == 1 + @test TaylorSeries.findlast(t^2) == 2 + @test TaylorSeries.findlast(ot) == 0 + @test TaylorSeries.findlast(zt) == -1 @test iszero(zero(t)) @test !iszero(one(t)) @test isinf(Taylor1([typemax(1.0)])) @@ -92,7 +115,7 @@ using Base.Test @test ot == 1 @test 0.0 == zt - @test get_coeff(tim,1) == complex(0,1) + @test getcoeff(tim,1) == complex(0,1) @test zt+1.0 == ot @test 1.0-ot == zt @test t+t == 2t @@ -119,19 +142,26 @@ using Base.Test @test (t+im)^3 == Taylor1([-1im,-3,3im,1],15) @test (t+im)^4 == Taylor1([1,-4im,-6,4im,1],15) @test imag(tsquare+2im*t-1) == 2t - @test (Rational(1,2)*tsquare)[3] == 1//2 + @test (Rational(1,2)*tsquare)[2] == 1//2 @test t^2/tsquare == ot - @test ((1+t)^(1/3))[3]+1/9 ≤ tol1 - @test 1-tsquare == (1+t)-t*(1+t) + @test ((1+t)^(1/3))[2]+1/9 ≤ tol1 + @test (1.0-tsquare)^3 == (1.0-t)^3*(1.0+t)^3 @test (1-tsquare)^2 == (1+t)^2.0 * (1-t)^2.0 - @test (sqrt(1+t))[3] == -1/8 + @test (sqrt(1+t))[2] == -1/8 @test ((1-tsquare)^(1//2))^2 == 1-tsquare - @test ((1-t)^(1//4))[15] == -4188908511//549755813888 - @test abs(((1+t)^3.2)[14] + 5.4021062656e-5) < tol1 + @test ((1-t)^(1//4))[14] == -4188908511//549755813888 + @test abs(((1+t)^3.2)[13] + 5.4021062656e-5) < tol1 @test Taylor1(BigFloat,5)/6 == 1im*Taylor1(5)/complex(0,BigInt(6)) @test Taylor1(BigFloat,5)/(6*Taylor1(3)) == 1/BigInt(6) @test Taylor1(BigFloat,5)/(6im*Taylor1(3)) == -1im/BigInt(6) + # These tests involve some sort of factorization + @test t/(t+t^2) == 1/(1+t) + @test sqrt(t^2+t^3) == t*sqrt(1+t) + @test (t^3+t^4)^(1/3) ≈ t*(1+t)^(1/3) + @test norm((t^3+t^4)^(1/3) - t*(1+t)^(1/3), Inf) < eps() + @test ((t^3+t^4)^(1/3))[15] ≈ -8617640/1162261467 + trational = ta(0//1) @inferred ta(0//1) == Taylor1{Rational{Int}} @test eltype(trational) == Rational{Int} @@ -141,15 +171,27 @@ using Base.Test @test trational^3/complex(7,1) == Taylor1([0,0,0,complex(7//50,-1//50)],15) @test sqrt(zero(t)) == zero(t) - @test isapprox( rem(4.1 + t,4)[1], 0.1 ) - @test isapprox( mod(4.1 + t,4)[1], 0.1 ) - @test isapprox( rem(1+Taylor1(Int,4),4.0)[1], 1.0 ) - @test isapprox( mod(1+Taylor1(Int,4),4.0)[1], 1.0 ) - @test isapprox( mod2pi(2pi+0.1+t)[1], 0.1 ) + @test isapprox( rem(4.1 + t,4)[0], 0.1 ) + @test isapprox( mod(4.1 + t,4)[0], 0.1 ) + @test isapprox( rem(1+Taylor1(Int,4),4.0)[0], 1.0 ) + @test isapprox( mod(1+Taylor1(Int,4),4.0)[0], 1.0 ) + @test isapprox( mod2pi(2pi+0.1+t)[0], 0.1 ) @test abs(ta(1)) == ta(1) @test abs(ta(-1.0)) == -ta(-1.0) + + @test taylor_expand(x->2x,order=10) == 2*Taylor1(10) + @test taylor_expand(x->x^2+1) == Taylor1(15)*Taylor1(15) + 1 + @test evaluate(taylor_expand(cos,0.)) == cos(0.) + @test evaluate(taylor_expand(tan,pi/4)) == tan(pi/4) + @test eltype(taylor_expand(x->x^2+1,1)) == eltype(1) + tsq = t^2 + update!(tsq,2.0) + @test tsq == (t+2.0)^2 + update!(tsq,-2.0) + @test tsq == t^2 + @test log(exp(tsquare)) == tsquare @test exp(log(1-tsquare)) == 1-tsquare @test log((1-t)^2) == 2*log(1-t) @@ -158,12 +200,60 @@ using Base.Test @test exp(conj(tim)) == cos(t)-im*sin(t) == exp(tim') @test (exp(t))^(2im) == cos(2t)+im*sin(2t) @test (exp(t))^Taylor1([-5.2im]) == cos(5.2t)-im*sin(5.2t) - @test get_coeff(convert(Taylor1{Rational{Int}},cos(t)),8) == 1//factorial(8) - @test abs((tan(t))[8]- 17/315) < tol1 - @test abs((tan(t))[14]- 21844/6081075) < tol1 - @test evaluate(exp(Taylor1([0,1],17)),1.0) == 1.0*e + @test getcoeff(convert(Taylor1{Rational{Int}},cos(t)),8) == 1//factorial(8) + @test abs((tan(t))[7]- 17/315) < tol1 + @test abs((tan(t))[13]- 21844/6081075) < tol1 + @test evaluate(exp(Taylor1([0,1],17)),1.0) == 1.0*eeuler @test evaluate(exp(Taylor1([0,1],1))) == 1.0 @test evaluate(exp(t),t^2) == exp(t^2) + #Test function-like behavior for Taylor1s + t17 = Taylor1([0,1],17) + myexpfun = exp(t17) + @test myexpfun(1.0) == 1.0*eeuler + @test myexpfun() == 1.0 + @test myexpfun(t17^2) == exp(t17^2) + @test exp(t17^2)(t17) == exp(t17^2) + p = cos(t17) + q = sin(t17) + @test cos(-im*t)(1)+im*sin(-im*t)(1) == exp(-im*t)(im) + @test p(-im*t17)(1)+im*q(-im*t17)(1) == exp(-im*t17)(im) + cossin1 = x->p(q(x)) + @test evaluate(p, evaluate(q, pi/4)) == cossin1(pi/4) + cossin2 = p(q) + @test evaluate(evaluate(p,q), pi/4) == cossin2(pi/4) + @test evaluate(p, q) == cossin2 + @test p(q)() == evaluate(evaluate(p, q)) + @test evaluate(p, q) == p(q) + @test evaluate(q, p) == q(p) + cs = x->cos(sin(x)) + csdiff = (cs(t17)-cossin2(t17)).(-2:0.1:2) + @test norm(csdiff, 1) < 5e-15 + a = [p,q] + @test a(0.1) == evaluate.([p,q],0.1) + @test a.(0.1) == a(0.1) + @test a.() == evaluate.([p, q]) + @test a.() == [p(), q()] + @test a.() == a() + vr = rand(2) + @test p.(vr) == evaluate.(p, vr) + Mr = rand(3,3,3) + @test p.(Mr) == evaluate.(p, Mr) + mytaylor1 = Taylor1(rand(20)) + vr = rand(5) + @test p(vr) == p.(vr) + @test p(vr) == evaluate.(p,vr) + @test p(Mr) == p.(Mr) + @test p(Mr) == evaluate.(p,Mr) + taylor_a = Taylor1(Int64,10) + taylor_x = exp(Taylor1(Float64,13)) + @test taylor_x(taylor_a) == evaluate(taylor_x, taylor_a) + A_T1 = [t 2t 3t; 4t 5t 6t ] + @test evaluate(A_T1,1.0) == [1.0 2.0 3.0; 4.0 5.0 6.0] + @test evaluate(A_T1,1.0) == A_T1(1.0) + @test evaluate(A_T1) == A_T1() + @test A_T1(tsquare) == [tsquare 2tsquare 3tsquare; 4tsquare 5tsquare 6tsquare] + @test view(A_T1, :, :)(1.0) == A_T1(1.0) + @test view(A_T1, :, 1)(1.0) == A_T1[:,1](1.0) @test sin(asin(tsquare)) == tsquare @test tan(atan(tsquare)) == tsquare @@ -182,62 +272,107 @@ using Base.Test ut = 1.0*t tt = zero(ut) + TaylorSeries.one!(tt, ut, 0) + @test tt[0] == 1.0 + TaylorSeries.one!(tt, ut, 1) + @test tt[1] == 0.0 + TaylorSeries.abs!(tt, 1.0+ut, 0) + @test tt[0] == 1.0 TaylorSeries.add!(tt, ut, ut, 1) - @test tt[2] == 2.0 + @test tt[1] == 2.0 + TaylorSeries.add!(tt, -3.0, 0) + @test tt[0] == -3.0 + TaylorSeries.add!(tt, -3.0, 1) + @test tt[1] == 0.0 TaylorSeries.subst!(tt, ut, ut, 1) - @test tt[2] == 0.0 + @test tt[1] == 0.0 + TaylorSeries.subst!(tt, -3.0, 0) + @test tt[0] == 3.0 + TaylorSeries.subst!(tt, -2.5, 1) + @test tt[1] == 0.0 iind, cind = TaylorSeries.divfactorization(ut, ut) @test iind == 1 @test cind == 1.0 TaylorSeries.div!(tt, ut, ut, 0, iind) - @test tt[1] == cind + @test tt[0] == cind TaylorSeries.div!(tt, 1+ut, 1+ut, 0) + @test tt[0] == 1.0 + TaylorSeries.div!(tt, 1, 1+ut, 0) + @test tt[0] == 1.0 + TaylorSeries.pow!(tt, 1.0+t, 1.5, 0, 0) + @test tt[0] == 1.0 + TaylorSeries.pow!(tt, 1.0+t, 1.5, 0) + @test tt[0] == 1.0 + TaylorSeries.pow!(tt, 1.0+t, 0.5, 1) + @test tt[1] == 0.5 + TaylorSeries.pow!(tt, 1.0+t, 0, 0) + @test tt[0] == 1.0 + TaylorSeries.pow!(tt, 1.0+t, 1, 1) @test tt[1] == 1.0 - TaylorSeries.pow!(tt, 1+t, 1.5, 0, 0) - @test tt[1] == 1.0 - TaylorSeries.pow!(tt, 1+t, 1.5, 0) - @test tt[1] == 1.0 - TaylorSeries.sqrt!(tt, 1+t, 0, 0) - @test tt[1] == 1.0 - TaylorSeries.sqrt!(tt, 1+t, 0) - @test tt[1] == 1.0 - TaylorSeries.exp!(tt, t, 0) - @test tt[1] == exp(t[1]) + tt = zero(ut) + TaylorSeries.pow!(tt, 1.0+t, 2, 0) + @test tt[0] == 1.0 + TaylorSeries.pow!(tt, 1.0+t, 2, 1) + @test tt[1] == 2.0 + TaylorSeries.pow!(tt, 1.0+t, 2, 2) + @test tt[2] == 1.0 + TaylorSeries.sqrt!(tt, 1.0+t, 0, 0) + @test tt[0] == 1.0 + TaylorSeries.sqrt!(tt, 1.0+t, 0) + @test tt[0] == 1.0 + TaylorSeries.exp!(tt, 1.0*t, 0) + @test tt[0] == exp(t[0]) TaylorSeries.log!(tt, 1.0+t, 0) - @test tt[1] == 0.0 + @test tt[0] == 0.0 ct = zero(ut) - TaylorSeries.sincos!(tt, ct, t, 0) - @test tt[1] == sin(t[1]) - @test ct[1] == cos(t[1]) - TaylorSeries.tan!(tt, t, ct, 0) - @test tt[1] == tan(t[1]) - @test ct[1] == tan(t[1])^2 - TaylorSeries.asin!(tt, t, ct, 0) - @test tt[1] == asin(t[1]) - @test ct[1] == sqrt(1.0-t[1]^2) - TaylorSeries.acos!(tt, t, ct, 0) - @test tt[1] == acos(t[1]) - @test ct[1] == sqrt(1.0-t[1]^2) + TaylorSeries.sincos!(tt, ct, 1.0*t, 0) + @test tt[0] == sin(t[0]) + @test ct[0] == cos(t[0]) + TaylorSeries.tan!(tt, 1.0*t, ct, 0) + @test tt[0] == tan(t[0]) + @test ct[0] == tan(t[0])^2 + TaylorSeries.asin!(tt, 1.0*t, ct, 0) + @test tt[0] == asin(t[0]) + @test ct[0] == sqrt(1.0-t[0]^2) + TaylorSeries.acos!(tt, 1.0*t, ct, 0) + @test tt[0] == acos(t[0]) + @test ct[0] == sqrt(1.0-t[0]^2) TaylorSeries.atan!(tt, ut, ct, 0) - @test tt[1] == atan(t[1]) - @test ct[1] == 1.0+t[1]^2 + @test tt[0] == atan(t[0]) + @test ct[0] == 1.0+t[0]^2 TaylorSeries.sinhcosh!(tt, ct, ut, 0) - @test tt[1] == sinh(t[1]) - @test ct[1] == cosh(t[1]) + @test tt[0] == sinh(t[0]) + @test ct[0] == cosh(t[0]) TaylorSeries.tanh!(tt, ut, ct, 0) - @test tt[1] == tanh(t[1]) - @test ct[1] == tanh(t[1])^2 + @test tt[0] == tanh(t[0]) + @test ct[0] == tanh(t[0])^2 v = [sin(t), exp(-t)] - vv = Vector{Float64}(2) + @compat vv = Vector{Float64}(undef, 2) @test evaluate!(v, zero(Int), vv) == nothing @test vv == [0.0,1.0] @test evaluate(v) == vv @test evaluate(v, complex(0.0,0.2)) == [complex(0.0,sinh(0.2)),complex(cos(0.2),sin(-0.2))] + + @test derivative(exp(ta(1.0)), 0) == exp(ta(1.0)) + expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0))[0:10])) + @test derivative(exp(ta(1.0)), 5) ≈ expected_result_approx atol=eps() rtol=0.0 + expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0pi))[0:12]),15) + @test derivative(exp(ta(1.0pi)), 3) ≈ expected_result_approx atol=eps(16.0) rtol=0.0 + expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0pi))[0:5]),15) + @test derivative(exp(ta(1.0pi)), 10) ≈ expected_result_approx atol=eps(64.0) rtol=0.0 + + + + @test derivative(exp(ta(1.0)), 5)() == exp(1.0) + @test derivative(exp(ta(1.0pi)), 3)() == exp(1.0pi) + @test isapprox(derivative(exp(ta(1.0pi)), 10)() , exp(1.0pi) ) + @test derivative(5, exp(ta(1.0))) == exp(1.0) @test derivative(3, exp(ta(1.0pi))) == exp(1.0pi) @test isapprox(derivative(10, exp(ta(1.0pi))) , exp(1.0pi) ) + @test integrate(derivative(exp(t)),1) == exp(t) @test integrate(cos(t)) == sin(t) @@ -245,24 +380,100 @@ using Base.Test @test norm((inverse(exp(t)-1) - log(1+t)).coeffs) < 2tol1 cfs = [(-n)^(n-1)/gamma(n+1) for n = 1:15] - @test norm(inverse(t*exp(t))[2:end]./cfs-1) < 4tol1 + @test norm(inverse(t*exp(t))[1:end]./cfs .- 1) < 4tol1 @test_throws ArgumentError Taylor1([1,2,3], -2) @test_throws ArgumentError abs(ta(big(0))) @test_throws ArgumentError 1/t @test_throws ArgumentError zt/zt @test_throws ArgumentError t^1.5 - @test_throws DomainError t^(-2) + if VERSION < v"0.7.0-DEV" + @test_throws DomainError t^(-2) + else + @test_throws ArgumentError t^(-2) + end @test_throws ArgumentError sqrt(t) @test_throws ArgumentError log(t) @test_throws ArgumentError cos(t)/sin(t) @test_throws AssertionError derivative(30, exp(ta(1.0pi))) @test_throws ArgumentError inverse(exp(t)) + @test_throws ArgumentError abs(t) + displayBigO(false) + @test string(ta(-3)) == " - 3 + 1 t " + @test string(ta(0)^3-3) == " - 3 + 1 t³ " + @test TaylorSeries.pretty_print(ta(3im)) == " ( 3 im ) + ( 1 ) t " + @test string(Taylor1([1,2,3,4,5], 2)) == string(Taylor1([1,2,3])) + displayBigO(true) @test string(ta(-3)) == " - 3 + 1 t + 𝒪(t¹⁶)" @test string(ta(0)^3-3) == " - 3 + 1 t³ + 𝒪(t¹⁶)" @test TaylorSeries.pretty_print(ta(3im)) == " ( 3 im ) + ( 1 ) t + 𝒪(t¹⁶)" @test string(Taylor1([1,2,3,4,5], 2)) == string(Taylor1([1,2,3])) + + a = collect(1:12) + t_a = Taylor1(a,15) + t_C = complex(3.0,4.0) * t_a + rnd = rand(10) + @test typeof( norm(Taylor1(rnd)) ) == Float64 + @test norm(Taylor1(rnd)) > 0 + @test norm(t_a) == norm(a) + @test norm(Taylor1(a,15),3) == sum((a.^3))^(1/3) + @test norm(t_a,Inf) == 12 + @test norm(t_C) == norm(complex(3.0,4.0)*a) + + @test TaylorSeries.rtoldefault(Taylor1{Int64}) == 0 + @test TaylorSeries.rtoldefault(Taylor1{Float64}) == sqrt(eps(Float64)) + @test TaylorSeries.rtoldefault(Taylor1{BigFloat}) == sqrt(eps(BigFloat)) + @test TaylorSeries.real(Taylor1{Float64}) == Taylor1{Float64} + @test TaylorSeries.real(Taylor1{Complex{Float64}}) == Taylor1{Float64} + @test isfinite(t_C) + @test isfinite(t_a) + @test !isfinite( Taylor1([0, Inf]) ) + @test !isfinite( Taylor1([NaN, 0]) ) + b = convert(Vector{Float64}, a) + b[3] += eps(10.0) + b[5] -= eps(10.0) + t_b = Taylor1(b,15) + t_C2 = t_C+eps(100.0) + t_C3 = t_C+eps(100.0)*im + @test isapprox(t_C, t_C) + @test t_a ≈ t_a + @test t_a ≈ t_b + @test t_C ≈ t_C2 + @test t_C ≈ t_C3 + @test t_C3 ≈ t_C2 + t = Taylor1(25) + p = sin(t) + q = sin(t+eps()) + @test t ≈ t + @test t ≈ t+sqrt(eps()) + @test isapprox(p, q, atol=eps()) + + t = Taylor1(35) + @test Taylor1([180.0, rad2deg(1.0)], 35) == rad2deg(pi+t) + @test sin(pi/2+deg2rad(1.0)t) == sin(deg2rad(90+t)) + a = Taylor1(rand(10)) + b = Taylor1(rand(10)) + c = deepcopy(a) + TaylorSeries.deg2rad!(b, a, 0) + @test a == c + @test a[0]*(pi/180) == b[0] + TaylorSeries.deg2rad!.(b, a, [0,1,2]) + @test a == c + for i in 0:2 + @test a[i]*(pi/180) == b[i] + end + a = Taylor1(rand(10)) + b = Taylor1(rand(10)) + c = deepcopy(a) + TaylorSeries.rad2deg!(b, a, 0) + @test a == c + @test a[0]*(180/pi) == b[0] + TaylorSeries.rad2deg!.(b, a, [0,1,2]) + @test a == c + for i in 0:2 + @test a[i]*(180/pi) == b[i] + end end @testset "Matrix multiplication for Taylor1" begin @@ -281,10 +492,10 @@ end B = Taylor1{Float64}[Taylor1(collect(B1[i,1:i]),i) for i=1:n1] Y = Taylor1{Float64}[Taylor1(collect(Y1[k,1:k]),k) for k=1:k1] Bcopy = deepcopy(B) - A_mul_B!(Y,A,B) + @compat mul!(Y,A,B) # do we get the same result when using the `A*B` form? - @test A*B==Y + @test A*B≈Y # Y should be extended after the multilpication @test reduce(&, [y1.order for y1 in Y] .== Y[1].order) # B should be unchanged @@ -292,7 +503,7 @@ end # is the result compatible with the matrix multiplication? We # only check the zeroth order of the Taylor series. - y1=sum(Y)[1] + y1=sum(Y)[0] Y=A*B1[:,1] y2=sum(Y) @@ -300,9 +511,9 @@ end # multiplication and the specialized version @test abs(y1-y2) < n1*(eps(y1)+eps(y2)) - @test_throws DimensionMismatch A_mul_B!(Y,A[:,1:end-1],B) - @test_throws DimensionMismatch A_mul_B!(Y,A[1:end-1,:],B) - @test_throws DimensionMismatch A_mul_B!(Y,A,B[1:end-1]) - @test_throws DimensionMismatch A_mul_B!(Y[1:end-1],A,B) + @compat @test_throws DimensionMismatch mul!(Y,A[:,1:end-1],B) + @compat @test_throws DimensionMismatch mul!(Y,A[1:end-1,:],B) + @compat @test_throws DimensionMismatch mul!(Y,A,B[1:end-1]) + @compat @test_throws DimensionMismatch mul!(Y[1:end-1],A,B) end end diff --git a/test/runtests.jl b/test/runtests.jl index 05a68071..4abbb861 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,8 +6,10 @@ testfiles = ( "onevariable.jl", "manyvariables.jl", "mixtures.jl", + "mutatingfuncts.jl", "identities_Euler.jl", - "fateman40.jl" + "fateman40.jl", + "intervals.jl" ) for file in testfiles