Taylor models functions

#include <audi/taylor_model_functions.hpp>

long long audi::factorial2(int n)

Overload for the double factorial function.

Computes the double factorial of an integer n, denoted as \(n!!\). The double factorial is defined as:

\[\begin{split}n!! = \begin{cases} 1 & \text{if } n = 0 \text{ or } n = -1, \\ n \cdot (n-2) \cdot (n-4) \cdot \ldots \cdot 2 & \text{if } n > 0 \text{ and even}, \\ n \cdot (n-2) \cdot (n-4) \cdot \ldots \cdot 1 & \text{if } n > 0 \text{ and odd}. \end{cases} \end{split}\]

Parameters:

n – the integer argument (must satisfy \( n \geq -1 \))

Throws:

std::invalid_argument – if n < -1

Returns:

the double factorial of n

inline audi::taylor_model audi::pow(const audi::taylor_model &tm, int n)

Overload for the power of a taylor_model with integer exponent.

Computes the power of a taylor_model with an integer exponent. The Taylor polynomial is computed as \((T_f)^n\) using repeated multiplication, while the remainder is propagated according to the interval arithmetic rules of the taylor_model class.

Parameters:

• tm – the taylor_model base

• n – integer exponent

Returns:

a new taylor_model representing \(tm^n\)

template<typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
inline audi::taylor_model audi::pow(const audi::taylor_model &tm, T n)

This is an overloaded member function, provided for convenience. It differs from the above function only in what argument(s) it accepts.

Overload for the sqrt and inverse sqrt of a taylor_model with a real-valued arithmetic exponent. A different, generic real-valued exponent is not implemented.

The remainder for the square root is:

\[R = (-1)^k \sqrt(f_0) \cdot \frac{(2k - 1)!}{(k + 1)!2^{k+1}} \frac{\hat{f}^{k+1}}{f_0^{k+1}} \frac{1}{(1 + \theta \cdot \frac{\hat{f}}{f_0})^{k + \frac{1}{2}}} \]

And for the inverse square root:

\[R = (-1)^{k+1} \frac{1}{\sqrt(f_0)} \cdot \frac{(2k + 1)!!}{(k + 1)!2^{k+1}} \frac{\hat{f}^{k+1}}{f_0^{k+1}} \frac{1}{(1 + \theta \cdot \frac{\hat{f}}{f_0})^{k + \frac{3}{2}}} \]

where \( \forall \vec{x} \in [\vec{a}, \vec{b}] \text{, } B(\hat{f}) + I_f \subset (0, \infty) \)

Template Parameters:

T – arithmetic type of the exponent (e.g. double, float)

Parameters:

• tm – the taylor_model base

• n – real-valued exponent

Returns:

a new taylor_model representing \(tm^n\)

inline audi::taylor_model audi::exp(const audi::taylor_model &tm)

Overload for the exponential of a taylor_model.

Computes the exponential of a taylor_model. For this, the Taylor polynomial is calculated using audi::exp(const gdual<T, M> &d). The remainder term is represented as follows:

\[R = \exp(f_0) \frac{\hat{f}^{k+1}}{(k+1)!} \exp(\theta \cdot \hat{f}) \]

for any \( x \in [\vec{a}, \vec{b}] \text{, } \hat{f} \in B(\hat{f}) \text{, and } 0 < \theta < 1 \).

Parameters:

tm – the taylor_model whose exponential is to be computed

Returns:

a new taylor_model representing \( exp(tm) \)

inline audi::taylor_model audi::log(const audi::taylor_model &tm)

Overload for the log of a taylor_model.

Function calculating the log of a taylor_model.

The remainder for the square root is:

\[R = (-1)^{k+2} \frac{1}{k+1} \cdot \frac{\hat{f}^{k+1}}{f_0^{k+1}} \frac{1}{(1 + \theta \cdot \frac{\hat{f}}{f_0})^{k + 1}} \]

where \( \forall \vec{x} \in [\vec{a}, \vec{b}] \text{, } B(\hat{f}) + I_f \subset (0, \infty) \)

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(log(tm)\)

inline audi::taylor_model audi::sqrt(const audi::taylor_model &tm)

Overload for the sqrt of a taylor_model. This function runs audi::pow(tm, 1/2).

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(sqrt(tm)\)

inline audi::taylor_model audi::sin(const audi::taylor_model &tm)

Overload for the sine of a taylor_model.

Function that calculates the sine of a taylor_model.

The remainder for the sine is:

\[R = \frac{1}{(k+1)!} \bigl(\bar f(\vec{x})\bigr)^{k+1} \cdot J , \]

where

\[\begin{split}J = \begin{cases} -J_0 & \text{if } \operatorname{mod}(k, 4) = 1, 2, \\ J_0 & \text{else}, \end{cases} \end{split}\]

and

\[\begin{split}J_0 = \begin{cases} \cos\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{if } k \text{ is even}, \\ \sin\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{else}. \end{cases} \end{split}\]

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(sin(tm)\)

inline audi::taylor_model audi::cos(const audi::taylor_model &tm)

Overload for the cos of a taylor_model.

Function that calculates the cos of a taylor_model.

The remainder for the cosine is:

\[R = \frac{1}{(k+1)!} \bigl(\bar f(\vec{x})\bigr)^{k+1} \cdot J , \]

where

\[\begin{split}J = \begin{cases} -J_0 & \text{if } \operatorname{mod}(k, 4) = 0, 1, \\ J_0 & \text{else}, \end{cases} \end{split}\]

and

\[\begin{split}J_0 = \begin{cases} \sin\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{if } k \text{ is even}, \\ \cos\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{else}. \end{cases} \end{split}\]

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(cos(tm)\)

inline std::array<audi::taylor_model, 2> audi::sin_and_cos(const audi::taylor_model &tm)

Function for both sine and cosine of a taylor_model.

Function that calculates both the sine and cosine of a taylor_model. Most computations are overlapping, so this saves time. This structure is inspired by sin_and_cos(const gdual<T, M> &d).

Parameters:

tm – the taylor_model base

Returns:

an array of two new taylor_models representing \({sin(tm), cos(tm)}\)

inline audi::taylor_model audi::tan(const audi::taylor_model &tm)

Overload for tangent of taylor_model.

Function that calculates the tangent of a taylor_model. This just uses the usual identity of tan(x) = sin(x) / cos(x)

Parameters:

tm – the taylor_model base

Returns:

an array of two new taylor_models representing \(tan(tm)\)

inline audi::taylor_model audi::sinh(const audi::taylor_model &tm)

Overload for the sinh of a taylor_model.

Function that calculates the sinh of a taylor_model.

The remainder for the hyperbolic sine is:

\[R = \frac{1}{(k+1)!} \bigl(\bar f(\vec{x})\bigr)^{k+1} \cdot J , \]

where

\[\begin{split}J = \begin{cases} \cosh\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{if } k \text{ is even}, \\ \sinh\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{else}. \end{cases} \end{split}\]

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(sinh(tm)\)

inline audi::taylor_model audi::cosh(const audi::taylor_model &tm)

Overload for the cosh of a taylor_model.

Function that calculates the cosh of a taylor_model.

The remainder for the hyperbolic cosine is:

\[R = \frac{1}{(k+1)!} \bigl(\bar f(\vec{x})\bigr)^{k+1} \cdot J , \]

where

\[\begin{split}J = \begin{cases} \sinh\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{if } k \text{ is even}, \\ \cosh\!\bigl(f_0 + \theta \cdot \bar f(\vec{x})\bigr) & \text{else}. \end{cases} \end{split}\]

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(cosh(tm)\)

inline std::array<audi::taylor_model, 2> audi::sinh_and_cosh(const audi::taylor_model &tm)

Function for both sine and cosine of a taylor_model.

Function that calculates both the hyperbolic sine and cosine of a taylor_model. Most computations are overlapping, so this saves time. This structure is inspired by sinh_and_cosh(const gdual<T, M> &d) and analogous to sinh_and_cosh(const audi::taylor_model &tm).

Parameters:

tm – the taylor_model base

Returns:

an array of two new taylor_models representing \({sinh(tm), cosh(tm)}\)

inline audi::taylor_model audi::tanh(const audi::taylor_model &tm)

Overload of the hyperbolic tangent of a taylor_model.

Function that calculates the hyperbolic tangent of a taylor_model. This just uses the usual identity of tanh(x) = sinh(x) / cosh(x)

Parameters:

tm – the taylor_model base

Returns:

an array of two new taylor_models representing \(tanh(tm)\)

int_d audi::asin_derivative(int_d a, unsigned int order)

Function for higher-order derivatives of arcsin.

Computes the \(n\)-th derivative of the arcsine function \(\arcsin(a)\) using a recurrence relation. This function is used internally in the implementation of asin() for taylor_model.

The recurrence is based on:

\[R = \frac{d}{dx} \arcsin(x) = \frac{1}{\sqrt{1 - x^2}} \]

and higher-order derivatives are obtained via:

\[f^{(k+2)}(x) = \frac{(2k+1)x f^{(k+1)}(x) + k^2 f^{(k)}(x)}{1 - x^2}, \quad k \geq 1. \]

Base cases:

\[f^{(1)}(x) = \frac{1}{\sqrt{1-x^2}}, \qquad f^{(2)}(x) = \frac{x}{(1-x^2)^{3/2}}. \]

Parameters:

• a – the interval argument \(a\)

• order – the derivative order (must be \(\geq 1\))

Throws:

std::out_of_range – if order < 1

Returns:

the order -th derivative of \(\arcsin(a)\)

inline audi::taylor_model audi::asin(const audi::taylor_model &tm)

Overload for the asin of a taylor_model.

Function that calculates the asin of a taylor_model.

Here, some custom arithmetic is used from Makino (1998):

\( \arcsin(f) = \arcsin(f_0) + \arcsin(g) \)

where \( g = f \cdot \sqrt{1 - f_0^2} - f_0 \cdot \sqrt{1 - f^2} \)

and where \( \forall \vec{x} \in [\vec{a}, \vec{b}] \text{, } B(f) + I_f \subset (-1, 1) \) applies.

The remainder for the arcsine is:

\[R = \frac{1}{(k+1)!} g^{k+1} \cdot \arcsin^{(k+1)}(\theta \cdot g) \]

Here, \( \arcsin^{(k+1)}(\theta \cdot g) \) is the higher order derivative for an arcsine: see Makino (1998).

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(asin(tm)\)

inline audi::taylor_model audi::acos(const audi::taylor_model &tm)

Overload for the acos of a taylor_model.

Function that calculates the acos of a taylor_model.

Here, we use the following identity utilizing the arcsine. So see asin(const audi::taylor_model &tm) for the implementation.

\( \arccos(f) = \frac{\pi}{2} - \arcsin(f) \)

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(acos(tm)\)

inline audi::taylor_model audi::atan(const audi::taylor_model &tm)

Overload for the atan of a taylor_model.

Function that calculates the atan of a taylor_model.

Here, some custom arithmetic is used from Makino (1998):

\( \arctan(f) = \arctan(f_0) + \arctan(g) \)

where \( g = \frac{\hat{f}}{1 + f_0 \cdot f} \)

The remainder for the arcsine is:

\[R = \frac{1}{(k+1)!} g^{k+1} \cdot \arccos^{(k+1)}(\arctan(\theta \cdot g)) \cdot \sin((k + 1) \cdot (\arctan(\theta \cdot g) + \frac{\pi}{2})) \]

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(atan(tm)\)

inline audi::taylor_model audi::asinh(const audi::taylor_model &tm)

Overload for the asinh of a taylor_model.

Function that calculates the asinh of a taylor_model.

Here, the following identity is used:

\( arcsinh(f) = \log(f + \sqrt(f^2 + 1)) \)

For the remainder, see log(const audi::taylor_model &tm), sqrt(const audi::taylor_model &tm)

, and pow(const

audi::taylor_model &tm).

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(arcsinh(tm)\)

inline audi::taylor_model audi::acosh(const audi::taylor_model &tm)

Overload for the acosh of a taylor_model.

Function that calculates the acosh of a taylor_model.

Here, the following identity is used:

\( arccosh(f) = \log(f + \sqrt(f^2 - 1)) \)

where \( B(f) + I_f \subset (1, \infty) \) to prevent complex values from the sqrt.

For the remainder, see log(const audi::taylor_model &tm), sqrt(const audi::taylor_model &tm)

, and pow(const

audi::taylor_model &tm).

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(arccosh(tm)\)

inline audi::taylor_model audi::atanh(const audi::taylor_model &tm)

Overload for the atanh of a taylor_model.

Function that calculates the atanh of a taylor_model.

Here, the following identity is used:

\( arctanh(f) = \frac{1}{2} \cdot \log(\frac{1 + f}{1 - f}) \)

For the remainder, see log(const audi::taylor_model &tm).

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(arctanh(tm)\)

inline audi::taylor_model audi::abs(const audi::taylor_model &tm)

Overload for the abs of a taylor_model.

UNVERIFIED. Function that calculates the absolute value of a taylor_model.

For the remainder, it is assumed that one just takes the absolute value of the remainder bound.

Parameters:

tm – the taylor_model base

Returns:

a new taylor_model representing \(abs(tm)\)