Library

TaylorSeries

A Julia package for Taylor expansions in one or more independent variables.

The basic constructors are Taylor1 and TaylorN; see also HomogeneousPolynomial.

Taylor1{T<:Number} <: AbstractSeries{T}

DataType for polynomial expansions in one independent variable.

Fields:

• coeffs :: FixedSizeVectorDefault{T} Expansion coefficients; the $i$-th component is the coefficient of degree $i-1$ of the expansion.

Note that Taylor1 variables are callable. For more information, see evaluate.

HomogeneousPolynomial{T<:Number} <: AbstractSeries{T}

DataType for homogeneous 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 of the independent variables are specified by coeff_table[order+1][i].

• order :: Int order (degree) of the homogeneous polynomial.

Note that HomogeneousPolynomial variables are callable. For more information, see evaluate.

TaylorN{T<:Number} <: AbstractSeries{T}

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$.

• space :: JetSpace multivariate algebra associated with the polynomial.

Use order to obtain the maximum order of the polynomial expansion.

Note that TaylorN variables are callable. For more information, see evaluate.

AbstractSeries{T<:Number} <: Number

Parameterized abstract type for Taylor1, HomogeneousPolynomial and TaylorN.

JetSpace

Explicit multivariate Taylor algebra. A space owns the truncation order, variable metadata, and the lookup tables used by HomogeneousPolynomial and TaylorN arithmetic.

Fields

• order: maximum total degree.
• num_vars: number of variables.
• variable_names: variable names used for printing and compatibility APIs.
• variable_symbols: variable symbols corresponding to variable_names.
• coeff_table: exponent vectors grouped by homogeneous degree.
• index_table: hashed exponent-vector labels grouped by homogeneous degree.
• size_table: number of monomials for each homogeneous degree.
• pos_table: maps hashed exponent labels to coefficient positions.
• mul_table: product-table cache indexed by the positive degrees of the left and right homogeneous factors. Only products that fit within the truncation order are stored: row degree_a contains entries for degree_b = 1:(order-degree_a). Degree-zero products are scalar shortcuts and are not cached. Entries start as empty placeholders; input_positions is filled on first product-table use, while output_offsets and output_pairs are filled only if an output-grouped multiplication routine requests them.
• mul_table_lock: lock guarding lazy initialization of mul_table entries.

Taylor1([T::Type=Float64], order::Int)

Shortcut to define the independent variable of a Taylor1{T} polynomial of given order. The default type for T is Float64.

julia> Taylor1(16)
 1.0 t + 𝒪(t¹⁷)

julia> Taylor1(Rational{Int}, 4)
 1//1 t + 𝒪(t⁵)

HomogeneousPolynomial([T::Type=Float64], nv::Int])

Shortcut to define the nv-th independent HomogeneousPolynomial{T}. The default type for T is Float64.

julia> HomogeneousPolynomial(1)
 1.0 x₁

julia> HomogeneousPolynomial(Rational{Int}, 2)
 1//1 x₂

TaylorN([T::Type=Float64], nv::Int; [order::Int=order()])

Shortcut to define the nv-th independent TaylorN{T} variable as a polynomial. The order is defined through the keyword parameter order, whose default corresponds to order(). The default of type for T is Float64.

julia> TaylorN(1)
 1.0 x₁ + 𝒪(‖x‖⁷)

julia> TaylorN(Rational{Int},2)
 1//1 x₂ + 𝒪(‖x‖⁷)

variables!([T::Type], names::String; [order=order(), numvars=-1, nowarn=false])

Return a TaylorN{T} vector with each entry representing an independent variable. names defines the output for each variable (separated by a space). The default type T is Float64, and the default for order is the current default-space order.

This function updates the compatibility default JetSpace. Set nowarn=true to suppress the warning emitted when the default space changes.

If numvars is not specified, it is inferred from names. If only one variable name is defined and numvars>1, it uses this name with subscripts for the different variables.

julia> variables!(Int, "x y z", order=4)
3-element Array{TaylorSeries.TaylorN{Int},1}:
  1 x + 𝒪(‖x‖⁵)
  1 y + 𝒪(‖x‖⁵)
  1 z + 𝒪(‖x‖⁵)

julia> variables!("α", numvars=2)
2-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 α₁ + 𝒪(‖x‖⁵)
  1.0 α₂ + 𝒪(‖x‖⁵)

julia> variables!("x", order=6, numvars=2)
2-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 x₁ + 𝒪(‖x‖⁷)
  1.0 x₂ + 𝒪(‖x‖⁷)

variables([T::Type=Float64], [order::Int=order()])
variables([T::Type=Float64], space::JetSpace; order=order(space))

Return a TaylorN{T} vector with each entry representing an independent variable. Without an explicit space, variables uses the current compatibility default space. The order can be explicitly set without changing the default-space values for num_vars or variable_names.

With an explicit space, variables(space) returns the independent variables for that space. Omitting T defaults to Float64.

order()
order(a)
order(space::JetSpace)

Return the maximum expansion order of the current default multivariate algebra, of a Taylor series or homogeneous polynomial a, or of an explicit JetSpace.

show_params_TaylorN()

Display the current parameters for TaylorN and HomogeneousPolynomial types.

show_monomials(ord::Int) --> nothing
show_monomials(space::JetSpace, ord::Int) --> nothing

List the indices and corresponding monomials of a HomogeneousPolynomial of degree ord, using either the default JetSpace or the given explicit space.

getcoeff(a, n)

Return the coefficient of order n::Int of a a::Taylor1 polynomial; the constant term corresponds to n=0.

getcoeff(a, v)

Return the coefficient of a::HomogeneousPolynomial, specified by v, which is a tuple (or vector) with the indices of the specific monomial.

getcoeff(a, v)

Return the coefficient of a::TaylorN, specified by v, which is a tuple (or vector) with the indices of the specific monomial.

evaluate(a, [dx])

Evaluate a Taylor1 polynomial using Horner's rule (hand coded). If dx is omitted, 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).

evaluate(x, δt)

Evaluates each element of x::AbstractArray{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).

evaluate(a::Taylor1, x::Taylor1)

Substitute x::Taylor1 as independent variable in a a::Taylor1 polynomial. Note that the syntax a(x) is equivalent to evaluate(a, x).

evaluate(a, [vals])

Evaluate a HomogeneousPolynomial polynomial at vals. If vals is omitted, it's evaluated at zero. Note that the syntax a(vals) is equivalent to evaluate(a, vals); and a() is equivalent to evaluate(a).

evaluate(a, [vals]; sorting::Bool=true)

Evaluate the TaylorN polynomial a at vals. If vals is omitted, it's evaluated at zero. The keyword parameter sorting can be used to avoid sorting (in increasing order by abs2) the terms that are added.

Note that the syntax a(vals) is equivalent to evaluate(a, vals); and a() is equivalent to evaluate(a); use a(b::Bool, x) corresponds to evaluate(a, x, sorting=b).

evaluate!(x, δt, dest)
evaluate!(x, vals, dest; sorting=...)
evaluate!(x, δt, dest, aux)
evaluate!(x, vals, dest, valscache, aux; sorting=false)

Evaluate a polynomial, or an array of polynomials, and write the result into dest, a pre-allocated destination. All these methods return nothing.

For TaylorN evaluation, sorting=true sorts the contributions by magnitude before adding them, to reduce rounding error. sorting=false skips sorting. The default, given by _defaultsorting, is true when both the coefficients and evaluation values are usual numbers (other than Taylor series or intervals), and false otherwise. Passing pre-allocated auxiliaries does not disable sorting; choose sorting=false to avoid the temporary result used by sorted evaluation.

The methods for a single polynomial are also used by the corresponding array methods. The allocating evaluate methods do not all call evaluate!; some call evaluate separately for each array element.

Examples

Evaluate two polynomials at a number:

julia> using TaylorSeries

julia> t = Taylor1(3);

julia> polys = [1 + 2t, 3 + 4t];

julia> dest = zeros(2);

julia> evaluate!(polys, 0.5, dest);

julia> dest == [2.0, 5.0]
true

Evaluate a multivariable polynomial at Taylor series. Allocate the destination and auxiliaries once; they can be reused in subsequent calls. The docstrings below explain which orders and JetSpaces are allowed, and which objects must be kept separate.

julia> using TaylorSeries

julia> space = JetSpace(order=2, variables=[:x, :y]);

julia> x, y = variables(space);

julia> vals = (x + 1, y + 2);

julia> dest = zero(x); valscache = [zero(x), zero(y)]; aux = zero(x);

julia> evaluate!(x + y^2, vals, dest, valscache, aux);

julia> dest == (x + 1) + (y + 2)^2
true

evaluate!(a, δt, dest, aux)
evaluate!(x, δt, dest, aux)

Evaluate a Taylor1 polynomial, or an array of them, at a Taylor series δt. Write the result into dest and reuse the pre-allocated auxiliary aux for the calculations. The method does not create aux on each call.

The coefficients of a, δt, dest, and aux must have the same series type, including their coefficient type, and the same JetSpace for TaylorN. Keep δt, dest, and aux separate from one another. Neither dest nor aux may share coefficient storage with a coefficient of a. dest and aux must have the same order, no greater than the order of δt. Each coefficient of a must have at least that order, except that order-zero coefficients are treated as exact constants.

evaluate!(a, vals, dest, valscache, aux; sorting=false)
evaluate!(a_array, vals, dest_array, valscache, aux; sorting=false)

Evaluate a TaylorN polynomial, or an array of them, at the Taylor series in vals and write into dest. valscache provides one pre-allocated auxiliary per evaluation value, and aux is another pre-allocated auxiliary used for the arithmetic. Both are overwritten during evaluation. Every evaluation value and auxiliary must have the same JetSpace and order as dest; the source polynomial must share that JetSpace but may have a different order. Keep each entry of valscache and aux separate from the inputs, destination, and other auxiliaries. The destination must also be separate from the inputs.

With sorting=false, these methods reuse the supplied valscache and aux. With sorting=true, each polynomial is evaluated with its contributions sorted by magnitude before addition, and the result is copied into dest. This sorted evaluation may allocate an intermediate result.

evaluate!(a, vals, dest, valscache, aux; sorting=false)

Evaluate each TaylorN polynomial in a, writing the results into dest. Reuse the same pre-allocated auxiliaries valscache and aux for each polynomial in turn. Their requirements and the sorting behavior are the same as for evaluating a single polynomial.

Check the shared evaluation values and pre-allocated auxiliaries once per call. Then check every source polynomial and destination before evaluating any of them.

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.

update!(a, x0)

Takes a <: Union{Taylo1,TaylorN} and expands it around the coordinate x0.

differentiate(a)

Return the Taylor1 polynomial of the differential of a::Taylor1. The order of the result is order(a)-1.

The function derivative is an exact synonym of differentiate.

differentiate(a, n)

Compute recursively the Taylor1 polynomial of the n-th derivative of a::Taylor1. The order of the result is order(a)-n.

differentiate(n, a)

Return the value of the n-th differentiate of the polynomial a.

differentiate(a, r)

Partial differentiation of a::HomogeneousPolynomial series with respect to the r-th variable.

differentiate(a, r)

Partial differentiation of a::TaylorN series with respect to the r-th variable. The r-th variable may be also specified through its symbol.

differentiate(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.

differentiate(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.

derivative

An exact synonym of differentiate.

integrate(a, [x])

Return the integral of a::Taylor1. The constant of integration (0-th order coefficient) is set to x, which is zero if omitted. Note that the order of the result is order(a)+1.

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 order(), a zero HomogeneousPolynomial of 0-th order is returned.

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 omitted, it is taken as zero.

    gradient(f)
    ∇(f)

Compute the gradient of the polynomial f::TaylorN.

    jacobian(vf)
    jacobian(vf, [vals])

Compute the jacobian matrix of vf, a vector of TaylorN polynomials, evaluated at the vector vals. If vals is omitted, it is evaluated at zero.

    jacobian!(jac, vf)
    jacobian!(jac, vf, [vals])

Compute the jacobian matrix of vf, a vector of TaylorN polynomials evaluated at the vector vals, and write results to jac. If vals is omitted, it is evaluated at zero.

    jacobianmatrix(vf)

Compute the jacobian matrix of vf, a vector of TaylorN polynomials.

    hessian(f)
    hessian(f, [vals])

Return the hessian matrix (jacobian of the gradient) of f::TaylorN, evaluated at the vector vals. If vals is omitted, it is evaluated at zero.

    hessian!(hes, f)
    hessian!(hes, f, [vals])

Return the hessian matrix (jacobian of the gradient) of f::TaylorN, evaluated at the vector vals, and write results to hes. If vals is omitted, it is evaluated at zero.

    hessianmatrix(f)

Return the hessian matrix (jacobian of the gradient) of f::TaylorN.

constant_term(a)

Return the constant value (zero order coefficient) for Taylor1 and TaylorN. The fallback behavior is to return a itself if a::Number, or a[1] when a::Vector.

linear_polynomial(a)

Returns the linear part of a as a polynomial (Taylor1 or TaylorN), without the constant term. The fallback behavior is to return a itself.

nonlinear_polynomial(a)

Returns the nonlinear part of a. The fallback behavior is to return zero(a).

inverse(f)

Return the Taylor expansion of $f^{-1}(t)$, of order N = order(f), for f::Taylor1 polynomial, assuming the first coefficient of f is zero. Otherwise, a DomainError is thrown.

The algorithm implements Lagrange inversion at $t=0$ if $f(0)=0$:

\[\begin{equation*} f^{-1}(t) = \sum_{n=1}^{N} \frac{t^n}{n!} \left. \frac{{\rm d}^{n-1}}{{\rm d} z^{n-1}}\left(\frac{z}{f(z)}\right)^n \right\vert_{z=0}. \end{equation*}\]

inverse_map(f)

Return the Taylor expansion of $f^{-1}(t)$, of order N = order(f), for Taylor1 or TaylorN polynomials, assuming the first coefficient of f is zero. Otherwise, a DomainError is thrown.

This method is based in the algorithm by M. Berz, Modern map methods in Particle Beam Physics, Academic Press (1999), Sect 2.3.1. See inverse (for f::Taylor1).

abs(a::Taylor1{T})
abs(a::TaylorN{T})

For T<:Real type returns a if constant_term(a) > 0 and -a if constant_term(a) < 0. For other Number types, it returns sqrt(abs2(a)).

norm(x::AbstractSeries, p::Real)

Returns the p-norm of an x::AbstractSeries, defined by

\[\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.

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.

isless(a::Taylor1{<:Real}, b::Real)
isless(a::TaylorN{<:Real}, b::Real)

Compute isless by comparing the constant_term(a) and b. If they are equal, returns a[nz] < 0, with nz the first non-zero coefficient after the constant term. This defines a total order.

For many variables, the ordering includes a lexicographical convention in order to be total. We have opted for the simplest one: the larger variable appears before for the TaylorN variables are defined (e.g., through variables!).

For nested Taylor1{Taylor1{...}}s, the ordering is established by which one is a constant_term of the other.

Refs:

• M. Berz, AIP Conference Proceedings 177, 275 (1988); https://doi.org/10.1063/1.37800
• M. Berz, "Automatic Differentiation as Nonarchimedean Analysis", Computer Arithmetic and Enclosure Methods, (1992), Elsevier, 439-450.

isless(a::Taylor1{T}, b::Taylor1{T})
isless(a::TaylorN{T}, b::Taylor1{T})

Returns isless(a - b, zero(T)).

isfinite(x::AbstractSeries) -> Bool

Test whether the coefficients of the polynomial x are finite.

displayBigO(d::Bool) --> nothing

Set/unset displaying of the big 𝒪 notation in the output of Taylor1 and TaylorN polynomials. The initial value is true.

use_Base_show(d::Bool) --> nothing

Use Base.show_default method (default show method in Base), or a custom display. The initial value is false, so customized display is used.

set_taylor1_varname(var::String)
set_taylor1_varname(numvars::Int, names::String)
set_taylor1_varname(numvars::Int, names::Vector{String})

Change the displayed variable(s) for Taylor1 objects; the names of variables cannot include spaces

space(a::Union{HomogeneousPolynomial,TaylorN})

Return the JetSpace associated with the multivariate Taylor object a.

ParamsTaylor1

DataType holding the current variable name(s) for Taylor1. The parameters can be changed using set_taylor1_varname

HomogeneousProductTable

Precomputed input-pair and output-major schedules for multiplying two homogeneous polynomials in a single JetSpace.

For multiplying homogeneous polynomials, say, c = a * b, each coefficient of a multiplies each coefficient of b. So the code walks pairs like::

a[1] * b[1]
a[1] * b[2]
a[1] * b[3]
a[2] * b[1]
...

For each pair of indices, the field input_positions[pair] tells us which output coefficient c[pos] receives that product (in general, many pairs contribute to the same output coefficient). The input-pair traversing order schedule maps each pair of input coefficients to the output coefficient it contributes to. The output-major schedule stores the same products grouped by output coefficient, so dense kernels (e.g., _mul_output_major_unchecked!) can accumulate one output coefficient at a time.

At initialization, only the input-pair table input_positions is filled. Later, only if a kernel needs the output-major layout, the fields output_offsets and output_pairs are filled, i.e., output-major data is filled lazily only when a kernel requests it. Currently, this lazy loading is realized via a mutable struct, so that fields can be assigned when appropriate, but it is not strictly necessary.

Fields

• input_positions: output coefficient position for each input coefficient pair.
• output_offsets: start positions of each output coefficient's group in output_pairs.
• output_pairs: input-pair identifiers grouped by output coefficient.
• num_right: number of coefficients in the right input polynomial.

_InternalMutFuncs

Contains parameters and expressions that allow a simple programmatic construction for calling the internal mutating functions.

generate_tables(num_vars, order)

Return the lookup tables coeff_table, index_table, size_table and pos_table used by a JetSpace.

Hash tables

coeff_table :: Array{Array{Array{Int,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$.

index_table :: Array{Array{Int,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$.

size_table :: Array{Int,1}

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{Int,Int},1}

The $i+1$-th component maps the hash index to the (lexicographic) position of the corresponding monomial in coeffs_table.

Return the input-position schedule for multiplying two homogeneous degrees.

Initialize and return the output-major product schedule for two positive degrees.

Return empty valid positive-degree product-table placeholders for lazy initialization.

Return the cached product table for two positive degrees, initializing it if needed.

Initialize and cache the input-position product table for two positive degrees.

generate_index_vectors(num_vars, degree)

Return a vector of index vectors with num_vars (number of variables) and degree.

in_base(order, v)

Convert vector v of non-negative integers to base oorder, where oorder is the next odd integer of order.

make_inverse_dict(v)

Return a Dict with the enumeration of v: the elements of v point to the corresponding index.

It is used to construct pos_table from index_table.

Set the compatibility default JetSpace.

The global default_space binding is a constant Ref, but its contents are mutable. Module loading initializes default_space[]; later compatibility calls such as variables! replace default_space[] with space. Existing TaylorN and HomogeneousPolynomial objects keep their original spaces, while future default-space constructors use the new default algebra.

Use nowarn=true to suppress the warning emitted when replacing the default space.

_coeffsHP(x::T, order::Int) where {T<:Number}
_coeffsHP(coeffs::AbstractArray{T,1}, order::Int) where {T<:Number}
_coeffsHP(space::JetSpace, x::T, order::Int) where {T<:Number}
_coeffsHP(space::JetSpace, coeffs::AbstractArray{T,1}, order::Int) where {T<:Number}

Returns a FixedSizeVectorDefault of size space.size_table[order+1] to be used in the construction of a HomogeneousPolynomial of order order. The returned vector has the first entries of coeffs, and then is filled with zeros.

_coeffsTN(v::AbstractArray{T,1}, order::Int) where {T<:Number}

Returns a FixedSizeVectorDefault{HomogeneousPolynomial{T}} of size order+1, to be used in the construction of a TaylorN{T} of order order. The returned vector has the entries of v at the proper location according to their order, and otherwise it is filled with the corresponding zeros.

_check_same_space(space_a::JetSpace, space_b::JetSpace)
_check_same_space(a::Taylor1, b::Taylor1)
_check_same_space(a::Taylor1, b::Taylor1, c::Taylor1)
_check_same_space(a::Union{HomogeneousPolynomial,TaylorN},
    b::Union{HomogeneousPolynomial,TaylorN})
_check_same_space(a::Union{HomogeneousPolynomial,TaylorN},
    b::Union{HomogeneousPolynomial,TaylorN},
    c::Union{HomogeneousPolynomial,TaylorN})
_check_same_space(space::JetSpace, v::AbstractVector{<:HomogeneousPolynomial})

Throw an ArgumentError unless all arguments belong to the same JetSpace by object identity.

numtype(a::AbstractSeries)

Returns the type of the elements of the coefficients of a.

mul!(c, a, b, k::Int) --> nothing

Update the k-th expansion coefficient c[k] of c = a * b, where all c, a, and b are either Taylor1 or TaylorN. Note that for TaylorN the result of a * b is accumulated in c[k].

The coefficients are given by

\[c_k = \sum_{j=0}^k a_j b_{k-j}.\]

mul!(c, a, b) --> nothing

Accumulates in c the result of a*b with minimum allocation. Arguments c, a and b are HomogeneousPolynomial.

mul_scalar!(c, scalar, a, b) --> nothing

Accumulates in c the result of scalar*a*b with minimum allocation. Arguments c, a and b are HomogeneousPolynomial; scalar is a NumberNotSeries.

mul!(Y, A, B)

Multiply A*B and save the result in Y.

div!(c, a, b, k::Int)

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

\[c_k = \frac{1}{b_0} \big(a_k - \sum_{j=0}^{k-1} c_j b_{k-j}\big).\]

For Taylor1 polynomials, a similar formula is implemented which exploits k_0, the order of the first non-zero coefficient of a.

pow!(c, a, aux, r::Real, k::Int)

Update the k-th expansion coefficient c[k] of c = a^r, for both c, a and aux either Taylor1 or TaylorN.

The coefficients are given by

\[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, a similar formula is implemented which exploits k_0, the order of the first non-zero coefficient of a.

square(a::AbstractSeries) --> typeof(a)

Return a^2; see TaylorSeries.sqr!.

sqr!(c, a, aux, k::Int) --> nothing

Update the k-th expansion coefficient c[k] of c = a^2, for both c and a either Taylor1{T} or TaylorN{T}; aux::T is an auxiliary.

The coefficients are given by

\[\begin{aligned} c_k &= 2 \sum_{j=0}^{(k-1)/2} a_{k-j} a_j, \text{ if $k$ is odd,} \\ c_k &= 2 \sum_{j=0}^{(k-2)/2} a_{k-j} a_j + (a_{k/2})^2, \text{ if $k$ is even.} \end{aligned}\]

accsqr!(c, a)

Returns c += a*a with no allocation; all parameters are HomogeneousPolynomial.

sqrt!(c, a, aux, k::Int, k0::Int=0)

Compute the k-th expansion coefficient c[k] of c = sqrt(a) for bothc and a either Taylor1 or TaylorN.

The coefficients are given by

\[\begin{aligned} 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{aligned}\]

For Taylor1 polynomials, k0 is the order of the first non-zero coefficient, which must be even.

exp!(c, a, k) --> nothing

Update the k-th expansion coefficient c[k+1] of c = exp(a) for both c and a either Taylor1 or TaylorN.

The coefficients are given by

\[\begin{equation*} c_k = \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} c_j. \end{equation*}\]

log!(c, a, k) --> nothing

Update the k-th expansion coefficient c[k+1] of c = log(a) for both c and a either Taylor1 or TaylorN.

The coefficients are given by

\[\begin{equation*} c_k = \frac{1}{a_0} \big(a_k - \frac{1}{k} \sum_{j=0}^{k-1} j a_{k-j} c_j \big). \end{equation*}\]

sincos!(s, c, a, k) --> nothing

Update the k-th expansion coefficients s[k+1] and c[k+1] of s = sin(a) and c = cos(a) simultaneously, for s, c and a either Taylor1 or TaylorN.

The coefficients are given by

\[\begin{aligned} s_k &= \frac{1}{k}\sum_{j=0}^{k-1} (k-j) a_{k-j} c_j ,\\ c_k &= -\frac{1}{k}\sum_{j=0}^{k-1} (k-j) a_{k-j} s_j. \end{aligned}\]

tan!(c, a, p, k::Int) --> nothing

Update the k-th expansion coefficients c[k+1] of c = tan(a), for c and a either Taylor1 or TaylorN; p = c^2 and is passed as an argument for efficiency.

The coefficients are given by

\[\begin{equation*} c_k = a_k + \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} p_j. \end{equation*}\]

asin!(c, a, r, k)

Update the k-th expansion coefficients c[k+1] of c = asin(a), for c and a either Taylor1 or TaylorN; r = sqrt(1-c^2) and is passed as an argument for efficiency.

\[\begin{equation*} c_k = \frac{1}{ \sqrt{r_0} } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big). \end{equation*}\]

acos!(c, a, r, k)

Update the k-th expansion coefficients c[k+1] of c = acos(a), for c and a either Taylor1 or TaylorN; r = sqrt(1-c^2) and is passed as an argument for efficiency.

\[\begin{equation*} c_k = - \frac{1}{ r_0 } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big). \end{equation*}\]

atan!(c, a, r, k)

Update the k-th expansion coefficients c[k+1] of c = atan(a), for c and a either Taylor1 or TaylorN; r = 1+a^2 and is passed as an argument for efficiency.

\[\begin{equation*} c_k = \frac{1}{r_0}\big(a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j\big). \end{equation*}\]

sinhcosh!(s, c, a, k)

Update the k-th expansion coefficients s[k+1] and c[k+1] of s = sinh(a) and c = cosh(a) simultaneously, for s, c and a either Taylor1 or TaylorN.

The coefficients are given by

\[\begin{aligned} s_k = \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} c_j, \\ c_k = \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} s_j. \end{aligned}\]

tanh!(c, a, p, k)

Update the k-th expansion coefficients c[k+1] of c = tanh(a), for c and a either Taylor1 or TaylorN; p = a^2 and is passed as an argument for efficiency.

\[\begin{equation*} c_k = a_k - \frac{1}{k} \sum_{j=0}^{k-1} (k-j) a_{k-j} p_j. \end{equation*}\]

asinh!(c, a, r, k)

Update the k-th expansion coefficients c[k+1] of c = asinh(a), for c and a either Taylor1 or TaylorN; r = sqrt(1-c^2) and is passed as an argument for efficiency.

\[\begin{equation*} c_k = \frac{1}{ \sqrt{r_0} } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big). \end{equation*}\]

acosh!(c, a, r, k)

Update the k-th expansion coefficients c[k+1] of c = acosh(a), for c and a either Taylor1 or TaylorN; r = sqrt(c^2-1) and is passed as an argument for efficiency.

\[\begin{equation*} c_k = \frac{1}{ r_0 } \big( a_k - \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j \big). \end{equation*}\]

atanh!(c, a, r, k)

Update the k-th expansion coefficients c[k+1] of c = atanh(a), for c and a either Taylor1 or TaylorN; r = 1-a^2 and is passed as an argument for efficiency.

\[\begin{equation*} c_k = \frac{1}{r_0}\big(a_k + \frac{1}{k} \sum_{j=1}^{k-1} j r_{k-j} c_j\big). \end{equation*}\]

differentiate!(res, a) --> nothing

In-place version of differentiate. Compute the Taylor1 polynomial of the differential of a::Taylor1 and return it as res (order of res remains unchanged).

differentiate!(p, a, k) --> nothing

Update in-place the k-th expansion coefficient p[k] of p = differentiate(a) for both p and a Taylor1.

The coefficients are given by

\[p_k = (k+1) a_{k+1}.\]

constant_term!(a, c)

Update the constant term (zero order coefficient) of a to c, leaving higher order coefficients unchanged.

_isthinzero(x)

Generic wrapper to function iszero, which allows using the correct function for Intervals

_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 and _dict_binary_calls. The call contains the prefix TaylorSeries..

_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, an Expr defining its calling pattern is added as the last entry of the vector.

_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.

_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.

_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.

_populate_dicts!()

Function that populates the internal dictionaries _dict_unary_calls and _dict_binary_calls

_evaluate(a::Taylor1, dx::NumberNotSeries)

Evaluate a at a number dx using Horner's rule. Neither the coefficients nor dx are Taylor series. Both evaluate and the array form of evaluate! use this method; it does not create temporary arrays.

_evaluate(a::HomogeneousPolynomial, vals)

Evaluate one homogeneous polynomial at vals, computing each term directly without first collecting powers or terms in arrays. Supply one value per variable. If the values are TaylorN series, their JetSpaces must match the polynomial's.

_evaluate(a::TaylorN, vals::NTuple, ::Val{true})
_evaluate(a::TaylorN, vals::NTuple, ::Val{false})

Evaluate a at vals and return the sum of its contributions. Val(true) sorts the contributions from smallest to largest magnitude before adding them, to reduce rounding error. Val(false) skips this sorting step.

The third argument is needed because _evaluate(a::TaylorN, vals::NTuple) returns an array with one contribution per polynomial degree, rather than their sum. The sorted method needs those separate values so it can reorder them before adding them. Using Val(true) or Val(false) selects the method for the requested summation without changing the two-argument method's result.

The sorted method currently allocates an array to hold the contributions before sorting and adding them.

_evaluate(a::TaylorN, vals::NTuple)

Return an array containing the evaluated contribution from each polynomial degree of a, starting with its constant term. The entries have not been added together. For the complete result, use the three-argument form with Val(true) to sort before adding, or Val(false) to add without sorting.

_evaluate(a::TaylorN, vals::AbstractVector, ::Val{false})

Evaluate a at vals and add the contributions in increasing polynomial degree, without sorting by their values. This vector method adds each contribution directly, without allocating an array to hold them first.

The third argument distinguishes this complete result from the two-argument _evaluate(a::TaylorN, vals::NTuple), which returns the separate contributions.

_evaluation_order(a, x)

Return the smallest positive order among x and the coefficients of a. Return zero if all those orders are zero. Order-zero series are treated as exact constants, so they do not limit the order of the evaluated result.

_evaluate!(x, vals, dest, ::Val{sorting})

Evaluate each polynomial in x at vals and write its result into dest. Val(false) adds contributions directly without a temporary result array. Val(true) sorts each polynomial's contributions by magnitude before adding them, and may allocate temporary storage for that sorting.

_check_series_evaluation(a, δt, dest, aux)

Check that the evaluation value, destination, and pre-allocated auxiliary aux have compatible orders and, for TaylorN, the same JetSpace. Also check that updating the destination or auxiliary will not overwrite an input or the other output. These checks run before changing dest.

_check_taylorN_evaluation(a, vals, dest, valscache, aux)

Check the number of evaluation values and the orders and JetSpaces of the polynomial, destination, values, and pre-allocated auxiliaries. Also check that the destination and auxiliaries are separate from the inputs and from one another, since evaluation overwrites them. These checks run before changing dest.

_defaultsorting(T, S)

Function to set defaults for sorting, useful for evaluating with Intervals

@isonethread (expr)

Internal macro used to check the number of threads in use, to prevent a data race that modifies coefficient tables when using differentiate or integrate; see https://github.com/JuliaDiff/TaylorSeries.jl/issues/318.

This macro is inspired by the macro @threaded; see https://github.com/trixi-framework/Trixi.jl/blob/main/src/auxiliary/auxiliary.jl; and https://github.com/trixi-framework/Trixi.jl/pull/426/files.