block = [draw, fill=white, rectangle,
minimum height=3em, minimum width=6em]
sum = [draw, fill=white, circle, node distance=1cm]
input = [coordinate]
output = [coordinate]
pinstyle = [pin edge=to-,thin,black]

Exact Real Arithmetic with

Perturbation Analysis and Proof of Correctness

Sarmen Keshishzadeh
Jan Friso Groote
Department of Mathematics and Computer Science

Eindhoven University of Technology

Den Dolech 2, 5612 AZ Eindhoven, The Netherlands

In this article, we consider a simple representation for real numbers and propose top-down procedures to approximate various algebraic and transcendental operations with arbitrary precision. Detailed algorithms and proofs are provided to guarantee the correctness of the approximations. Moreover, we develop and apply a perturbation analysis method to show that our approximation procedures only recompute expressions when unavoidable.

In the last decade, various theories have been developed and implemented to realize real computations with arbitrary precision. Proof of correctness for existing approaches typically consider basic algebraic operations, whereas detailed arguments about transcendental operations are not available. Another important observation is that in each approach some expressions might require iterative computations to guarantee the desired precision. However, no formal reasoning is provided to prove that such iterative calculations are essential in the approximation procedures. In our approximations of real functions, we explicitly relate the precision of the inputs to the guaranteed precision of the output, provide full proofs and a precise analysis of the necessity of iterations.

§ INTRODUCTION

Various scientific disciplines use computations involving real numbers to model and reason about different phenomena in the world. Real numbers are typically approximated by floating point numbers in scientific calculations. Round-off errors are inevitable in such approximations and they might build up into catastrophic errors in some cases. Exact real arithmetic approaches address this issue by devising computation procedures that given an expression and a precision requested by the user produce an output that is guaranteed to meet the precision requirement.

Several approaches [16, 13] to exact real arithmetic are based on iterative bottom-up calculations. Given an expression and a desired precision, bottom-up approaches typically start with calculating the inputs with an arbitrary precision higher than the requested precision. Then, the sub-expressions are evaluated in a bottom-up way. After evaluating the sub-expressions of each level, the guaranteed precision is passed to the higher level. These calculations proceed until the main expression is calculated and its guaranteed precision is determined. If the precision obtained for the expression is not adequate, the computation restarts with increased precisions for the inputs.

In contrast, various authors [3, 15] have proposed top-down approaches to exact real arithmetic. In top-down approaches, the required precision of each sub-expression is determined based on the precision required for its immediate parent expression. For certain types of expressions, the required precision of the sub-expressions can be calculated immediately. However, some expressions may require to first obtain additional information about the magnitude of the values of their sub-expressions before determining their required precision. Thus, in general, it might be necessary to recompute certain expressions.

The main benefit of top-down approaches is that they exploit the structure of a given expression to estimate the required precision of the sub-expressions. In this context, one would ideally like to have top-down approximations for algebraic and transcendental functions such that

* the approximations are proven to be correct and

* iterative calculations are avoided unless they are proven to be necessary.

In several studies [3, 17], proofs of correctness for algebraic operations are available. However, the arguments about transcendental functions provide little insight about the correctness of the approximations and the effect of these operations on precision. Taylor expansions are the most prominent way to approximate transcendental functions. Calculations with Taylor expansions are typically restricted to a base interval; range reduction identities are used to extend the computations to the complete domain of a function. Proofs of correctness for transcendental functions are limited to the base interval [18, 3, 17], whereas little attention is given to the general case where the computations introduced by range reduction identities influence the output precision.

The second desired property for a top-down approach is related to the iterative nature of the computations. As discussed above, bottom-up and top-down approaches rely on iterative computation schemes. However, no formal reasoning is provided to prove that such iterative calculations are essential in the approximation procedures.

In this article, we consider a simple representation for real numbers and propose a top-down approach to approximate various algebraic and transcendental functions with arbitrary precision. For each operation, we describe an approximation procedure and relate the precision of the inputs to the guaranteed precision for the output. To guarantee the correctness of our top-down approach, we provide detailed proofs of correctness for the proposed approximations.

To identify computational problems that require iterative calculations in our top-down approach, we have developed a perturbation analysis method. Our analysis describes the influence of errors in the inputs of a computational problem on the output precision. We apply perturbation analysis to show that our approximation procedures only recompute expressions when this is unavoidable.

We discuss different approaches to defining computability of functions in Section <ref>. In Section <ref> we introduce a representation for real numbers and specify the syntax of the expressions that we consider in our computations. To analyze computational problems in this setting, a perturbation analysis method is introduced in Section <ref>. In Section <ref> we discuss our approximations of algebraic operations. We approximate transcendental functions using Riemann sums and Taylor expansions in Section <ref> and <ref>, respectively. Section <ref> contains discussions about related work. In Section <ref> we draw some conclusions and suggest directions for future research.

§ COMPUTABLE REAL FUNCTIONS

Real arithmetic is concerned with performing computations on  real numbers. In order to do calculations with real numbers, it is necessary to define what it means for an  operation to be computable. In this section we briefly discuss different approaches to defining computability.

Since real numbers are infinite objects, one can use infinite streams from a finite alphabet $\Sigma$ to represent them. This gives rise to a definition of computability called Type-$2$ Theory of Effectivity (TTE, [21]). In TTE computable operations are defined in terms of functions $f:\Sigma^{\omega}\rightarrow \Sigma^{\omega}$ that receive infinite words as input and produce infinite words as output. An essential property of a Type-$2$ computable function is the finiteness property [21]. This property indicates that for a computable function $f$, any finite prefix of the output $f(x)$ is computable by some finite portion of the input $x$.

An alternative definition of computability has been introduced by the Russian school of constructive analysis [14, 12]. In this definition, computable operations are defined in terms of functions $f:\Sigma^*\rightarrow \Sigma^*$ that receive finite words as input and produce finite words as output. This approach is sometimes referred to as Type-$1$ computability. A function $f$ is Type-$1$ computable if there is a Turing machine that transforms any finite input string $x\in \Sigma^*$ to the finite output strings $f(x)\in \Sigma^*$. Type-$1$ machines provide a natural way to define computability on, for instance, rational numbers and finite graphs.

In both Type-$1$ and Type-$2$ approaches, $\Sigma$ depends on the concrete representation that we use for input/output objects. For instance, one can use the binary signed-digit representation to represent real numbers in the inputs/outputs of computations.

The relation between Type-$1$ and Type-$2$ computable functions has been investigated in various studies. It is known that restricting the domain of a Type-$2$ computable function to finite streams results in a Type-$1$ computable function [4, 19]. However, not every Type-$1$ function can be obtained by restricting some Type-$2$ computable function [19, 11]. To illustrate this, we consider a function $f_1: \{0,1\}^*\rightarrow \{0,1\}^*$ defined as follows:

\begin{align*}
\begin{cases}
0& \text{if } s=0^k\\
1& \text{if } s=0^k1s'
\end{cases}
\end{align*}

where $k\in \nat$ is the length of the longest prefix of zeros in $s$ and $s'\in \Sigma^*$ is a finite suffix of $s$ in the second case. The function $f_1$ performs computations on finite strings and one can construct a Type-$1$ Turing machine to compute this function. By extending the domain of $f_1$ to infinite strings, we obtain $f_2: \Sigma^{\omega} \rightarrow \Sigma^{\omega}$ such that:

\begin{align*}
\begin{cases}
0& \text{if } s=0^\omega\\
1& \text{if } s=0^k1s'
\end{cases}
\end{align*}

The function $f_2$ is not computable with a Type-$2$ Turing machine; it is not possible to write $0$ in the output after reading a finite prefix from the input. The interested reader can refer to [4] for more details about the relation between Type-$1$ and Type-$2$ computable functions.

In addition to Type-$1$ and Type-$2$, one can also consider a third approach to defining computability based on certain finite structures that provide precise descriptions for specific classes of real numbers. For instance, Lagrange's theorem on continued fractions indicates that the real numbers whose continued fraction is periodic are the quadratic irrationals. Based on this observation, one can define computability in terms of functions $f$ that given a finite and precise representation of $x$ produce a finite and precise representation of $f(x)$.

The exact real arithmetic approach that we introduce in this article is based on Type-$2$ computability. In Section <ref> we discuss a representation for real numbers in terms of rational numbers that are coupled with a notion of precision. Our approximations for arithmetic operations rely on the finiteness property of computable functions. Thus, for a given computational problem, a desired precision for the output is obtained based on sufficiently good approximations of the inputs.

§ REAL NUMBERS: REPRESENTATION & OPERATIONS

In this section we first discuss our representation of real numbers and then describe the syntax of the expressions that we aim to calculate in our setting.

Since real numbers are infinite objects, a finite representation of an arbitrary real number $x$ can only approximate $x$ with a certain precision. In scientific measurements and calculations, the amount of error that we commit in approximations is measured by an absolute or relative error. In practice, an absolute error is of little use. Since numbers tend to have very different orders of magnitude, it is the relative error that shows the significance of the lost digits in measurements or calculations. Hence, in our setting we use a representation based on the relative error.

A real number $x$ is represented by a tuple $(m,n,p)$ such that:

\begin{align*}
|x-\frac{m}{n}| < |\frac{m}{n}|\frac{1}{2^p}
\end{align*}

where $m,n\in \nonzerointeger, p\in \nat$.

The representation $(m,n,p)$ for $x$ means that $\frac{m}{n}$ approximates $x$ and the relative error of this approximation does not exceed $\frac{1}{2^p}$.

In this article, we focus on calculating expressions that can be described with the following grammar:

\begin{equation}
\begin{aligned}
~E\cdot E~|
&	~\ln(E)~|
\end{aligned}
\end{equation}

where $c$ represents a rational constant.

The following identities show that other interesting operations can be described in terms of the operations of this grammar:

\begin{align*}
\tan(x)&=\frac{\sin(x)}{\cos(x)}\\
\cot(x)&=\frac{1}{\tan(x)}\\
\arcsin(x)&=\arctan(\frac{x}{\sqrt{1-x^2}})\\
\arccos(x)&=\arctan(\frac{\sqrt{1-x^2}}{x})\\
\arccot(x)&=\arccos(\frac{x}{\sqrt{1+x^2}})
\end{align*}

§ SENSITIVITY OF OPERATIONS TO PERTURBATIONS IN THE ARGUMENTS

Our goal is to develop a top-down exact real arithmetic approach based on the representation of Definition <ref>. Thus, for a given computational problem, it is essential to estimate the required precision of the inputs based on the desired precision in the output. Moreover, we would like to investigate to what extend the operations of  grammar (<ref>) can be calculated in a top-down manner without iterations.

To analyze the operations of grammar (<ref>), we introduce a pertubation analysis method for measuring the sensitivity of the operations to perturbations in their arguments. We consider two general cases in our analysis. First, we consider a function $f(x)$ in one variable and show how errors in the input influence the output (Section <ref>). Then, we consider a function $f(x,y)$ with two arguments and investigate the effect of errors in the inputs on the output (Section <ref>).

§.§ Perturbation Analysis for Unary Functions

Let $f(x)$ be a differentiable function that we want to calculate in point $x=a$. Suppose that $\dx$ is a perturbation in the argument $a$. The relative error in the calculation of $f(a)$ caused by this perturbation is:

\begin{align*}
\end{align*}

We want to relate this relative error to the relative error of the argument, namely $|\frac{\dx}{a}|$. To this end, we use the following approximation of the function $f$ in point $x=a+\dx$:

\begin{align*}
f(a+\dx)\approx f(a)+f'(a)\dx
\end{align*}

We can approximate the relative error of $f$ as follows:

\begin{align}
|\frac{f(a+\dx)-f(a)}{f(a)}|\approx |\frac{f'(a)\dx}{f(a)}|=|\frac{af'(a)}{f(a)}||\frac{\dx}{a}|\label{eqn:condition_single}
\end{align}

From equality (<ref>) one can see that the quantity $|\frac{af'(a)}{f(a)}|$ determines the effect of the relative error $|\frac{\dx}{a}|$ on the output. In numerical analysis and linear algebra the quantity $|\frac{xf'(x)}{f(x)}|$ is usually referred to as the condition number of $f(x)$ [5, 6].

§.§ Perturbation Analysis for Binary Functions

Let $f(x,y)$ be a differentiable function that we want to calculate in point $(x,y)=(a,b)$. Suppose $\dx$ and $\dy$ are perturbations in the arguments $a$ and $b$, respectively. The relative error of $f(a,b)$ caused by these perturbations can be calculated as follows:

\begin{align}
\end{align}

To find an upper-bound for (<ref>), we use the following first-order approximation of the function $f$ in point $(x,y)=(a+\dx, b+\dy)$:

\begin{align*}
f(a+\dx,b+\dy)\approx f(a,b)+f_x(a,b)\dx+f_y(a,b)\dy
\end{align*}

where $f_x(a,b)$ and $f_y(a,b)$ are the partial derivatives of $f$ with respect to $x$ and $y$ in $(a,b)$, respectively. We relate the relative error in the calculation of $f$ to the relative errors $|\frac{\dx}{a}|$ and $|\frac{\dy}{b}|$ as follows:

\begin{align}
|\frac{f(a+\dx,b+\dy)-f(a,b)}{f(a,b)}|&\approx |\frac{f_x(a,b)\dx +f_y(a,b)\dy}{f(a,b)}|\leq \nonumber\\
& |\frac{af_x(a,b)}{f(a,b)}||\frac{\dx}{a}|+|\frac{bf_y(a,b)}{f(a,b)}||\frac{\dy}{b}|\leq \nonumber\\
& (|\frac{af_x(a,b)}{f(a,b)}|+|\frac{bf_y(a,b)}{f(a,b)}|)\max\{|\frac{\dx}{a}|,|\frac{\dy}{b}|\}\label{eqn:condition_multiple}
\end{align}

The quantity $|\frac{af_x(a,b)}{f(a,b)}|+|\frac{bf_y(a,b)}{f(a,b)}|$ determines the upper-bound calculated in inequality (<ref>) and we use this quantity to measure the effect of erroneous arguments on the output. It should be noted that in inequality (<ref>) we have considered $|\frac{\dx}{a}|$ and $|\frac{\dy}{b}|$ as independent factors that can influence the relative error of $f$. This way of reasoning about the sensitivity of $f(x,y)$ is related to componentwise analysis of perturbation in numerical analysis and linear algebra [7], which we use in this article.

Another possibility is to relate the relative error of (<ref>) to the quantity:

\begin{align*}
\frac{
\begin{Vmatrix}
\begin{bmatrix}
\dx \\ \dy
\end{bmatrix}
\end{Vmatrix}  }
\begin{bmatrix}
a \\ b
\end{bmatrix}
\end{Vmatrix}}
\end{align*}

This type of analysis is usually referred to as normwise analysis of perturbation [7].

In the following sections, we provide a top-down approach for approximating various algebraic and transcendental functions. Perturbation analysis will be used to show that  in our approximations we only recompute expressions when essential.

§ APPROXIMATING ALGEBRAIC OPERATIONS

In this section we calculate the algebraic operations of grammar (<ref>) using a top-down approach. We formulate and prove theorems that allow us to calculate expressions involving unary negation (Section <ref>), multiplication (Section <ref>), inverse (Section <ref>), addition (Section <ref>), and square root (Section <ref>). Based on the theorems, we provide different implementations of Compute$(expr,p)$ to calculate algebraic operations. These implementations receive an algebraic expression $expr$ and a desired precision $p$ and produce an output with the desired precision. In each case, we also apply the perturbation analysis of Section <ref> and show that we avoid unnecessary iterations in our approximations.

§.§ Unary Negation

Let $x$ be a real number represented by $(m,n,p)$. Then $-x$ can be represented by $(-m,n,p)$.

Since $x$ is represented by $(m,n,p)$ we can write:

\begin{align*}
\frac{m}{n}-|\frac{m}{n}|\frac{1}{2^p}&<x<\frac{m}{n}+|\frac{m}{n}|\frac{1}{2^p}\\
\end{align*}

Thus, we can represent $-x$ by $(-m,n,p)$.

Algorithm <ref> applies Theorem <ref> and approximates $-x$ based on a representation $(m,n,p)$ of $x$. To confirm that $-x$ can be approximated with arbitrary precision in one pass, we calculate $|\frac{xf'(x)}{f(x)}|$ for the function $f(x)=-x$:

\begin{align*}
\end{align*}

The quantity $|\frac{xf'(x)}{f(x)}|$ is small and independent of the argument $x$ and hence the amount of precision that we lose in unary negation (which is $0$ in the case of Theorem <ref>) can be calculated independently of $x$.

Unary Negation
[1]
$\expr$ has the shape $-x$

$\frac{m}{n} \gets \Call{Compute}{x,p}$

§.§ Multiplication $x\cdot y$

Let $x$ and $y$ be two real numbers represented by $(m,n,p)$ and $(m',n',p)$, respectively. Then $x\cdot y$ can be represented by $(mm',nn',p-2)$.

From Definition <ref> we can write:

\begin{align}
\frac{m}{n}-|\frac{m}{n}|\frac{1}{2^{p}}<x<\frac{m}{n}+|\frac{m}{n}|\frac{1}{2^{p}} \label{eqn:mult_x}\\
\frac{m'}{n'}-|\frac{m'}{n'}|\frac{1}{2^{p}}<y<\frac{m'}{n'}+|\frac{m'}{n'}|\frac{1}{2^{p}} \label{eqn:mult_y}
\end{align}

We consider three cases:

* Suppose $\frac{mm'}{nn'}>0$. We can multiply inequalities (<ref>) and (<ref>) as follows:

\begin{align*}
\frac{mm'}{nn'}(1-\frac{1}{2^{p}})^2<x\cdot y <\frac{mm'}{nn'}(1+\frac{1}{2^{p}})^2\
\end{align*}

If $x\cdot y$ can be represented by $(mm',nn',p-2)$ then it must be the case that:

\begin{align*}
\frac{mm'}{nn'}(1-\frac{1}{2^{p-2}})<x\cdot y<\frac{mm'}{nn'}(1+\frac{1}{2^{p-2}})
\end{align*}

To see that this is valid, we need to show that:

\begin{align*}
\frac{mm'}{nn'}(1+\frac{1}{2^{p}})^2\leq \frac{mm'}{nn'}(1+\frac{1}{2^{p-2}})\\
\frac{mm'}{nn'}(1-\frac{1}{2^{p-2}}) \leq \frac{mm'}{nn'}(1-\frac{1}{2^{p}})^2
\end{align*}

But $\frac{mm'}{nn'}>0$ and from Proposition <ref> (see <ref>) we know that both inequalities hold.

* Suppose $\frac{m}{n}>0, \frac{m'}{n'}<0$. We can rewrite inequalities (<ref>) and (<ref>) as follows:

\begin{align}
\frac{m}{n}(1-\frac{1}{2^{p}})<x<\frac{m}{n}(1+\frac{1}{2^{p}})\label{eqn:mult_x_pos} \\
\frac{m'}{n'}(1+\frac{1}{2^{p}})<y<\frac{m'}{n'}(1-\frac{1}{2^{p}}) \label{eqn:mult_y_neg}
\end{align}

Multiplying inequalities (<ref>) and (<ref>) we get:

\begin{align*}
\frac{mm'}{nn'}(1+\frac{1}{2^{p}})^2<x\cdot y<\frac{mm'}{nn'}(1-\frac{1}{2^{p}})^2
\end{align*}

If $x\cdot y$ is representable by $(mm',nn',p-2)$ then it must be the case that:

\begin{align*}
\frac{mm'}{nn'}(1+\frac{1}{2^{p-2}})<x\cdot y<\frac{mm'}{nn'}(1-\frac{1}{2^{p-2}})
\end{align*}

To show that this is valid, it suffices to prove that:

\begin{align*}
\frac{mm'}{nn'}(1-\frac{1}{2^{p}})^2 \leq \frac{mm'}{nn'}(1-\frac{1}{2^{p-2}}) \\
\frac{mm'}{nn'}(1+\frac{1}{2^{p-2}})\leq \frac{mm'}{nn'}(1+\frac{1}{2^{p}})^2
\end{align*}

But $\frac{mm'}{nn'}<0$ and from Proposition <ref> (see <ref>) we know that both inequalities hold.

* Suppose $\frac{m}{n}<0, \frac{m'}{n'}>0$. This case can be proved similarly to the second case.

[1]
$\expr$ has the shape $x\cdot y$

$\frac{m}{n} \gets \Call{Compute}{x,p+2}$
$\frac{m'}{n'} \gets \Call{Compute}{y,p+2}$

Algorithm <ref> depicts an approximation of $x\cdot y$ based on Theorem <ref>. Given approximations of $x$ and $y$, this algorithm approximates $x\cdot y$ in one pass; the loss of precision is predictable without additional information about $x$ and $y$. To confirm this claim, we apply the perturbation analysis of Section <ref> on the function $f(x,y)=x\cdot y$:

\begin{align*}
\end{align*}

The sensitivity measure is small and independent of the arguments. Thus, in a top-down approach, $x\cdot y$ can be approximated in one pass without iterative computations of $x$ and $y$.

§.§ Inverse

Let $x$ be a real number represented by $(m,n,p)$. Then $\frac{1}{x}$ can be represented by $(n,m,p-1)$.

We consider two cases:

* Suppose $\frac{m}{n}>0$. Since $x$ is represented by $(m,n,p)$ we can write:

\begin{align*}
\frac{m}{n}(1-\frac{1}{2^{p}})<x<\frac{m}{n}(1+\frac{1}{2^{p}})\\
\frac{n}{m}(\frac{2^{p}}{2^{p}+1})<\frac{1}{x}<\frac{n}{m}(\frac{2^{p}}{2^{p}-1})
\end{align*}

If $\frac{1}{x}$ is representable by $(n,m,p-1)$, then it must be the case that:

\begin{align*}
\frac{n}{m}(1-\frac{1}{2^{p-1}})<\frac{1}{x}<\frac{n}{m}(1+\frac{1}{2^{p-1}})
\end{align*}

To see that this is valid, it suffices to show:

\begin{align*}
\frac{n}{m}(\frac{2^{p}}{2^{p}-1})\leq \frac{n}{m}(1+\frac{1}{2^{p-1}})\\
\frac{n}{m}(1-\frac{1}{2^{p-1}})\leq \frac{n}{m}(\frac{2^{p}}{2^{p}+1})
\end{align*}

But $\frac{n}{m}>0$ and hence both inequalities follow from Proposition <ref> (see <ref>).

* Suppose $\frac{m}{n}<0$. Since $x$ is representable by $(m,n,p)$ we have:

\begin{align*}
\frac{m}{n}(1+\frac{1}{2^{p}})<x<\frac{m}{n}(1-\frac{1}{2^{p}})\\
\frac{n}{m}(\frac{2^{p}}{2^{p}-1})<\frac{1}{x}<\frac{n}{m}(\frac{2^{p}}{2^{p}+1})
\end{align*}

If $\frac{1}{x}$ is representable by $(n,m,p-1)$ then it must be the case that:

\begin{align*}
\frac{n}{m}(1+\frac{1}{2^{p-1}})<\frac{1}{x}<\frac{n}{m}(1-\frac{1}{2^{p-1}})
\end{align*}

To see that this is valid, we need to show:

\begin{align*}
\frac{n}{m}(\frac{2^{p}}{2^{p}+1})\leq \frac{n}{m}(1-\frac{1}{2^{p-1}}) \\
\frac{n}{m}(1+\frac{1}{2^{p-1}})\leq \frac{n}{m}(\frac{2^{p}}{2^{p}-1})
\end{align*}

But $\frac{n}{m}<0$ and hence the inequalities follow from Proposition <ref> (see <ref>).

[1]
$expr$ has the shape $\frac{1}{x}$

$\frac{m}{n} \gets \Call{Compute}{x,p+1}$

Algorithm <ref> approximates $\frac{1}{x}$ with precision $p$ based on Theorem <ref>. Given an approximation of $x$ with precision $p$, the algorithm allows us to approximate $\frac{1}{x}$ in one pass. We use perturbation analysis and calculate the quantity $|\frac{xf'(x)}{f(x)}|$ for $f(x)=\frac{1}{x}$ to show that loss of precision in the inverse can be estimated independently of the argument:

\begin{align*}
\end{align*}

The quantity $|\frac{xf'(x)}{f(x)}|$ is a constant and hence iterative computations can be avoided when calculating the inverse.

§.§ Addition

Let $x$ and $y$ be two real numbers represented by $(m,n,p)$ and $(m',n',p)$, respectively. The value of $x+y$ can be approximated as follows:

If $\frac{mm'}{nn'}>0$, then $x+y$ can be represented by $(mn'+m'n,nn',p)$.

If $\frac{mm'}{nn'}<0$ and $i\in \nat^+$ is the smallest natural number such that $i\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}}{1-\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}})$, then $x+y$ can be represented by $(mn'+m'n,nn',p-i)$.

* For numbers $x$ and $y$ we can write:

\begin{align}
\frac{m}{n}-|\frac{m}{n}|\frac{1}{2^{p}}&<x<\frac{m}{n}+|\frac{m}{n}|\frac{1}{2^{p}} \label{eqn:add_x}\\
\frac{m'}{n'}-|\frac{m'}{n'}|\frac{1}{2^{p}}&<y<\frac{m'}{n'}+|\frac{m'}{n'}|\frac{1}{2^{p}} \label{eqn:add_y}
\end{align}

From inequalities (<ref>) and (<ref>) we can write:

\begin{align}
\end{align}

If $x+y$ is representable by $(mn'+m'n,nn',p)$, then it must be the case that:

\begin{align*}
\end{align*}

To show that this is valid, we need to prove that:

\begin{align}
(\frac{m}{n}+\frac{m'}{n'})+\frac{1}{2^{p}}(|\frac{m}{n}|+|\frac{m'}{n'}|)\leq (\frac{m}{n}+\frac{m'}{n'})+\frac{1}{2^{p}}|\frac{m}{n}+\frac{m'}{n'}| \label{eqn:add_samesign_pos} \\
(\frac{m}{n}+\frac{m'}{n'})-\frac{1}{2^{p}}|\frac{m}{n}+\frac{m'}{n'}| \leq (\frac{m}{n}+\frac{m'}{n'})-\frac{1}{2^{p}}(|\frac{m}{n}|+|\frac{m'}{n'}|) \label{eqn:add_samesign_neg}
\end{align}

The rational numbers $\frac{m}{n}$ and $\frac{m'}{n'}$ have the same sign. Therefore, $|\frac{m}{n}+\frac{m'}{n'}|=|\frac{m}{n}|+|\frac{m'}{n'}|$ holds and inequalities (<ref>) and (<ref>) are valid.

* If $x+y$ is representable by $(mn'+m'n,nn',p-i)$, then it must be the case that:

\begin{align*}
\end{align*}

To show that this holds, we need to prove the following (see inequality (<ref>)):

\begin{align*}
(\frac{m}{n}+\frac{m'}{n'})+\frac{1}{2^{p}}(|\frac{m}{n}|+|\frac{m'}{n'}|)\leq (\frac{m}{n}+\frac{m'}{n'})+\frac{1}{2^{p-i}}|\frac{m}{n}+\frac{m'}{n'}|\\
(\frac{m}{n}+\frac{m'}{n'})-\frac{1}{2^{p-i}}|\frac{m}{n}+\frac{m'}{n'}|\leq (\frac{m}{n}+\frac{m'}{n'})-\frac{1}{2^{p}}(|\frac{m}{n}|+|\frac{m'}{n'}|)
\end{align*}

For both inequalities, it boils down to proving the following:

\begin{align}
\frac{1}{2^i}(|\frac{m}{n}|+|\frac{m'}{n'}|)\leq |\frac{m}{n}+\frac{m'}{n'}|\label{eqn:add_diff_sgn}
\end{align}

To prove inequality (<ref>), we consider the following two cases:

* Suppose $\frac{m}{n}>0$ and $\frac{m'}{n'}<0$. We can rewrite inequality (<ref>) as follows:

\begin{align*}
|\frac{m}{n}+\frac{m'}{n'}| \geq \frac{1}{2^i}(\frac{m}{n}-\frac{m'}{n'}) \Leftrightarrow & (\frac{m}{n}+\frac{m'}{n'})\geq \frac{1}{2^i}(\frac{m}{n}-\frac{m'}{n'}) \vee \\
&(\frac{m}{n}+\frac{m'}{n'})\leq \frac{1}{2^i}(\frac{m'}{n'}-\frac{m}{n})
\end{align*}

In other words, we should show $\frac{\frac{m}{n}}{-\frac{m'}{n'}}\geq \frac{2^i+1}{2^i-1}$ or $\frac{-\frac{m'}{n'}}{\frac{m}{n}}\geq \frac{2^i+1}{2^i-1}$. Depending on the values of $|\frac{m}{n}|$ and $|\frac{m'}{n'}|$, both cases follow from $\frac{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}\geq \frac{2^i+1}{2^i-1}$. Equivalently, we should have $i\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}}{1-\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}})$.

* Suppose $\frac{m}{n}<0$ and $\frac{m'}{n'}>0$. We can rewrite inequality (<ref>) as follows:

\begin{align*}
|\frac{m}{n}+\frac{m'}{n'}|\geq \frac{1}{2^i}(\frac{m'}{n'}-\frac{m}{n}) \Leftrightarrow & (\frac{m}{n}+\frac{m'}{n'})\geq \frac{1}{2^i}(\frac{m'}{n'}-\frac{m}{n}) \vee \\
& (\frac{m}{n}+\frac{m'}{n'}) \leq \frac{1}{2^i}(\frac{m}{n}-\frac{m'}{n'})
\end{align*}

In other word, we need to prove $\frac{\frac{m'}{n'}}{-\frac{m}{n}}\geq \frac{2^i+1}{2^i-1}$ or $\frac{-\frac{m}{n}}{\frac{m'}{n'}}\geq \frac{2^i+1}{2^i-1}$. Depending on the values of $|\frac{m}{n}|$ and $|\frac{m'}{n'}|$, both cases follow from $\frac{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}\geq \frac{2^i+1}{2^i-1}$. Equivalently, we should have $i\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}}{1-\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}})$.

[1]
$\expr$ has the shape $x+y$

$dp \gets p$

$\frac{m}{n} \gets \Call{Compute}{x,dp}$
$\frac{m'}{n'} \gets \Call{Compute}{y,dp}$
$\frac{mm'}{nn'}>0$ Theorem <ref>.<ref>
Theorem <ref>.<ref>
$i\gets \lceil \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}}{1-\frac{\min(|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}}) \rceil$
$dp-i\geq p$
$dp\gets dp+1$

Algorithm <ref> applies Theorem <ref> to approximate $x+y$ with precision $p$. If approximations $\frac{m}{n}$ and $\frac{m'}{n'}$ have the same sign, we do not lose precision by calculating $x+y$. On the other hand, if $\frac{m}{n}$ and $\frac{m'}{n'}$ have different signs, the amount of precision that is lost depends on the magnitude of $\frac{\min (|\frac{m}{n}|,|\frac{m'}{n'}|)}{\max(|\frac{m}{n}|,|\frac{m'}{n'}|)}$. This indicates that if $\frac{mm'}{nn'}<0$ and $|\frac{m}{n}|\approx |\frac{m'}{n'}|$, a significant amount of precision can be lost in $x+y$. Thus, if the guaranteed precision for $x+y$ (i.e., $p-i$) is not sufficient, $x$ and $y$ must be recomputed with higher precisions (see Line <ref>,<ref> in Algorithm <ref>).

To confirm this observation, we apply the perturbation analysis of Section <ref> on $f(x,y)=x+y$:

\begin{align*}
\end{align*}

The quantity $|\frac{x}{x+y}|+|\frac{y}{x+y}|$ is $1$ when $xy>0$. Thus, we can estimate the loss of precision in $x+y$ independently of the arguments when $xy>0$.

However, the quantity $|\frac{x}{x+y}|+|\frac{y}{x+y}|$ can be arbitrarily large when $xy<0$ and $|x+y|\approx 0$. This confirms that a significant amount of precision can be lost in $x+y$.

Perturbation analysis shows that in general, we cannot estimate the amount of precision that is lost in $x+y$ independently of the arguments. Hence, if approximation $\frac{m}{n}$ and $\frac{m'}{n'}$ have different signs, recomputing $x$ and $y$  might be essential to obtain the desired precision for $x+y$. Loss of precision in $x+y$ is sometimes referred to as loss of significance [9] or catastrophic cancellation [2].

It should be noted that Theorem <ref> does not imply that $x+y$ is always fundamentally problematic when $xy<0$ and $|x+y|\approx 0$. In certain cases, the calculation can be adjusted in such a way that loss of significance can be avoided and the expression can be calculated in one pass.

Suppose we want to calculate $\sqrt{x+1}-\sqrt{x}$ for a relatively large $x$. Since $\sqrt{x+1}\approx\sqrt{x}$, we will lose a significant amount of precision if we directly calculate $\sqrt{x+1}-\sqrt{x}$. However, we can change the calculation algorithm by rewriting the expression as follows:

\begin{align*}
\sqrt{x+1}-\sqrt{x}=(\sqrt{x+1}-\sqrt{x})\times \frac{\sqrt{x+1}+\sqrt{x}}{\sqrt{x+1}+\sqrt{x}}=\frac{1}{\sqrt{x+1}+\sqrt{x}}
\end{align*}

In the new expression, all the operations can be approximated with a desired precision in one pass (see Section <ref> on calculating square root). Hence, we can approximate the new expression without recomputing the sub-expressions with higher precisions. Applying the perturbation analysis of Section <ref> also shows that $f(x)=\sqrt{x+1}-\sqrt{x}$ is not fundamentally problematic:

\begin{align*}
\end{align*}

Using perturbation analysis, we can identify instances of $x+y$ where adjustments in the algorithm can avoid loss of significance. However, to our knowledge a general scheme for making such adjustments does not exist.

§.§ Square Root

In this section, we calculate $\sqrt{x}$ by approximating the root of $f(y)=y^2-x$ using the Newton-Raphson method [9].
The Newton-Raphson method starts with an initial approximation $y_0$ for $\sqrt{x}$ and iteratively generates a sequence of approximations. Assuming that the precise value of $x$ is available, the sequence of approximations is generated by:

\begin{align*}
\end{align*}

In what follows, we prove a theorem for approximating $\sqrt{x}$ by the Newton-Raphson method when an approximation $(m,n,p)$ of $x$ is available.

Let $x$ be a real number represented by $(m,n,p)$ such that:

\begin{align*}
\frac{m}{n}=0.b_1\ldots b_k\times 2^a
\end{align*}

where $b_i\in\{0,1\}$ for $1\leq i\leq k$ and $b_1=1, a\in \integer$. Then $\sqrt{x}$ can be represented by $(m',n',p-4N)$ where $N=\lceil \sqrtexp\rceil+1$ and $\frac{m'}{n'}$ is the $N$-th term of the following sequence:

\begin{align*}
y_{n+1}=\frac{y_n^2+\frac{m}{n}}{2y_n},~y_0=2^{\lceil\frac{a}{2} \rceil +1}
\end{align*}

From Definition <ref> we write:

\begin{align*}
\frac{m}{n}(1-\frac{1}{2^{p}})<x<\frac{m}{n}(1+\frac{1}{2^{p}})
\end{align*}

To prove the theorem, it suffices to show:

\begin{align}
\end{align}

We rewrite the left hand side of inequality (<ref>):

\begin{align}
|\sqrt{x}-y_N|=|\sqrt{x}-z_N+z_N-y_N|\leq |\sqrt{x}-z_N|+|z_N-y_N|\label{eqn:sqrt_rewrite}
\end{align}

In inequality (<ref>), $z_N$ is the $N$-th term of the following sequence:

\begin{align*}
z_{n+1}=\frac{z_n^2+x}{2z_n},~z_0=2^{\lceil\frac{a}{2} \rceil+1}
\end{align*}

To prove inequality (<ref>), it suffices to show that the following inequalities hold:

\begin{align}
\end{align}

Proof for inequality (<ref>):
First, we show that $y_n>\sqrt{\frac{x}{2}}$ for all $n\in \nat$:

\begin{align*}
y_n-\sqrt{\frac{x}{2}}=\frac{y_{n-1}^2+\frac{m}{n}}{2y_{n-1}}-\sqrt{\frac{x}{2}}&>\frac{y_{n-1}^2+\frac{m}{n}}{2y_{n-1}}-\sqrt{\frac{m}{2n}(1+\frac{1}{2^p})}\geq \frac{y_{n-1}^2+\frac{m}{n}}{2y_{n-1}}-\sqrt{\frac{m}{n}}\\
&=\frac{y_{n-1}^2+\frac{m}{n}-2y_{n-1}\sqrt{\frac{m}{n}}}{2y_{n-1}}=\frac{(y_{n-1}-\sqrt{\frac{m}{n}})^2}{2y_{n-1}}\geq 0
\end{align*}

Observe that as $y_n>\sqrt{\frac{x}{2}}$, inequality (<ref>) is valid, if the following inequality holds:

\begin{align}
\end{align}

To prove inequality (<ref>), we find $N$ such that:

\begin{align}
\end{align}

We calculate the quantity  $\frac{z_N+\sqrt{x}}{z_N-\sqrt{x}}$:

\begin{align}
\frac{z_N+\sqrt{x}}{z_N-\sqrt{x}}=&\frac{\frac{z_{N-1}^2+x}{2z_{N-1}}+\sqrt{x}}{\frac{z_{N-1}^2+x}{2z_{N-1}}-\sqrt{x}}=\frac{z_{N-1}^2+x+2z_{N-1}\sqrt{x}}{z_{N-1}^2+x-2z_{N-1}\sqrt{x}}\nonumber\\
\end{align}

Suppose equality (<ref>) holds for $N$. We can rewrite the right hand side of equality (<ref>) as follows:

\begin{align}
\end{align}

The index $N$ that guarantees the required precision of equality (<ref>) can be calculated from equality (<ref>):

\begin{align}
&2^{N}\log_2 (\frac{z_{0}+\sqrt{x}}{z_{0}-\sqrt{x}})=\log_2 (2^{p+3}+1)\nonumber\\
&N=\log_2 (\log_2 (2^{p+3}+1))- \log_2 (\log_2 (\frac{z_{0}+\sqrt{x}}{z_{0}-\sqrt{x}}))\label{eqn:nr_iteration_rewrite}
\end{align}

To guarantee that $N$ is well-defined, we show that $z_0>x$. To prove the inequality, we use the assumptions $\frac{m}{n}=0.b_1\ldots b_k\times 2^a, b_i\in \{0,1\} \text{ for }i=1,\ldots , k$ and $b_0=1, z_0=2^{\lceil \frac{a}{2}\rceil +1}$:

\begin{align*}
\sqrt{x}<\sqrt{\frac{m}{n}(1+\frac{1}{2^p})}\leq\sqrt{0.b_1\ldots b_k}\times 2^{\frac{a}{2}}\times \sqrt{2}<2^{\frac{a}{2}}\times \sqrt{2}<z_0
\end{align*}

To obtain an estimation for $N$ from equality (<ref>), we calculate an upper bound for $\frac{z_0}{\sqrt{x}}$:

\begin{align*}
\frac{z_0}{\sqrt{x}}<\frac{2^{\lceil \frac{a}{2} \rceil+1}}{\sqrt{\frac{m}{n}(1-\frac{1}{2^{p}})}}<\frac{2^{\frac{a+1}{2} +1}}{\sqrt{0.b_1\ldots b_k}\times 2^{\frac{a}{2}}\times \sqrt{\frac{1}{2}}}\leq \frac{2\sqrt{2}}{\sqrt[4]{2}\times\sqrt{\frac{1}{2}}}=\frac{4}{\sqrt[4]{2}}
\end{align*}

In the worst case, the initial approximation $z_0$ differs from $\sqrt{x}$ by a factor $\frac{4}{\sqrt[4]{2}}$. We use this estimation in equality (<ref>) to calculate the number of iterations for the Newton-Raphson method:

\begin{align}
N=\log_2 (\log_2 (2^{p+3}+1))-& \log_2 (\log_2 (\frac{\frac{4}{\sqrt[4]{2}}+1}{\frac{4}{\sqrt[4]{2}}-1}))\nonumber\\
<& \lceil\sqrtexp \rceil+1\label{eqn:n_upper}
\end{align}

Proof for inequality (<ref>): To prove inequality (<ref>), we consider the calculations in $z_N=\frac{z_{N-1}^2+x}{2z_{N-1}}$ and estimate the amount of error that we commit in the approximation $y_N=\frac{y_{N-1}^2+\frac{m}{n}}{2y_{N-1}}$.

Let $P(k)$ denote the amount of precision that we lose when we approximate $z_k$ by $y_k$. Thus, we lose $P(N-1)$ units of precision if we approximate $z_{N-1}$ by $y_{N-1}$:

\begin{align*}
\end{align*}

The precision is reduced by $2$ units when $z_{N-1}^2$ is approximated by $y_{N-1}^2$ (see Theorem <ref>):

\begin{align*}
\end{align*}

We approximate $z_{N-1}^2+x$ by $y_{N-1}^2+\frac{m}{n}$. We do not lose precision in this approximation (see Theorem <ref>.<ref>):

\begin{align}
\end{align}

Given the approximation $y_{N-1}$ of $z_{N-1}$, one unit of precision is lost in the approximation of $\frac{1}{z_{N-1}}$ (see Theorem <ref>):

\begin{align}
|\frac{1}{z_{N-1}}-\frac{1}{y_{N-1}}|&<\frac{1}{2^{p-P(N-1)-1}}|\frac{1}{y_{N-1}}|\nonumber \\
\end{align}

Finally, we approximate $\frac{z_{N-1}^2+x}{2z_{N-1}}$ based on the approximations described in inequality (<ref>) and (<ref>) (see Theorem <ref>):

\begin{align}
\end{align}

From inequality (<ref>) we obtain the following recursive formula:

\begin{align*}
\end{align*}

Since $y_0=z_0=2^{\lceil\frac{a}{2} \rceil +1}$, we lose $P(N)=P(0)+4N=4N$ units of precision in our approximation of $z_N$. We apply the number of iterations calculated in inequality (<ref>) and obtain:

\begin{align*}
P(N)=4N<4\lceil\sqrtexp \rceil+4
\end{align*}

Square Root
[1]
$\expr$ has the shape $\sqrt{x}$

Choose $p_x$ such that $p_x\geq p +4\lceil\log_2\big(\log_2(2^{p_x+3}+1)\big) \rceil+4$
$N\gets \lceil\log_2\big(\log_2(2^{p_x+3}+1)\big) \rceil+1$
$\frac{m}{n}\gets \Call{Compute}{x,p_x}$
“Undefined operation”
$\frac{m}{n}$ can be represented as $0.b_1\ldots b_k\times 2^a$
$b_i\in\{0,1\}$ for $1\leq i\leq k, b_1=1$ and $a\in\integer$
$a\gets \lfloor\log_2(\frac{m}{n}) \rfloor +1$
$y_0\gets 2^{\lceil\frac{a}{2} \rceil +1}$
$i=1$ to $N$
$y_{i}\gets \frac{y_{i-1}^2+\frac{m}{n}}{2y_{i-1}}$
$\frac{m'}{n'}\gets y_N$

Algorithm <ref> applies Theorem <ref> to approximate $\sqrt{x}$ with precision $p$. As indicated in Theorem <ref>, loss of precision in the square root can be estimated independently of the argument $x$ and hence Algorithm <ref> approximates $\sqrt{x}$ in one pass. We apply perturbation analysis on $f(x)=\sqrt{x}$ to show this:

\begin{align*}
\end{align*}

The quantity $|\frac{xf'(x)}{f(x)}|$ is a small constant. Thus, in our top down approach, we can approximate $\sqrt{x}$ without iterative computations.

§ APPROXIMATING TRANSCENDENTAL FUNCTIONS BY RIEMANN SUMS

In this section we introduce approximations for $e^x$ (Section <ref>), $\ln(x)$ (Section <ref>), and $\arctan(x)$ (Section <ref>). We use Riemann sums to approximate these functions.

The Compute$(expr,p)$ function introduced in Section <ref> will be extended to allow the approximation of $e^x, \ln(x)$ and $\arctan(x)$ with a given precision $p$. Perturbation analysis will also be used to identify computational problems in which iterative computations are unavoidable.

§.§ Exponential Function

To approximate the exponential function by Riemann sums, we first provide a simple approximation for $e^x$ where we assume that $x$ is precise. Then we extend this calculation to approximate $e^x$ where $x$ is represented by $(m,n,p)$.

axis lines=middle,
every extra x tick/.style=
xticklabel style=
extra x ticks =-0.9,0.9,1,
extra x tick labels=$-x$,$x$,$t$,
axis on top
[smooth, thick,domain=-0.95:0.95,samples=40]e^x;
]e^x ;

[Right & Left Riemann Sums]
axis lines=middle,
every extra x tick/.style=
xticklabel style=
extra x ticks =-0.9,0.9,1,
extra x tick labels=$-x$,$x$,$t$,
axis on top
[draw=green, fill=green!10, ybar interval, samples=20, domain=0.9:-0.9]
[draw=red, fill=red!10, ybar interval, samples=20, domain=-0.9:0.9]
[smooth, thick,domain=-0.95:0.95,samples=40]e^x;

Approximating $e^x$

Suppose $x>0$. To calculate $e^x$ we consider the curve $y=e^t$ and calculate the area enclosed by this curve and the t-axis between $t=-x$ and $t=x$ as follows (see Fig. <ref>):

\begin{align*}
\int_{-x}^{x} e^t dt=e^x-e^{-x}
\end{align*}

We use Riemann sums to approximate this area; Fig. <ref> shows two approximations from above and below using rectangles. Thus, we get the following inequalities for $N$ rectangles:

\begin{align}
\sum_{i=0}^{N-1}\frac{2x}{N}e^{-x+\frac{2x}{N}i}\leq e^x-e^{-x}\leq \sum_{i=1}^{N}\frac{2x}{N}e^{-x+\frac{2x}{N}i}\label{eqn:exp_sum}
\end{align}

We rewrite inequality (<ref>) as follows:

\begin{align}
\sum_{i=0}^{N-1}(\frac{2x}{N})e^{\frac{2x}{N}i}&\leq e^{2x}-1\leq \sum_{i=1}^{N}(\frac{2x}{N})e^{\frac{2x}{N}i}\nonumber\\
(\frac{2x}{N})(\frac{e^{2x}-1}{e^{\frac{2x}{N}}-1})&\leq e^{2x}-1\leq (\frac{2x}{N})(\frac{e^{\frac{2x}{N}}(e^{2x}-1)}{e^{\frac{2x}{N}}-1})\nonumber\\
\frac{2x}{N}&\leq e^{\frac{2x}{N}}-1 \leq \frac{2x}{N}e^{\frac{2x}{N}}\label{eqn:exp_sum_rewrite}
\end{align}

We assume that $N>2x$ and calculate an upper bound and a lower bound for $e^x$ from inequality (<ref>):

\begin{align}
(1+\frac{2x}{N})^{\frac{N}{2}} \leq e^x\leq (\frac{N}{N-2x})^{\frac{N}{2}}\label{eqn:exp_approx}
\end{align}

We can estimate the precision of the approximations calculated in inequality (<ref>). For example, we can approximate $e^x$ by $(\frac{N}{N-2x})^{\frac{N}{2}}$ and the absolute error of this approximation can be calculated as follows:

\begin{align}
|(\frac{N}{N-2x})^{\frac{N}{2}}-e^x|&\leq |(\frac{N}{N-2x})^{\frac{N}{2}}-(1+\frac{2x}{N})^{\frac{N}{2}}|\nonumber\\
&\leq \frac{4x^2}{N(N-2x)}\sum_{i=0}^{\frac{N}{2}-1}(\frac{N}{N-2x})^{\frac{N}{2}-1}=\frac{2x^2}{N}(\frac{N}{N-2x})^{\frac{N}{2}}\label{eqn:exp_riemann_prec}
\end{align}

For the last inequality we apply $\frac{N}{N-2x} \geq \frac{N+2x}{N}$.

In the discussion above, we have treated $x$ as a precise value. In the following theorem, we extend this calculation and describe an approximation of $e^x$ that relies on a representation $(m,n,p)$ of $x$. To simplify our approximations, we first assume that $|\frac{m}{n}|<1$. Afterwards, we extend our approximations to an arbitrary $(m,n,p)$.

Let $x$ be a real number represented by $(m,n,p)$ and $|\frac{m}{n}|<1$. Suppose $N$ is a natural number such that $N>2^{\frac{p+11}{3}}$. The value of $e^x$ can be approximated as follows:

If $0<\frac{m}{n}<1$ then $e^x$ can be represented by $(m',n',p-2\lceil\log_2(\frac{N}{2}) \rceil-3)$ where $\frac{m'}{n'}=(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}$.

If $-1<\frac{m}{n}<0$ then $e^x$ can be represented by $(m',n',p-2\lceil\log_2(\frac{N}{2}) \rceil-4)$ where $\frac{m'}{n'}=\frac{1}{(\frac{N}{N+\frac{2m}{n}})^{\frac{N}{2}}}$.

* Suppose $0<\frac{m}{n}<1$. Since $x$ is represented by $(m,n,p)$, we have:

\begin{align}
0<\frac{m}{n}(1-\frac{1}{2^p})<x<\frac{m}{n}(1+\frac{1}{2^p})\leq \frac{2m}{n}<2\label{eqn:exp_bound}
\end{align}

To prove the theorem, it suffices to show that:

\begin{align}
|e^x-(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|<\frac{1}{2^{p-2\lceil \log_2(\frac{N}{2})\rceil-3}}|(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|\label{eqn:exp}
\end{align}

We rewrite the left hand side of inequality (<ref>) as follows:

\begin{align*}
&\leq |e^x-(\frac{N}{N-2x})^{\frac{N}{2}}|+|(\frac{N}{N-2x})^{\frac{N}{2}}-(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|
\end{align*}

To prove inequality (<ref>), it suffices to show that the following inequalities hold:

\begin{align}
&|e^x-(\frac{N}{N-2x})^{\frac{N}{2}}|<\frac{1}{2^{p-2\lceil \log_2(\frac{N}{2})\rceil-2}}|(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|\label{eqn:exp_part_1} \\
&|(\frac{N}{N-2x})^{\frac{N}{2}}-(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|<\frac{1}{2^{p-2\lceil\log_2(\frac{N}{2}) \rceil-2}}|(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|\label{eqn:exp_part_2}
\end{align}

Proof for inequality (<ref>):

From inequality (<ref>), we obtain an upper-bound for the left hand side of inequality (<ref>):

\begin{align}
|e^x-(\frac{N}{N-2x})^{\frac{N}{2}}|\leq \frac{2x^2}{N}(\frac{N}{N-2x})^{\frac{N}{2}}\label{eqn:exp_part1_upper}
\end{align}

To calculate an upper bound for the right hand side of inequality (<ref>), we consider the function $f(x)=\frac{2x^2}{N}(\frac{N}{N-2x})^{\frac{N}{2}}$ and calculate its derivative:

\begin{align*}
\end{align*}

From inequality (<ref>) we obtain $x\in(0,2)$. We choose:

\begin{align}
N> 4>2x\label{eqn:exp_bound1}
\end{align}

to ensure that $f(x)$ is increasing for $x\in (0,2)$, i.e., $f'(x)>0$. We rewrite inequality (<ref>) as follows:

\begin{align}
|e^x-(\frac{N}{N-2x})^{\frac{N}{2}}|\leq f(x) \leq f(2)=(\frac{8}{N})(\frac{N}{N-4})^{\frac{N}{2}}\label{eqn:exp_part_1_upper1}
\end{align}

We calculate a lower bound for the right hand side of inequality (<ref>) as follows:

\begin{align}
\frac{1}{2^{p-2\lceil\log_2(\frac{N}{2}) \rceil-2}}|(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|>\frac{1}{2^{p-2\log_2(\frac{N}{2})-2}}\label{eqn:exp_part_1_lower}
\end{align}

Thus, to prove inequality (<ref>) it suffices to show that the following inequality holds (see inequality (<ref>), (<ref>)):

\begin{align*}
\end{align*}

This is equivalent to the following:

\begin{align}
\end{align}

We choose $N>8$ and apply Proposition <ref> (see <ref>) to obtain an upper bound for $\log_2(1+\frac{4}{N-4})$:

\begin{align}
\log_2(1+\frac{4}{N-4})<\frac{4}{(N-4)\ln(2)}<\frac{8}{N-4}\label{eqn:exp_n_upper}
\end{align}

Based on inequality (<ref>),(<ref>), it is sufficient to find an $N>8$ satisfying:

\begin{align}
\end{align}

Inequality (<ref>) is equivalent to the following:

\begin{align*}
\end{align*}

From $N>8$, we conclude $\frac{16}{N-4}<4$. Thus, we choose $N$ such that $N>\max(2^{\frac{p+11}{3}},8)=2^{\frac{p+11}{3}}$.

Proof for inequality (<ref>): To prove the inequality, we estimate the amount of precision that is lost when we approximate $(\frac{N}{N-2x})^{\frac{N}{2}}$ by $(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}$.

The number $x$ is represented by $(m,n,p)$. Thus, we have:

\begin{align*}
\end{align*}

Since $N>8$, one unit of precision is lost when we approximate $N-2x$ by $N-\frac{2m}{n}$ (see Theorem <ref>.<ref>):

\begin{align*}
\end{align*}

Approximating $\frac{1}{N-2x}$ by $\frac{1}{N-\frac{2m}{n}}$ reduces the precision by one unit (see Theorem <ref>):

\begin{align*}
\end{align*}

Finally, approximating $(\frac{N}{N-2x})^{\frac{N}{2}}$ reduces the precision by $2\lceil \log_2(\frac{N}{2}) \rceil$ units (see Lemma <ref> in <ref>):

\begin{align*}
|(\frac{N}{N-2x})^{\frac{N}{2}}-(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|<\frac{1}{2^{p-2\lceil\log_2(\frac{N}{2}) \rceil -2}}|(\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}|
\end{align*}

* Suppose $\frac{m}{n}<0$. We use the following identity to calculate $e^x$:

\begin{align*}
\end{align*}

We represent $-x$ by $(-m,n,p)$ (see Theorem <ref>). Then, we apply the first part of the theorem and Theorem <ref> to approximate $e^{-x}$ and $\frac{1}{e^{-x}}$, respectively.

In what follows, we extend the approximations of Theorem <ref> and calculate the exponential function for $x$ represented by $(m,n,p)$ where $|\frac{m}{n}|\geq 1$.

Let $x$ be a real number represented by $(m,n,p)$ and $|\frac{m}{n}|\geq 1$. Suppose $k$ and $N$ are natural numbers such that:

\begin{align*}
\end{align*}

The value of $e^x$ can be approximated as follows:

If $\frac{m}{n}>0$ then $e^x$ can be represented by $(m',n',p-2\lceil\log_2(\frac{N}{2}) \rceil-2k-3)$ where $\frac{m'}{n'}=(\frac{N}{N-\frac{2m}{n}})^{N\cdot 2^{k-1}}$.

If $\frac{m}{n}<0$ then $e^x$ can be represented by $(m',n',p-2\lceil\log_2(\frac{N}{2}) \rceil-2k-4)$ where $\frac{m'}{n'}=\frac{1}{(\frac{N}{N+\frac{2m}{n}})^{N\cdot 2^{k-1}}}$.

Since $|\frac{m}{n}|\geq1$, we choose $k\in\nat$ such that $|\frac{m}{2^kn}|<1$. We use the following identity to calculate $e^x$:

\begin{align}
e^x=(e^{\frac{x}{2^k}})^{2^k} \label{eqn:exp_range}
\end{align}

We approximate $\frac{x}{2^k}$ by $\frac{m}{2^kn}$. Since $2^k$ is a constant, we do not lose precision in this approximation. We apply Theorem <ref> to approximate $e^\frac{x}{2^k}$. This approximation reduces the precision by:

* $2\lceil\log_2(\frac{N}{2}) \rceil +3$ units, if $0<\frac{m}{2^kn}<1$;

* $2\lceil\log_2(\frac{N}{2}) \rceil +4$ units, if $-1<\frac{m}{2^kn}<0$.

Suppose $\frac{m'}{n'}$ is the approximation obtained for $e^{\frac{x}{2^k}}$ from Theorem <ref>. We approximate $(e^{\frac{x}{2^k}})^{2^k}$ by $(\frac{m'}{n'})^{2^k}$; we lose $2k$ units of precision in this calculation (see Lemma <ref> in <ref>).

Exponential Function
[1]
$\expr$ has the shape $e^x$
Choose $N$ such that $N>2^{\frac{p+11}{3}}$
$p_x\gets p+2\lceil\log_2(\frac{N}{2}) \rceil+4$

$\frac{m}{n}\gets \Call{Compute}{x,p_x}$
$0<\frac{m}{n}<1$ Theorem <ref>.<ref>
$\frac{m'}{n'}\gets (\frac{N}{N-\frac{2m}{n}})^{\frac{N}{2}}$
$-1<\frac{m}{n}<0$ Theorem <ref>.<ref>
$\frac{m'}{n'}\gets \frac{1}{(\frac{N}{N+\frac{2m}{n}})^{\frac{N}{2}}}$
Choose $k\in\nat$ such that $|\frac{m}{2^kn}|<1$
$(\frac{m}{n}>0) \wedge$ Theorem <ref>.<ref>

elseif$(p_x-2\lceil\log_2(\frac{N}{2}) \rceil-2k-3\geq p)$
$\frac{m'}{n'}\gets (\frac{N}{N-\frac{2m}{n}})^{N\cdot 2^{k-1}}$
$(\frac{m}{n}<0) \wedge$ Theorem <ref>.<ref>

elseelse if$(p_x-2\lceil\log_2(\frac{N}{2}) \rceil-2k-4\geq p)$
$\frac{m'}{n'}\gets \frac{1}{(\frac{N}{N+\frac{2m}{n}})^{N\cdot 2^{k-1}}}$
$p_x\gets p_x +1$

Algorithm <ref> implements the approximations described by Theorem <ref> and <ref> to calculate $e^x$ with arbitrary precision. Observe that when $|\frac{m}{n}|<1$, $e^x$ can be approximated in one pass. However, when $|\frac{m}{n}|\geq 1$, loss of precision depends on the magnitude of $|\frac{m}{n}|$. Thus, recomputing $x$ with higher precisions might be necessary to compensate for the  loss of precision caused by applying equality (<ref>) (see Line <ref>,<ref> in Algorithm <ref>).

To confirm that iterative computations are unavoidable, we apply perturbation analysis on $f(x)=e^x$:

\begin{align*}
\end{align*}

The quantity $|x|$ can become arbitrarily large and hence approximating $e^x$ with precision $p$ in one pass is not always possible.

§.§ Natural Logarithm

In this section, we first discuss an approximation for $\ln(x)$ based on Riemann sums where we assume that $x$ is precise. Then we extend this calculation to approximate $\ln(x)$ where $x$ is represented by $(m,n,p)$.

axis lines=middle,
every extra x tick/.style=
xticklabel style=
extra x ticks =1,3,3.7,
extra x tick labels=$1$,$x$,$t$,
axis on top
[smooth, thick,domain=0.3:3.7,samples=40]x^-1;
]x^-1 ;

[Right & Left Riemann Sums]
axis lines=middle,
every extra x tick/.style=
xticklabel style=
extra x ticks =1,3,3.7,
extra x tick labels=$1$,$x$,$t$,
axis on top
[draw=red, fill=red!10, ybar interval, samples=9, domain=1:3]
[draw=green, fill=green!10, ybar interval, samples=9, domain=3:1]
[smooth, thick,domain=0.3:3.7,samples=40]x^-1;

Approximating $\ln(x)$

Suppose $x>1$ is a real number and we want to approximate $\ln(x)$. We consider the curve $y=\frac{1}{t}$ (see Fig. <ref>) and calculate the area enclosed by this curve and the t-axis between $t=1$ and $t=x$. This area can be calculated as follows:

\begin{align*}
\int_1^{x} \frac{dt}{t}=\ln(x)
\end{align*}

We use Riemann sums to approximate this area; Fig. <ref> shows how the area can be approximated from below and above using rectangles. Thus, we get the following inequalities for $N$ rectangles:

\begin{align*}
\lnsuml{x}{N} \leq \ln(x) \leq \lnsumu{x}{N}
\end{align*}

This gives us an upper bound and a lower bound for $\ln(x)$ and by increasing $N$ we get more precise approximations.

We can estimate the precision of our approximations. For instance, if we approximate $\ln(x)$ by $\lnsumu{x}{N}$ the absolute error can be estimated as follows:

\begin{align}
&\phantom{\frac{x-1}{N}\sum}= \frac{x-1}{N}(1-\frac{1}{x})=\frac{(x-1)^2}{Nx}\label{eqn:ln_riemann_prec}
\end{align}

Up to this point, we have assumed that the precise value of $x$ is available. In what follows, we formulate a theorem to describe an approximation of $\ln(x)$ based on a representation $(m,n,p)$ of $x$.

Let $x$ be a real number represented by $(m,n,p)$ such that $p\geq 1$.

* If $\frac{m}{n}>1$, then $\ln(x)$ can be represented by $(m',n',p-j-4)$ where $\frac{m'}{n'}=\lnsumu{(\frac{m}{n})}{N}$, $N=\lceil 2^{p-2}\frac{(\frac{m}{n})^2}{\frac{m}{n}-1}\rceil$, and $j$ is the smallest natural number such that $j\geq \log_2(\frac{1+\frac{n}{m}}{1-\frac{n}{m}})$ holds.

* If $0<\frac{m}{n}<1$, then $\ln(x)$ can be represented by $(m',n',p-j-5)$ where $\frac{m'}{n'}=-\lnsumu{(\frac{n}{m})}{N}$, $N=\lceil 2^{p-2}\frac{(\frac{n}{m})^2}{\frac{n}{m}-1}\rceil$, and $j$ is the smallest natural number such that $j\geq \log_2(\frac{1+\frac{m}{n}}{1-\frac{m}{n}})$ holds.

* Suppose $\frac{m}{n}>1$. The number $x$ is represented by $(m,n,p)$ and hence we can write:

\begin{align}
\frac{1}{2}<\frac{m}{2n}\leq \frac{m}{n}(1-\frac{1}{2^p})&<x<\frac{m}{n}(1+\frac{1}{2^p})\leq \frac{3m}{2n} \label{eqn:ln_x_bound}
\end{align}

To prove the theorem, we need to prove the following inequality:

\begin{align}
\end{align}

We rewrite the left hand side of inequality (<ref>) as follows:

\begin{align*}
+&\lnsumu{x}{N}-\lnsumu{\frac{m}{n}}{N}|\leq \\
\end{align*}

To prove inequality (<ref>), it suffices to show that the following inequalities hold:

\begin{align}
&|\ln(x)-\lnsumu{x}{N}|<\frac{1}{2^{p-j-3}}|\lnsumu{\frac{m}{n}}{N}|\label{eqn:ln_part_1} \\
\end{align}

Proof for inequality (<ref>):
We use inequality (<ref>) and calculate an upper bound for the left hand side of inequality (<ref>):

\begin{align}
|\ln(x)-\lnsumu{x}{N}|\leq \frac{(x-1)^2}{Nx}\label{eqn:ln_part_1_upper}
\end{align}

To calculate an upper bound for $\frac{(x-1)^2}{Nx}$, we consider the function $f(x)=\frac{(x-1)^2}{Nx}$ and calculate its derivative:

\begin{align*}
\end{align*}

From inequality (<ref>) we conclude that $x\in [\frac{1}{2},\frac{3m}{2n}]$. The function $f(x)$ is decreasing ($f'(x)\leq 0$) in the interval $[\frac{1}{2},1]$ and increasing ($f'(x)\geq 0$) in the interval $[1,\frac{3m}{2n}]$. Thus, the maximum of $f(x)$ for $x\in [\frac{1}{2},\frac{3m}{2n}]$ is $\max (f(\frac{1}{2}),f(\frac{3m}{2n}))$. We use this to rewrite inequality (<ref>):

\begin{align}
|\ln(x)-\lnsumu{x}{N}|&\leq f(x)\leq \max (f(\frac{1}{2}),f(\frac{3m}{2n}))\nonumber \\
&\leq \max (\frac{1}{2N},\frac{(\frac{3m}{2n})^2}{N(\frac{3m}{2n})})=\frac{3m}{2Nn}\label{eqn:ln_part_1_upper_1}
\end{align}

We also calculate a lower bound for the right hand side of inequality (<ref>):

\begin{align}
\frac{1}{2^{p-j-3}}|\lnsumu{\frac{m}{n}}{N}|>&\frac{1}{2^{p-3}}|\frac{(\frac{m}{n}-1)}{N}\sum_{i=0}^{N-1}\frac{1}{1+(\frac{N-1}{N})(\frac{m}{n}-1)}|\nonumber\\
\end{align}

To show that inequality (<ref>) holds, it suffices to prove the following inequality (see inequality (<ref>),(<ref>)):

\begin{align*}
\frac{3m}{2Nn}<\frac{(\frac{m}{n}-1)}{2^{p-3}(\frac{m}{n})}
\end{align*}

Thus, it suffices to choose $N\geq 2^{p-2}\frac{(\frac{m}{n})^2}{\frac{m}{n}-1}$.

Proof for inequality (<ref>): To prove inequality (<ref>) we consider the calculations in $\lnsumu{x}{N}$ and estimate the amount of precision that we lose in the approximation $\lnsumu{\frac{m}{n}}{N}$.

From Theorem <ref>.<ref> we conclude that $j$ units of precision is lost by subtracting $1$ from $x$ where $j\geq \log_2(\frac{1+\frac{n}{m}}{1-\frac{n}{m}})$. We obtain the following inequalities:

\begin{gather*}
\end{gather*}

The approximation $\frac{i}{N}(\frac{m}{n}-1)$ of $\frac{i}{N}(x-1)$ is positive. Hence, we do not lose precision by adding $\frac{i}{N}(x-1)$ and $1$ (see Theorem <ref>.<ref>).

\begin{align*}
\end{align*}

By approximating the inverse of $1+\frac{i}{N}(x-1)$, we lose $1$ unit of precision (see Theorem <ref>):

\begin{align*}
\end{align*}

We lose $2$ units of precision by multiplying $\frac{x-1}{N}$ and $\frac{1}{1+\frac{i}{N}(x-1)}$ (see Theorem <ref>):

\begin{align*}
\end{align*}

Since the numbers $\frac{x-1}{N}.\frac{1}{1+\frac{i}{N}(x-1)}$ are approximated by the positive numbers $\frac{\frac{m}{n}-1}{N}.\frac{1}{1+\frac{i}{N}(\frac{m}{n}-1)}$, calculating the summation $\lnsumu{x}{N}$ does not affect the precision:

\begin{align*}
\end{align*}

* Suppose $0<\frac{m}{n}<1$. We use the following identity to approximate $\ln(x)$:

\begin{align*}
\ln(x)=-\ln(\frac{1}{x})
\end{align*}

We approximate $\frac{1}{x}$ by $(n,m,p-1)$ (see Theorem <ref>). Then, we apply the first part of the theorem and Theorem <ref> to approximate $-\ln(\frac{1}{x})$.

Natural Logarithm
[1]
$\expr$ has the shape $\ln(x)$
$p_x \gets p+5$

$\frac{m}{n} \gets \Call{Compute}{x,p_x}$
$\textit{arg}\gets \frac{m}{n}$
$\ell \gets 4$
$\textit{arg}\gets \frac{n}{m}$
$\ell \gets 5$
“Undefined operation”

$j\gets \lceil\log_2(\frac{1+\frac{1}{\textit{arg}}}{1-\frac{1}{\textit{arg}}}) \rceil$
$p_x-j-\ell\geq p$
$N\gets \lceil 2^{p_x-2}\frac{\textit{arg}^2}{\textit{arg}-1}\rceil$
$\frac{m'}{n'}\gets \lnsumu{\textit{arg}}{N}$
$p_x\gets p_x+1$

Algorithm <ref> applies Theorem <ref> to approximate $\ln(x)$ with arbitrary precision. Note that when the approximation $\frac{m}{n}$ is close to $1$ the amount of precision that we lose in the calculations depends on the magnitude of $\frac{m}{n}$. Loss of precision for $x\approx 1$ in Algorithm <ref> is due to our approximation formula, $\lnsumu{x}{N}$. We divide the interval between $1$ and $x$ into $N$ subintervals and approximate the area under the curve $f(t)=\frac{1}{t}$. The length of the interval $[1,x]$ is crucial in our approximation. Thus, recomputing $x$ with higher precisions is necessary when a significant amount of precision is lost in $x-1$ (see Line <ref>,<ref> in Algorithm <ref>).

To show that approximating $\ln(x)$ for $x\approx 1$ is fundamentally problematic we apply perturbation analysis  on $f(x)=\ln(x)$:

\begin{align*}
\end{align*}

When $x\approx 1$ the quantity $\ln(x)$ is a small value and hence committing a small error in the approximation of $x$ causes a significant error in calculating $\ln(x)$. Thus, for $x\approx 1$, iterative computations  are unavoidable.

§.§ Arctangent

We first introduce an approximation of $\arctan(x)$ using Riemann sums. The assumption is that the precise value of $x$ is available. Afterwards, we extend our calculations to introduce an approximation of $\arctan(x)$ based on a representation $(m,n,p)$ of $x$.

axis lines=middle,
every extra x tick/.style=
xticklabel style=
extra x ticks =3,3.5,
extra x tick labels=$x$,$t$,
axis on top
[smooth, thick,domain=0:3,samples=40]1/(1+x^2);
]1/(1+x^2) ;

[Right & Left Riemann Sums]
axis lines=middle,
every extra x tick/.style=
xticklabel style=
extra x ticks =3,3.5,
extra x tick labels=$x$,$t$,
axis on top

[draw=red, fill=red!10, ybar interval, samples=14, domain=0:3]
[draw=green, fill=green!10, ybar interval, samples=14, domain=3:0]

[smooth, thick,domain=0:3,samples=40]1/(1+x^2);

Approximating $\arctan(x)$

Suppose that $x>0$ is a real number and we want to approximate $\arctan(x)$. We consider the curve $y=\frac{1}{1+t^2}$; see Fig. <ref>. We calculate the area enclosed by this curve and the t-axis between $t=0$ and $t=x$:

\begin{align*}
\int_0^{x} \frac{dt}{1+t^2}=\arctan(x)
\end{align*}

We approximate this area using Riemann sums. Fig. <ref> shows approximations from above and below for the integral using rectangles. From Fig. <ref> we can derive the following inequalities for $N$ rectangles:

\begin{align*}
\arctansuml{x}{N} \leq \arctan(x) \leq \arctansumu{x}{N}
\end{align*}

We want to estimate the precision of our approximations. Suppose we approximate $\arctan(x)$ by $\arctansumu{x}{N}$. The absolute error of this approximation can be estimated as follows:

\begin{align}
|\arctansumu{x}{N}-\arctan(x)|&\leq |\arctansumu{x}{N}-\arctansuml{x}{N}|\nonumber\\
\end{align}

Up to this point, we have assumed that $x$ is precisely calculated. In what follows, we assume that $x$ is approximated by $(m,n,p)$. We extend the Riemann sum calculation to compute $\arctan(x)$ using the given approximation of $x$.

Let $x$ be a real number represented by $(m,n,p)$. Then $\arctan(x)$ can be represented by $(m',n',p-6)$ where $\frac{m'}{n'}=\arctansumu{(\frac{m}{n})}{N}$ and $N= 2^{p-2}\lceil(\frac{m}{n})^2\rceil$.

We consider two cases:

* Suppose $\frac{m}{n}>0$. Since $x$ is represented by $(m,n,p)$ we can write:

\begin{align}
\frac{m}{n}(1-\frac{1}{2^p})<x<\frac{m}{n}(1+\frac{1}{2^p})\leq \frac{2m}{n}\label{eqn:arctan_x_bound}
\end{align}

To prove the theorem, we should show that:

\begin{align}
\end{align}

We rewrite the left hand side of inequality (<ref>) as follows:

\begin{align*}
&\phantom{|\arctan(x)}+\arctansumu{x}{N}-\arctansumu{(\frac{m}{n})}{N}|\leq \\
\end{align*}

Thus, to prove inequality (<ref>), it suffices to prove the following inequalities:

\begin{align}
\end{align}

Proof for inequality (<ref>):
We first consider the left hand side of inequality (<ref>) and calculate an upper bound for it. From inequality (<ref>) we can write:

\begin{align}
|\arctan(x)-\arctansumu{x}{N}|\leq \frac{x^3}{N(1+x^2)}\label{eqn:arctan_part_1_upper}
\end{align}

To obtain an upper-bound for $\frac{x^3}{N(1+x^2)}$, we consider the function $f(x)=\frac{x^3}{N(1+x^2)}$ and calculate its derivative:

\begin{align*}
\end{align*}

Thus, $f(x)$ is an increasing function and its maximum occurs when $x$ gets its maximum value. Inequality (<ref>) implies that $\frac{2m}{n}$ is an upper bound for $x$ and hence we can rewrite inequality (<ref>) as follows:

\begin{align}
|\arctan(x)-\arctansumu{x}{N}|\leq \frac{x^3}{N(1+x^2)}&\leq \frac{8(\frac{m}{n})^3}{N(1+4(\frac{m}{n})^2)}\nonumber\\
\end{align}

We also calculate a lower bound for the right hand side of inequality (<ref>):

\begin{align}
\frac{1}{2^{p-5}}|\arctansumu{(\frac{m}{n})}{N}|>&(\frac{1}{2^{p-5}})\frac{(\frac{m}{n})}{N}\sum_{i=0}^{N-1}\frac{1}{1+(\frac{N-1}{N})^2(\frac{m}{n})^2}\nonumber\\
\end{align}

To prove inequality (<ref>), it suffices to show that the following inequality holds (see inequality (<ref>),(<ref>)):

\begin{align*}
\frac{8(\frac{m}{n})^3}{N(1+(\frac{m}{n})^2)}<\frac{(\frac{m}{n})}{2^{p-5}(1+(\frac{m}{n})^2)}
\end{align*}

We divide all the components by $\frac{\frac{m}{n}}{(1+(\frac{m}{n})^2)}$, obtaining:

\begin{align*}
\frac{8(\frac{m}{n})^2}{N}<\frac{1}{2^{p-5}}
\end{align*}

Thus, it suffices to take $N= 2^{p-2}\lceil(\frac{m}{n})^2\rceil$.

Proof for inequality (<ref>):
To prove the inequality, we consider the calculations in $\arctansumu{x}{N}$ and estimate the amount of error that we commit in the approximation $\arctansumu{(\frac{m}{n})}{N}$.

From Theorem <ref> and inequality (<ref>), we conclude that $2$ units of precision is lost by approximating $x^2$ by $(\frac{m}{n})^2$ and hence we obtain:

\begin{gather*}
\end{gather*}

The approximation $(\frac{i}{N})^2(\frac{m}{n})^2$ of $(\frac{i}{N})^2x^2$ is positive. Thus, we do not lose precision by adding $(\frac{i}{N})^2x^2$ and $1$ (see Theorem <ref>.<ref>):

\begin{align*}
\end{align*}

We lose $1$ unit of precision by approximating the inverse of $1+(\frac{i}{N})^2x^2$ (see Theorem <ref>):

\begin{align*}
\end{align*}

We lose $2$ units of precision in the approximation of $\frac{x}{N}.\frac{1}{1+(\frac{i}{N})^2x^2}$ (see Theorem <ref>):

\begin{align*}
\end{align*}

The approximations $\frac{\frac{m}{n}}{N}.\frac{1}{1+(\frac{i}{N})^2(\frac{m}{n})^2}$ are positive for $0\leq i\leq N-1$. Thus, we do not lose precision by calculating the summation $\arctansumu{x}{N}$ (see Theorem <ref>.<ref>):

\begin{align*}
&\phantom{|\arctansumu{x}{N}-\frac{(\frac{m}{n})}{N}\sum \frac{m}{n}\frac{m}{n}}<\frac{1}{2^{p-5}}|\arctansumu{(\frac{m}{n})}{N}|
\end{align*}

* Suppose $\frac{m}{n}<0$. We can write:

\begin{align*}
\frac{m}{n}(1+\frac{1}{2^p})<x<\frac{m}{n}(1-\frac{1}{2^p})
\end{align*}

Observe that $x<0$; we can use the following identity for $\arctan(x)$:

\begin{align*}
\arctan(x)=-\arctan(-x)
\end{align*}

where $-x>0$. We first approximate $\arctan(-x)$ using the first part of the proof. Afterwards, we apply Theorem <ref> to approximate $-\arctan(-x)$. Unary negation does not influence the precision.

[1]
$\expr$ has the shape $\arctan(x)$

$\frac{m}{n} \gets \Call{Compute}{x,p+6}$
$N\gets 2^{p+4}\lceil(\frac{m}{n})^2\rceil$
$\frac{m'}{n'}\gets \arctansumu{(\frac{m}{n})}{N}$

Algorithm <ref> applies Theorem <ref> to approximate $\arctan(x)$ with arbitrary precision. Note that  Theorem <ref> predicts the amount of precision that is lost by calculating $\arctan(x)$ independently of the argument $x$ and hence Algorithm <ref> calculates $\arctan(x)$ in one pass.

To confirm that $\arctan(x)$ is computable in one pass, we apply perturbation analysis and calculate the quantity $|\frac{xf'(x)}{f(x)}|$ for $f(x)=\arctan(x)$:

\begin{align*}
\end{align*}

Proposition <ref> (see <ref>) shows that $|\frac{x}{(1+x^2)\arctan(x)}|<1$; iterative computations are not required for approximating $\arctan(x)$.

§ APPROXIMATING TRANSCENDENTAL FUNCTIONS BY TAYLOR EXPANSIONS

In this section we first briefly discuss the basics of approximating functions using Taylor expansions. Afterwards, we use Taylor expansions to approximate $\sin(x)$ and $\cos(x)$.

Suppose $f: D\rightarrow R$ is a function  and $I=(a,b)\subseteq D$ such that:

* $f$ has $n$ continuous derivatives on $I$ (denoted by $f^{(i)}(x)$ for $1\leq i\leq n$);

* $f^{(n+1)}$ exists on $I$;

* $x_0\in I$.
Taylor's theorem states that for every $x\in I$
there is a number $c_x$ between $x$ and $x_0$ such that $f(x)=P_n(x)+R_n(x)$ where:

\begin{align}
\end{align}

The formula $R_n(x)$ is called the Lagrange form of the remainder [20].

In the remaining of this section, we discuss approximations for $\sin(x)$ and $\cos(x)$. Our approximations are based on the following Taylor expansions around the point $x_0=0$:

\begin{align}
\sin(x)=&\taylorsin{0}{\infty}{x} \label{eqn:taylor_sin}\\
\cos(x)=&\taylorcos{0}{\infty}{x}\label{eqn:taylor_cos}
\end{align}

For each function, we first describe an approximation that is applicable to the base interval $I=(-1,1)$ (see Section <ref>). Afterwards, we extend our calculations to the complete domain of the functions using range reduction identities (see Section <ref>). Proof of correctness for the base and general cases are provided. Moreover, the Compute$(expr,p)$ function is extended to approximate $\sin(x)$ and $\cos(x)$ with arbitrary precision. We use perturbation analysis to confirm the observations obtained from the approximations.

§.§ Approximating Functions in the Base Interval

§.§.§ Sine

In this section we use the Taylor expansion from equality (<ref>) to approximate $\sin(x)$. We assume that the input argument $x$ is represented by $(m,n,p)$ and $|\frac{m}{n}|<1$.

Let $x$ be a real number represented by $(m,n,p)$ such that $-1<\frac{m}{n}<1$. Then $\sin(x)$ can be represented by $(m',n',p-2\lceil\log_2(2N-1)\rceil-3)$ where $\frac{m'}{n'}=\taylorsin{0}{N}{(\frac{m}{n})}$ and $N\in\nat$ is an odd number such that $(\frac{5}{6})(\frac{(2N+2)!(2N-1)^2}{2^{2N}})>2^p$.

The real number $x$ is given with precision $p$:

\begin{align}
\frac{m}{n}-|\frac{m}{n}|\leq \frac{m}{n}-\frac{1}{2^p}|\frac{m}{n}|<x<\frac{m}{n}+\frac{1}{2^p}|\frac{m}{n}|\leq \frac{m}{n}+|\frac{m}{n}|\label{eqn:sin_bound}
\end{align}

To prove the theorem, we should show that the following inequality holds:

\begin{align}
\end{align}

We rewrite the left hand side of inequality (<ref>) as follows:

\begin{align*}
&|\sin(x)-\taylorsin{0}{N}{x}+\taylorsin{0}{N}{x}-\taylorsin{0}{N}{(\frac{m}{n})}|\leq \\
\end{align*}

Thus, we prove inequality (<ref>) by showing that the following inequalities hold:

\begin{align}
\end{align}

Proof for inequality (<ref>): Based on Taylor's theorem (see equality (<ref>)) and inequality (<ref>), we can rewrite the left hand side of inequality (<ref>) as follows:

\begin{align}
|\sin(x)-\taylorsin{0}{N}{x}|=|\frac{\sin^{(2N+2)}(c_x)x^{2N+2}}{(2N+2)!}|\leq \frac{|\frac{m}{n}|^{2N+2}2^{2N+2}}{(2N+2)!}\label{eqn:sin_part_1_left}
\end{align}

We choose $N=2k+1\geq 1$. Applying Proposition <ref> (see <ref>) we can rewrite the right hand side of inequality (<ref>) as follows:

\begin{align}
\frac{1}{2^{p-2\lceil\log_2(2N-1)\rceil-2}}|\taylorsin{0}{N}{(\frac{m}{n})}|=&\bigg(\frac{sgn(\frac{m}{n})}{2^{p-2\lceil\log_2(2N-1)\rceil-2}}\bigg)\nonumber\\
\cdot &\bigg(\taylorsin{0}{N}{(\frac{m}{n})}\bigg)\label{eqn:sin_part_1_right}
\end{align}

where $sgn$ is the sign function. To prove inequality (<ref>), it suffices to show that the following inequality holds (see (<ref>),(<ref>)):

\begin{align*}
\frac{|\frac{m}{n}|^{2N+2}2^{2N+2}}{(2N+2)!}<\frac{sgn(\frac{m}{n})}{2^{p-2\lceil\log_2(2N-1)\rceil-2}}\taylorsin{0}{N}{(\frac{m}{n})}
\end{align*}

This is equivalent to:

\begin{align}
\end{align}

The quantity $sgn(\frac{m}{n})\taylorsin{2}{N}{(\frac{m}{n})}$ is positive (see Proposition <ref> in <ref>). Moreover, the approximation $\frac{m}{n}$ satisfies $|\frac{m}{n}|<1$. To prove inequality (<ref>) it is sufficient to show:

\begin{align}
\frac{sgn(\frac{m}{n})}{2^{p-2\lceil\log_2(2N-1)\rceil-2}}(\frac{m}{n}-(\frac{1}{3!})(\frac{m}{n})^3)-\frac{|\frac{m}{n}|^{2N+2}2^{2N+2}}{(2N+2)!}>\nonumber\\
\frac{sgn(\frac{m}{n})(2N-1)^2}{2^{p-2}}(\frac{m}{n}-(\frac{1}{3!})(\frac{m}{n})^3)-\frac{|\frac{m}{n}|^{3}2^{2N+2}}{(2N+2)!}>0\label{eqn:sin_part_1_rewrite_1}
\end{align}

Inequality (<ref>) is equivalent to:

\begin{align*}
\end{align*}

Since $sgn(\frac{m}{n})(\frac{m}{n})>0$ and $(\frac{m}{n})^2<1$, we should choose an $N$ such that:

\begin{align*}
\frac{1}{6}+\frac{2^{p+2N-2\lceil\log_2(2N-1) \rceil}}{(2N+2)!(2N-1)^2}<1
\end{align*}

Thus, $N$ should satisfy $(\frac{5}{6})\frac{(2N+2)!(2N-1)^2}{2^{2N}}>2^p$.

Proof for inequality (<ref>): We show that given an approximation $\frac{m}{n}$ of $x$ with precision $p$ we can approximate $\taylorsin{0}{N}{x}$ by $\taylorsin{0}{N}{(\frac{m}{n})}$ with precision $p-2\lceil\log_2(2N-1)\rceil -2$.

The sign of the terms $\frac{(-1)^i}{(2i+1)!}x^{2i+1}$ alternates between positive and negative. Thus, adding two arbitrary terms with different signs from the expansion can significantly reduce the precision (see Theorem <ref>.<ref>). To avoid this, we first consider specific pairs of terms for which loss of precision due to addition is bounded. Then, we calculate the summation of these pairs. The following identity shows the way we calculate $\taylorsin{0}{N}{x}$:

\begin{align*}
\taylorsin{0}{N}{x}=\sum_{i=0}^{k}\frac{x^{4i+1}}{(4i+1)!}(1-\frac{x^2}{(4i+2)(4i+3)})
\end{align*}

Choosing $N=2k+1$ allows us to make pairs of terms.

The number $x$ is given with precision $p$. Thus, we can approximate $x^2$ with precision $p-2$ (see Theorem <ref>):

\begin{align*}
\end{align*}

and hence

\begin{align*}
\end{align*}

We approximate $1-\frac{x^2}{(4i+2)(4i+3)}$ by $1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)}$. Loss of precision in the approximation can be estimated by calculating the quantity $\log_2(\frac{1+\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)}}{1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)}})$ (see Theorem <ref>.<ref>):

\begin{align*}
\log_2(\frac{1+\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)}}{1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)}})\leq \log_2(\frac{1+\frac{1}{6}}{1-\frac{1}{6}})=\log_2(\frac{7}{5})<1
\end{align*}

Thus, we lose at most $1$ unit of precision:

\begin{align*}
\end{align*}

The powers $x^{4i+1}$ are approximated by $(\frac{m}{n})^{4i+1}$ for $0\leq i\leq k$ and $2\lceil\log_2(4i+1) \rceil$ units of precision is lost in this operation. Thus, in the worst case, the precision is reduced by $2\lceil\log_2(4k+1) \rceil$ units:

\begin{align*}
\end{align*}

Multiplying approximations of $x^{4i+1}$ by a constant does not influence the precision:

\begin{align*}
|\frac{x^{4i+1}}{(4i+1)!}-\frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}|<\frac{1}{2^{p-2\lceil\log_2(2N-1)\rceil}}| \frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}|
\end{align*}

The multiplication $\frac{x^{4i+1}}{(4i+1)!}(1-\frac{x^2}{(4i+2)(4i+3)})$ can be approximated by multiplying the approximations $\frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}$ and $1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)}$ (see Theorem <ref>):

\begin{align*}
&<\frac{1}{2^{p-\max\{2\lceil\log_2(2N-1) \rceil,3\}-2}}|\frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}(1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)})|\\
&=\frac{1}{2^{p-2\lceil\log_2(2N-1) \rceil-2}}|\frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}(1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)})|
\end{align*}

For the last equality, we use the assumptions $(\frac{5}{6})(\frac{(2N+2)!(N-1)^2}{2^{2N-2\lceil\log_2(2N-1) \rceil}})>2^p$ and $N=2k+1$; we can conclude that $N\geq 3$.

For $0\leq i\leq k$ the terms $\frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}(1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)})$ have the same sign which is determined by the sign of $\frac{m}{n}$ (see Proposition <ref> in <ref>). Thus, we can approximate $\sum_{i=0}^{k}\frac{x^{4i+1}}{(4i+1)!}(1-\frac{x^2}{(4i+2)(4i+3)})$ by $\sum_{i=0}^{k}\frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}(1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)})$ without losing precision (see Theorem <ref>.<ref>):

\begin{align*}
\frac{1}{2^{p-2\lceil\log_2(2N-1)\rceil -2}}|\sum_{i=0}^{k}\frac{(\frac{m}{n})^{4i+1}}{(4i+1)!}(1-\frac{(\frac{m}{n})^2}{(4i+2)(4i+3)})|
\end{align*}

Theorem <ref> provides a top-down approximation for $\sin(x)$ in the base interval. Loss of precision in this approximation is estimated independently of the argument $x$. To show that iterative calculations can be avoided in the base interval, we calculate $|\frac{xf'(x)}{f(x)}|$ for $f(x)=\sin(x)$:

\begin{align}
|\frac{xf'(x)}{f(x)}|=|\frac{x\cdot \cos(x)}{\sin(x)}|=|x\cdot \cot(x)|\label{eqn:pert_sin}
\end{align}

Proposition <ref> (see <ref>) shows that $|x\cdot \cot(x)|<1$ for $x\in (-1,1)$. Thus, we can approximate $\sin(x)$ in the base interval in one pass.

§.§.§ Cosine

In this section we introduce an approximation for $\cos(x)$ where $x$ is represented by $(m,n,p)$ and $|\frac{m}{n}|<1$. Our approximation is based on the Taylor expansion from equality (<ref>).

Let $x$ be a real number represented by $(m,n,p)$ such that $-1<\frac{m}{n}<1$. Then $\cos(x)$ can be represented by $(m',n',p-2\lceil \log_2(2N-2) \rceil -3)$ where $\frac{m'}{n'}=\taylorcos{0}{N}{(\frac{m}{n})}$ and $N\in\nat$ is an odd number such that $\frac{(2N+1)!(2N-2)^2}{2^{2N}}>2^p$.

The real number $x$ is represented by $(m,n,p)$ and hence we can write:

\begin{align}
\frac{m}{n}-|\frac{m}{n}|<\frac{m}{n}-\frac{1}{2^p}|\frac{m}{n}|<x<\frac{m}{n}+\frac{1}{2^p}|\frac{m}{n}|\leq \frac{m}{n}+|\frac{m}{n}|\label{eqn:cos_bound}
\end{align}

To prove the theorem, we need to show:

\begin{align}
|\cos(x)-\taylorcos{0}{N}{(\frac{m}{n})}|<\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -3}}|\taylorcos{0}{N}{(\frac{m}{n})}|\label{eqn:cos}
\end{align}

We rewrite the left hand side of inequality (<ref>) as follows:

\begin{align*}
\end{align*}

We prove inequality (<ref>) by showing that the following inequalities are valid:

\begin{gather}
|\cos(x)-\taylorcos{0}{N}{x}|<\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}|\taylorcos{0}{N}{(\frac{m}{n})}|\label{eqn:cos_part_1}\\
|\taylorcos{0}{N}{x}-\taylorcos{0}{N}{(\frac{m}{n})}|<\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}|\taylorcos{0}{N}{(\frac{m}{n})}|\label{eqn:cos_part_2}
\end{gather}

Proof for inequality (<ref>): We use Taylor's theorem (see equality (<ref>)) and the bounds calculated for $x$  in inequality (<ref>) to rewrite the left hand side of inequality (<ref>):

\begin{align}
|\cos(x)-\taylorcos{0}{N}{x}|=|\frac{\cos^{(2N+1)}(c_x)x^{2N+1}}{(2N+1)!}|\leq \frac{|\frac{m}{n}|^{2N+1}2^{2N+1}}{(2N+1)!}\label{eqn:cos_part_1_left}
\end{align}

We choose $N=2k+1\geq 1$ and apply Proposition <ref> (see <ref>) to rewrite the right hand side of inequality (<ref>):

\begin{align}
\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}|\taylorcos{0}{N}{(\frac{m}{n})}|=\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}\taylorcos{0}{N}{(\frac{m}{n})}\label{eqn:cos_part_1_right}
\end{align}

To show that inequality (<ref>) holds, it suffices to prove the following (see inequality (<ref>),(<ref>)):

\begin{align}
\frac{|\frac{m}{n}|^{2N+1}2^{2N+1}}{(2N+1)!}<\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}\taylorcos{0}{N}{(\frac{m}{n})}\label{eqn:cos_part_1_rewrite}
\end{align}

Inequality (<ref>) is equivalent to:

\begin{align*}
&\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}\taylorcos{0}{N}{(\frac{m}{n})}-\frac{|\frac{m}{n}|^{2N+1}2^{2N+1}}{(2N+1)!}=\\
&\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}(1-\frac{1}{2}(\frac{m}{n})^2+\taylorcos{2}{N}{(\frac{m}{n})})-\frac{|\frac{m}{n}|^{2N+1}2^{2N+1}}{(2N+1)!}>0
\end{align*}

From Proposition <ref> (see <ref>) we conclude that the quantity $\taylorcos{2}{N}{(\frac{m}{n})}$ is positive. Since $|\frac{m}{n}|<1$, it suffices to show:

\begin{align}
&\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil -2}}(1-\frac{1}{2}(\frac{m}{n})^2)-\frac{|\frac{m}{n}|^{2N+1}2^{2N+1}}{(2N+1)!}>\nonumber \\
&\frac{(2N-2)^2}{2^{p -2}}(1-\frac{1}{2}(\frac{m}{n})^2)-\frac{2^{2N+1}}{(2N+1)!}(\frac{m}{n})^2>0\label{eqn:cos_part_1_rewrite1}
\end{align}

Inequality (<ref>) is equivalent to:

\begin{align*}
\end{align*}

Since $(\frac{m}{n})^2<1$, it is sufficient to choose an $N$ that satisfies the following inequality:

\begin{align*}
\frac{1}{2}+\frac{2^{p+2N-1}}{(2N+1)!(2N-2)^2}<1
\end{align*}

Thus, we should choose an $N=2k+1$ that satisfies $\frac{(2N+1)!(2N-2)^2}{2^{2N}}>2^p$.

Proof for inequality (<ref>): To prove the inequality, we estimate the amount of precision that is lost when we approximate $\taylorcos{0}{N}{x}$ by $\taylorcos{0}{N}{(\frac{m}{n})}$. The sign of the terms $\frac{(-1)^i}{(2i)!}x^{2i}$ alternates between positive and negative. Thus, adding two arbitrary terms with different signs from the expansion can potentially cause significant loss of precision (see Theorem <ref>.<ref>). To avoid this, we first consider pairs of terms for which addition can be calculated with a bounded loss of precision. Afterwards, we calculate the summation of these pairs. The following identity shows our computation scheme:

\begin{align*}
\taylorcos{0}{N}{x}=\sum_{i=0}^{k}\frac{x^{4i}}{(4i)!}(1-\frac{x^2}{(4i+1)(4i+2)})
\end{align*}

Choosing $N=2k+1$ allows us to pair the terms of the summation.

Since $x$ is given with precision $p$, we can approximate $x^2$ with precision $p-2$ (see Theorem <ref>):

\begin{align*}
\end{align*}

Multiplying the approximation of $x^2$ by a constant does not influence the precision:

\begin{align*}
\end{align*}

We approximate $1-\frac{x^2}{(4i+1)(4i+2)}$ by $1-\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)}$. To estimate lose of precision in our approximation, we calculate the quantity $\log_2(\frac{1+\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)}}{1-\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)}})$ (see Theorem <ref>.<ref>):

\begin{align*}
\log_2(\frac{1+\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)}}{1-\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)}})\leq \log_2(\frac{1+\frac{1}{2}}{1-\frac{1}{2}})=\log_2 3<2
\end{align*}

Thus, we lose at most $2$ units of precision:

\begin{align*}
\end{align*}

Approximating $x^{4i}$ by $(\frac{m}{n})^{4i}$ reduces the precision by $2\lceil\log_2 (4i) \rceil$ units (see Lemma <ref> in <ref>); in the worst case we lose $2\lceil\log_2 (4k) \rceil$ units of precision:

\begin{align*}
|x^{4i}-(\frac{m}{n})^{4i}|<\frac{1}{2^{p-2\lceil\log_2 (4k) \rceil}}|(\frac{m}{n})^{4i}|=\frac{1}{2^{p-2\lceil\log_2 (2N-2) \rceil}}|(\frac{m}{n})^{4i}|
\end{align*}

Multiplying $(\frac{m}{n})^{4i}$ by a constant factor does not influence the precision of the calculation:

\begin{align*}
|\frac{x^{4i}}{(4i)!}-\frac{(\frac{m}{n})^{4i}}{(4i)!}|<\frac{1}{2^{p-2\lceil\log_2 (2N-2) \rceil}}|\frac{(\frac{m}{n})^{4i}}{(4i)!}|
\end{align*}

We approximate $\frac{x^{4i}}{(4i)!}(1-\frac{x^2}{(4i+1)(4i+2)})$ by multiplying the approximations $\frac{(\frac{m}{n})^{4i}}{(4i)!}$ and $(1-\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)})$ (see Theorem <ref>):

\begin{align*}
~~~~~= &\frac{1}{2^{p-2\lceil\log_2(2N-2)\rceil-2}}|\frac{(\frac{m}{n})^{4i}}{(4i)!}(1-\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)})|
\end{align*}

To obtain the last equality, we use the assumptions $\frac{(2N+1)!(2N-2)^2}{2^{2N}}>2^p$ and $N=2k+1$; we conclude that $N\geq 3$.

The terms $\frac{(\frac{m}{n})^{4i}}{(4i)!}(1-\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)})$ are positive for $0\leq i\leq k$ (see Proposition <ref> in <ref>) and hence we can approximate the summation $\sum_{i=0}^{k}\frac{x^{4i}}{(4i)!}(1-\frac{x^2}{(4i+1)(4i+2)})$ without losing precision (see Theorem <ref>.<ref>):

\begin{align*}
~~~~~&\frac{1}{2^{p-2\lceil\log_2(2N-2) \rceil -2}}|\sum_{i=0}^{k}\frac{(\frac{m}{n})^{4i}}{(4i)!}(1-\frac{(\frac{m}{n})^2}{(4i+1)(4i+2)})|
\end{align*}

Theorem <ref> estimates loss of precision in $\cos(x)$ in the base interval independently of the argument $x$. To show that iterative computations can be avoided in the base interval, we calculate $|\frac{xf'(x)}{f(x)}|$ for $f(x)=\cos(x)$:

\begin{align}
|\frac{xf'(x)}{f(x)}|=|\frac{-x\cdot \sin(x)}{\cos(x)}|=|x\cdot \tan(x)|\label{eqn:pert_cos}
\end{align}

From Proposition <ref> we conclude that $|x\cdot \tan(x)|<\tan(1)$ for $x\in (-1,1)$. Thus, we can approximate $\cos(x)$ in the base interval in one pass.

§.§ Extending Base Interval Approximations

In Section <ref> and <ref>, we discussed approximations for $\sin(x)$ and $\cos(x)$ in the base interval. In what follows, we show that range reduction identities can be used to extend these approximations to calculate sine and cosine for an argument $x$ represented by $(m,n,p)$ where $|\frac{m}{n}|\geq 1$.

In our calculations, we use a representation $(m',n',p)$ of $\pi$. This representation can be obtained based on our approximation for $\arctan(x)$ (see Section <ref>) and the following identity:

\begin{align*}
\pi=4\arctan(1)
\end{align*}

§.§.§ Sine

Let $x$ be a real number represented by $(m,n,p)$ such that $|\frac{m}{n}|\geq 1$ and $(m',n',p)$ be a representation for $\pi$. Suppose $\frac{\om}{\on}=\frac{m}{n}+k\frac{m'}{n'}$ where $k\in\integer$ and $0<\frac{\om}{\on}<\frac{m'}{n'}$. The value of $\sin(x)$ can be approximated as follows:

If $0<\frac{\om}{\on}<1$, then $\sin(x)$ can be represented by $(m_1,n_1,p-i_1)$ where:

\begin{align*}
\frac{m_1}{n_1}=&\taylorsin{0}{N_1}{(\frac{\om}{\on})}
\end{align*}

If $1\leq \frac{\om}{\on}<2$, then $\sin(x)$ can be represented by $(m_2,n_2,p-i_2)$ where:

\begin{align*}
\frac{m_2}{n_2}=&2\bigg(\taylorsin{0}{N_1}{(\frac{\om}{2\on})}\bigg)\cdot\bigg(\taylorcos{0}{N_2}{(\frac{\om}{2\on})}\bigg)
\end{align*}

If $2\leq \frac{\om}{\on}<\frac{m'}{n'}$, then $\sin(x)$ can be represented by $(m_3,n_3,p-i_3)$ where:

\begin{align*}
\frac{m_3}{n_3}=&8\bigg(\taylorsin{0}{N_1}{(\frac{\om}{4\on})}\bigg)\cdot\bigg(\taylorcos{0}{N_2}{(\frac{\om}{4\on})}\bigg)\\
& \cdot\bigg(\taylorsin{0}{N_3}{(\frac{m'}{4n'}-\frac{\om}{4\on})}\bigg)\cdot\bigg(\taylorcos{0}{N_4}{(\frac{m'}{4n'}-\frac{\om}{4\on})}\bigg)
\end{align*}

In the approximations above, $s_1,s_2,N_1,N_2,N_3,N_4\in \nat$ are the smallest natural numbers satisfying:

\begin{gather*}
s_1\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,k\frac{m'}{n'})}{\max(|\frac{m}{n}|,k\frac{m'}{n'})}}{1-\frac{\min(|\frac{m}{n}|,k\frac{m'}{n'})}{\max(|\frac{m}{n}|,k\frac{m'}{n'})}}) ~,~s_2\geq \log_2(\frac{1+\frac{\min(\frac{m'}{4n'},\frac{\om}{4\on})}{\max(\frac{m'}{4n'},\frac{\om}{4\on})}}{1-\frac{\min(\frac{m'}{4n'},\frac{\om}{4\on})}{\max(\frac{m'}{4n'},\frac{\om}{4\on})}})\\
(\frac{5}{6})(\frac{(2N_1+2)!(2N_1-1)^2}{2^{2N_1}})>2^{p-s_1} ,\frac{(2N_2+1)!(2N_2-2)^2}{2^{2N_2}}>2^{p-s_1}\\
\end{gather*}

and $t_1,t_2\in \nat$ are defined as follows:

\begin{align*}
t_1=&2\lceil\log_2(2N_1-1) \rceil+3\\
t_2=&\max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3)\\
t_3=&\max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3,\\
&\phantom{\max(} 2\lceil\log_2(2N_3-1) \rceil+3,2\lceil\log_2(2N_4-2) \rceil+3)
\end{align*}

Since $|\frac{m}{n}|\geq 1$, we can choose $k\in \integer$ such that $\frac{\om}{\on}=\frac{m}{n}+k\frac{m'}{n'}$ and $0< \frac{\om}{\on}<\frac{m'}{n'}$. Suppose $y=x+k\pi$. We use the following identity to calculate $\sin(x)$:

\begin{align}
\sin(x)=\sin(x+k\pi)=\sin(y)\label{eqn:red_sin}
\end{align}

We approximate $y$ by $\frac{\om}{\on}=\frac{m}{n}+k\frac{m'}{n'}$ (see Theorem <ref>.<ref>). By performing this approximation, we lose $s_1\in\nat$ units of precision where $s_1$ is the smallest number satisfying:

\begin{align*}
s_1\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,k\frac{m'}{n'})}{\max(|\frac{m}{n}|,k\frac{m'}{n'})}}{1-\frac{\min(|\frac{m}{n}|,k\frac{m'}{n'})}{\max(|\frac{m}{n}|,k\frac{m'}{n'})}})
\end{align*}

We consider three cases for calculating $\sin(y)$:

* Suppose $0<\frac{\om}{\on}<1$. We apply Theorem <ref>. Thus, $t_1=2\lceil\log_2(2N_1-1) \rceil+3$ units of precision is lost in the approximation of $\sin(y)$. A total of $s_1+t_1$ units of precision is lost in the approximation of $\sin(x)$.

* Suppose $1\leq \frac{\om}{\on}<2$. We use the following identity to bring the argument within the base interval:

\begin{align}
\sin(y)&=2\sin(\frac{y}{2})\cos(\frac{y}{2})\label{eqn:sin_id1}
\end{align}

We approximate $\frac{y}{2}$ by $\frac{\om}{2\on}$. Since $\frac{1}{2}\leq \frac{\om}{2\on}<1$, we apply Theorem <ref> and <ref> to approximate $\sin(\frac{y}{2})$ and $\cos(\frac{y}{2})$, respectively. Loss of precision in these calculations is as follows:

\begin{align*}
t_2=\max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3)
\end{align*}

The multiplication in equality (<ref>) reduces the precision by $2$ units (see Theorem <ref>). A total of $s_1+t_2+2$ units of precision is lost in the approximation of $\sin(x)$.

* Suppose $2\leq \frac{\om}{\on}<\frac{m'}{n'}$. We use the following identity to bring the argument within the base interval.

\begin{align}
\sin(y)&=2\sin(\frac{y}{2})\cos(\frac{y}{2})\nonumber\\
\end{align}

We approximate $\frac{y}{4}$ and $\frac{\pi}{4}-\frac{y}{4}$ by $\frac{\om}{4\on}$ and $\frac{m'}{4n'}-\frac{\om}{4\on}$, respectively. We lose $s_2\in \nat$ units of precision in this approximation (see Theorem <ref>.<ref>) where $s_2$ is the smallest number satisfying:

\begin{align*}
s_2\geq \log_2(\frac{1+\frac{\min(\frac{m'}{4n'},\frac{\om}{4\on})}{\max(\frac{m'}{4n'},\frac{\om}{4\on})}}{1-\frac{\min(\frac{m'}{4n'},\frac{\om}{4\on})}{\max(\frac{m'}{4n'},\frac{\om}{4\on})}})
\end{align*}

We apply Theorem <ref> and <ref> to approximate $\sin(\frac{y}{4}), \cos(\frac{y}{4}), \sin(\frac{\pi}{4}-\frac{y}{4})$, and $\cos(\frac{\pi}{4}-\frac{y}{4})$. Lose of precision in these approximations can be calculated as follows:

\begin{align*}
t_3=\max(&2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3,\\
&2\lceil\log_2(2N_3-1) \rceil+3,2\lceil\log_2(2N_4-2) \rceil+3)
\end{align*}

The three multiplications in $\sin(\frac{y}{4})\cos(\frac{y}{4})\sin(\frac{\pi}{4}-\frac{y}{4})\cos(\frac{\pi}{4}-\frac{y}{4})$ reduce the precision by $6$ units (see Theorem <ref>). A total of $s_1+s_2+t_3+6$ units of precision is lost in the approximation of $\sin(x)$.

[1]
$\expr$ has the shape $\sin{x}$
Choose an odd $N$ such that $(\frac{5}{6})(\frac{(2N+2)!(N-1)^2}{2^{2N}})>2^{p+2\lceil\log_2(2N-1) \rceil+3}$
$p_x\gets p+2\lceil\log_2(2N-1) \rceil+3$

$\frac{m}{n}\gets \Call{Compute}{x,p_x}$
$-1<\frac{m}{n}<1$ Theorem <ref>
$\frac{m_0}{n_0}\gets \taylorsin{0}{N}{(\frac{m}{n})}$
$\frac{m'}{n'}\gets \Call{Compute}{4\arctan(1),p_x}$
Choose $k\in\integer$ such that $0<\frac{m}{n}+k\frac{m'}{n'}<\frac{m'}{n'}$
$\frac{\om}{\on}\gets \frac{m}{n}+k\frac{m'}{n'}$
Choose $s_1\in\nat$ such that $s_1\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,k\frac{m'}{n'})}{\max(|\frac{m}{n}|,k\frac{m'}{n'})}}{1-\frac{\min(|\frac{m}{n}|,k\frac{m'}{n'})}{\max(|\frac{m}{n}|,k\frac{m'}{n'})}})$
Choose an odd $N_1$ such that $(\frac{5}{6})(\frac{(2N_1+2)!(2N_1-1)^2}{2^{2N_1}})>2^{p_x-s_1}$
Choose an odd $N_2$ such that $\frac{(2N_2+1)!(2N_2-2)^2}{2^{2N_2}}>2^{p_x-s_1}$
$0<\frac{\om}{\on}<1$ Theorem <ref>.<ref>
$t_1\gets 2\lceil\log_2(2N_1-1) \rceil+3$
$p_x-s_1-t_1\geq p$
$p_x\gets p_x+1$
$1\leq \frac{\om}{\on}<2$ Theorem <ref>.<ref>
$t_2\gets \max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3)$
$p_x-s_1-t_2-2\geq p$
$p_x\gets p_x+1$
Theorem <ref>.<ref>
Choose an odd $N_3$ such that  $(\frac{5}{6})(\frac{(2N_3+2)!(2N_3-1)^2}{2^{2N_3}})>2^{p_x-s_1-s_2}$
Choose an odd $N_4$ such that $\frac{(2N_4+1)!(2N_4-2)^2}{2^{2N_4}}>2^{p_x-s_1-s_2}$
[t]$t_3\gets \max(2\lceil\log_2(2N_1-1) \rceil+3, 2\lceil\log_2(2N_2-2) \rceil+3,$$\phantom{~~\max (}2\lceil\log_2(2N_3-1) \rceil+3, 2\lceil\log_2(2N_4-2) \rceil+3)$

Sine (Continued)
[1]
$p_x-s_1-s_2-t_3-6\geq p$
$\frac{m_3}{n_3}=8\big(\taylorsin{0}{N_1}{(\frac{\om}{4\on})}\big)\big(\taylorcos{0}{N_2}{(\frac{\om}{4\on})}\big)$$~~~~~~~~\phantom{\frac{m_3}{n_3}=8}\big(\taylorsin{0}{N_3}{(\frac{m'}{4n'}- \frac{\om}{4\on})}\big)\big(\taylorcos{0}{N_4}{(\frac{m'}{4n'}-\frac{\om}{4\on})}\big)$
$p_x\gets p_x+1$

Algorithm <ref> applies Theorem <ref> and <ref> to approximate $\sin(x)$ with arbitrary precision. In this algorithm, initially, we calculate $x$ with a precision that is adequate for calculations in the base interval. However, if the obtained approximation is outside the base interval, we use the half-angle formula or add rational multiples of $\pi$ to the argument (see equality (<ref>) and (<ref>)).

Observe that an arbitrary amount of precision can be lost in the approximation of $x+k\pi$ and $\frac{\pi}{4}-\frac{y}{4}$ when $x\approx -k\pi$. Hence, iterative computations might be necessary in the approximation (see Line <ref>,<ref>,<ref>,<ref> in Algorithm <ref>). To show that this is essential for sine, we reconsider the perturbation analysis in equality (<ref>) for $f(x)=\sin(x)$. The quantity $|x\cdot\cot(x)|$ can be arbitrary large for $x\approx -k\pi$. Thus, iterative computations are essential for approximating $\sin(x)$.

§.§.§ Cosine

Let $x$ be a real number represented by $(m,n,p)$ such that $|\frac{m}{n}|\geq 1$ and $(m',n',p)$ be a representation for $\pi$. Suppose $\frac{\om}{\on}=\frac{m}{n}+(2k+1)\frac{m'}{2n'}$ where $k\in\integer$ and $0<\frac{\om}{\on}<\frac{m'}{n'}$. The value of $\cos(x)$ can be approximated as follows:

If $0<\frac{\om}{\on}<1$, then $\cos(x)$ can be represented by $(m_1,n_1,p-i_1)$ where:

\begin{align*}
\frac{m_1}{n_1}=&(-1)^k\taylorsin{0}{N_1}{(\frac{\om}{\on})}
\end{align*}

If $1\leq \frac{\om}{\on}<2$, then $\cos(x)$ can be represented by $(m_2,n_2,p-i_2)$ where:

\begin{align*}
\frac{m_2}{n_2}=&2\cdot(-1)^k\bigg( \taylorsin{0}{N_1}{(\frac{\om}{2\on})}\bigg)\cdot\bigg(\taylorcos{0}{N_2}{(\frac{\om}{2\on})}\bigg)
\end{align*}

If $2\leq \frac{\om}{\on}<\frac{m'}{n'}$, then $\cos(x)$ can be represented by $(m_3,n_3,p-i_3)$ where:

\begin{align*}
\frac{m_3}{n_3}=&8\cdot (-1)^k \bigg(\taylorsin{0}{N_1}{(\frac{\om}{4\on})}\bigg)\cdot\bigg(\taylorcos{0}{N_2}{(\frac{\om}{4\on})}\bigg)\\
& \cdot\bigg(\taylorsin{0}{N_3}{(\frac{m'}{4n'}-\frac{\om}{4\on})}\bigg)\cdot\bigg(\taylorcos{0}{N_4}{(\frac{m'}{4n'}-\frac{\om}{4\on})}\bigg)
\end{align*}

In the approximations above, $s_1,s_2,N_1,N_2,N_3,N_4\in\nat$ are the smallest natural numbers satisfying:

\begin{align*}
s_1\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}{\max(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}}{1-\frac{\min(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}{\max(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}})~&,~s_2\geq \log_2(\frac{1+\frac{\min(\frac{m'}{4n'},\frac{\om}{4\on})}{\max(\frac{m'}{4n'},\frac{\om}{4\on})}}{1-\frac{\min(\frac{m'}{4n'},\frac{\om}{4\on})}{\max(\frac{m'}{4n'},\frac{\om}{4\on})}})\\
(\frac{5}{6})(\frac{(2N_1+2)!(2N_1-1)^2}{2^{2N_1}})>2^{p-s_1} ~&,~\frac{(2N_2+1)!(2N_2-2)^2}{2^{2N_2}}>2^{p-s_1}\\
\end{align*}

and $t_1,t_2\in \nat$ are defined as follows:

\begin{align*}
t_1=&2\lceil\log_2(2N_1-1) \rceil+3\\
t_2=&\max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3)\\
t_3=&\max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3,\\
&\phantom{\max(} 2\lceil\log_2(2N_3-1) \rceil+3,2\lceil\log_2(2N_4-2) \rceil+3)
\end{align*}

Since $|\frac{m}{n}|\geq 1$, we can choose $k\in \integer$ such that $\frac{\om}{\on}=\frac{m}{n}+(2k+1)\frac{m'}{2n'}$ and $0< \frac{\om}{\on}<\frac{m'}{n'}$. Suppose $y=x+(2k+1)\frac{\pi}{2}$. We use the following identity to calculate $\cos(x)$:

\begin{align*}
\cos(x)=(-1)^k\sin(x+(2k+1)\frac{\pi}{2})=(-1)^k\sin(y)
\end{align*}

We approximate $y$ by $\frac{\om}{\on}=\frac{m}{n}+(2k+1)\frac{m'}{2n'}$. From Theorem <ref>.<ref>, we can conclude that $s_1\in \nat$ units of precision will be lost in this approximation where $s_1$ is the smallest natural number satisfying:

\begin{align*}
s_1\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}{\max(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}}{1-\frac{\min(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}{\max(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}})
\end{align*}

We apply Theorem <ref> to approximate $\sin(y)$ and determine the loss of precision in the approximation. From Theorem <ref> we can conclude that the factor $(-1)^k$ does not influence the precision of the approximation.

[1]
$\expr$ has the shape $\cos{x}$
Choose an odd $N$ such that $\frac{(2N+1)!(2N-2)^2}{2^{2N}}>2^{p+2\lceil\log_2(2N-2) \rceil +2}$
$p_x\gets p+2\lceil\log_2(2N-2) \rceil+2$

$\frac{m}{n}\gets \Call{Compute}{x,p_x}$
$-1<\frac{m}{n}<1$ Theorem <ref>
$\frac{m'}{n'}\gets \Call{Compute}{4\arctan(1),p_x}$
Choose $k\in\integer$ such that $0<\frac{m}{n}+(2k+1)\frac{m'}{2n'}<\frac{m'}{n'}$
$\frac{\om}{\on}\gets \frac{m}{n}+(2k+1)\frac{m'}{2n'}$
Choose $s_1\in\nat$ such that $s_1\geq \log_2(\frac{1+\frac{\min(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}{\max(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}}{1-\frac{\min(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}{\max(|\frac{m}{n}|,(2k+1)\frac{m'}{2n'})}})$
Choose an odd $N_1$ such that $(\frac{5}{6})(\frac{(2N_1+2)!(2N_1-1)^2}{2^{2N_1}})>2^{p_x-s_1}$
Choose an odd $N_2$ such that $\frac{(2N_2+1)!(2N_2-2)^2}{2^{2N_2}}>2^{p_x-s_1}$
$0<\frac{\om}{\on}<1$ Theorem <ref>.<ref>
$t_1\gets 2\lceil\log_2(2N_1-1) \rceil +3$
$p_x-s_1-t_1\geq p$
$p_x\gets p_x +1$
$1\leq \frac{\om}{\on}<2$ Theorem <ref>.<ref>
$t_2\gets \max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3)$
$p_x-s_1-t_2-2\geq p$
$\frac{m_2}{n_2}=2\cdot (-1)^k\big(\taylorsin{0}{N_1}{(\frac{\om}{2\on})} \big)\big(\taylorcos{0}{N_2}{(\frac{\om}{2\on})}\big)$
$p_x\gets p_x +1$
Theorem <ref>.<ref>
Choose an odd $N_3$ such that $(\frac{5}{6})(\frac{(2N_3+2)!(2N_3-1)^2}{2^{2N_3}})>2^{p-s_1-s_2}$
Choose an odd $N_4$ such that $\frac{(2N_4+1)!(2N_4-2)^2}{2^{2N_4}}>2^{p-s_1-s_2}$
[t] $t_3\gets \max(2\lceil\log_2(2N_1-1) \rceil+3,2\lceil\log_2(2N_2-2) \rceil+3,$$\phantom{~~\max (} 2\lceil\log_2(2N_3-1) \rceil+3,2\lceil\log_2(2N_4-2) \rceil+3)$

Cosine (Continued)
[1]
$p_x-s_1-s_2-t_3\geq p$
[t]$\frac{m_3}{n_3}=8\cdot (-1)^k \big(\taylorsin{0}{N_1}{(\frac{\om}{4\on})}\big)\big(\taylorcos{0}{N_2}{(\frac{\om}{4\on})}\big)$

$p_x\gets p_x +1$

Algorithm <ref> applies Theorem <ref> and <ref> to approximate $\cos(x)$ with a desired precision. The algorithm first calculates $x$ with a precision that is sufficient for approximating $\cos(x)$ in the base interval. Similar to $\sin(x)$, we use range reduction identities to calculate $\cos(x)$ for an arbitrary $x$.

Observe that a significant amount of precision might be lost in the approximation when $x\approx \frac{-(2k+1)}{2}\pi$. Thus, iterative computations might be necessary in our approximation (see Line <ref>,<ref>,<ref>,<ref> in Algorithm <ref>). To show that these recomputations are essential, we reconsider the perturbation analysis of equality (<ref>). The quantity $|x\cdot \tan(x)|$ can be arbitrary large for $x\approx \frac{-(2k+1)}{2}\pi$ and hence iterative computations are unavoidable for $\cos(x)$.

§ RELATED WORK

An implementation for a bottom-up approach to exact real arithmetic is proposed in [16]. For a given expression, the inputs are computed with predefined precisions and a bottom-up scheme is used to determine the guaranteed precision of the output. Iterative computations are required if the obtained precision is not adequate. A formalization of a top-down approach in a theorem prover is proposed in [17]. The author first provides a definition for a metric space based on a ball relation. Afterwards, real numbers are defined as the completion of the metric space $\mathbb{Q}$. Rational operations are lifted to approximate operations on real numbers. This approach is optimized in [10].

Two closely-related top-down approaches based on absolute errors have been studied in [3, 15]. These approaches mainly differ in their approximations of the transcendental functions. In [3] the authors introduce a general way for calculating with Taylor expansions and apply this method to approximate the transcendental functions. In [15], the approximations of the transcendental functions are treated separately and in a more ad-hoc way. In contrast, our approach is based on relative errors. We provide detailed proofs of correctness for each operation and use perturbation analysis to identify essential recomputations.

In [1] the authors propose a layered framework for  computations with real numbers. The lowest layer is an implementation of floating point arithmetic. In the second layer, arithmetic operations are approximated using polynomial models. The highest layer supports more advanced features such as differential operations. Proof of correctness in Coq and an implementation based on [1] are also available. In this article, we have focused on approximating specific operations, whereas in [1] polynomial models are discussed in an abstract way without concrete examples from well-known arithmetic operations.

The approach introduced in [18] is based on linear fractional transformations (LFTs). Computations are encoded as trees of LFTs; various operations are defined to extract the result of a computation from the corresponding tree. However, this approach does not specify a top-down scheme to relate the desired precision in the output and the required operations on the expression tree. Expression trees are evaluated using lazy evaluation; computations terminate when adequate information is available in the root of a tree.

A symbolic approach to exact real arithmetic has been proposed in [8]. The author uses infinite binary sequences in the golden ratio base to represent real numbers. To calculate an expression, first the symbolic techniques available in Maple are applied to obtain a simplified expression. Additional Maple procedures are implemented by the author to extract binary sequences from simplified expressions. Performing operations on binary sequences is also possible. However, choosing a suitable balance between symbolic computations and direct manipulation of binary sequences depends on the given expression. As indicated in [8], using this approach to its full potentials requires expertise in Maple. Moreover, the procedure might need adaptations for each problem.

§ CONCLUSION

In this article, we proposed a simple representation for real numbers and discussed a top-down approach for approximating various arithmetic operations with arbitrary precision. The focus was on:

* providing complete algorithms and proofs of correctness for the approximations, and

* perturbation analysis to identify essential iterative computations.

Existing exact real arithmetic approaches have explored different representations for real numbers; approximations for algebraic operations and transcendental functions have also been proposed based on these representations. As far as we can see, proofs of correctness for existing approaches are restricted to basic operations. Moreover, no formal reasoning is provided to prove the necessity of iterative computations.

We envisage various extensions of the presented approach. From a practical point of view, some optimizations are essential. For example, the coefficients $m$ and $n$ in the representation $(m,n,p)$ can grow rapidly during computations. Thus, space efficiency is a relevant concern. One can consider an alternative representation in which large coefficients are represented in a more efficient way. Moreover, the computational efficiency of the transcendental functions can be improved by reducing the amount of required computations (i.e, number of rectangles in Riemann sums, number of terms in Taylor expansions) to guarantee the desired precision.

As discussed in Section <ref>, in certain computational problems, computing the expressions as they are would lead to loss of precision, whereas rewriting the expressions would allow us to compute them in one pass. Our top-down approach can be extended with a set of rewrite rules that transform problematic expressions into expressions that can be calculated in one pass.

This research was supported by the Dutch national program
COMMIT and carried out as part of the Allegio project.

[1]
P. Collins, M. Niqui, and N. Revol.
A validated real function calculus.
Mathematics in Computer Science, 5(4):437–467, 2011.

[2]
A. Feldstein and P. Turner.
Overflow, underflow, and severe loss of significance in
floating-point addition and subtraction.
IMA journal of numerical analysis, 6(2):241–251, 1986.

[3]
P. Gowland and D. Lester.
The correctness of an implementation of exact arithmetic.
In Proceedings of RNC 2000, volume 140, 2000.

[4]
P. Hertling.
Computable real functions: Type 1 computability versus type 2
In Proceedings of CCA 1996. Mathematik/Informatik,
Universität Trier, 1996.

[5]
N. Higham.
Accuracy and stability of numerical algorithms.
SIAM, 2002.

[6]
N. Higham.
Functions of matrices: theory and computation.
SIAM, 2008.

[7]
A. Hohmann and P. Deuflhard.
Numerical Analysis in Modern Scientific Computing: An
Introduction, volume 43.
Springer Science & Business Media, 2012.

[8]
T. Kelsey.
Exact numerical computation via symbolic computation.
In Proceedings of CCA 2000, pages 187–197. Springer, 2000.

[9]
D. Kincaid and E. Cheney.
Numerical analysis: mathematics of scientific computing,
volume 2.
American Mathematical Soc., 2002.

[10]
R. Krebbers and B. Spitters.
Type classes for efficient exact real arithmetic in Coq.
Logical Methods in Computer Science, 9(1:1):1–27, 2013.

[11]
G. Kreisel, D. Lacombe, and J. Shoenfield.
Partial recursive functionals and effective operations.
Constructivity in Mathematics, Studies in Logic and the
Foundations of Mathematics, pages 290–297, 2000.

[12]
B. Kushner and L. Lefman.
Lectures on constructive mathematical analysis, volume 60.
American Mathematical Soc., 1984.

[13]
T. M. library.
Visited: August 2015.

[14]
A. Markov.
On the continuity of constructive functions (in Russian).
Uspekhi Mat. Nauk (NS), 9:226–230, 1954.

[15]
V. Ménissier-Morain.
Arbitrary precision real arithmetic: design and algorithms.
The Journal of Logic and Algebraic Programming, 64(1):13–39,

[16]
N. Müller.
The iRRAM: Exact arithmetic in C++.
In Proceedings of CCA 2001, pages 222–252. Springer, 2001.

[17]
R. O'Connor.
Certified exact transcendental real number computation in Coq.
In Theorem Proving in Higher Order Logics, pages 246–261.
Springer, 2008.

[18]
P. Potts.
Exact real arithmetic using Mobius transformations.
PhD thesis, PhD-thesis, Imperial College London, 1998.

[19]
M. Richter and K. Wong.
Computable preference and utility.
Journal of Mathematical Economics, 32(3):339–354, 1999.

[20]
M. Spivak.
Benjamin, 1967.

[21]
K. Weihrauch.
Computable analysis: an introduction.
Springer Science & Business Media, 2012.

§ USEFUL PROPOSITIONS & LEMMAS

For any $p\in \nat$ the following inequalities hold:

\begin{align}
(1+\frac{1}{2^{p}})^2\leq 1+\frac{1}{2^{p-2}} \label{eqn:posprec2}\\
1-\frac{1}{2^{p-2}}\leq (1-\frac{1}{2^{p}})^2 \label{eqn:negprec2}
\end{align}

For inequality (<ref>) we can write:

\begin{align*}
(1+\frac{1}{2^{p}})^2=1+\frac{1}{2^{2p}}+\frac{1}{2^{p-1}}\leq 1+\frac{1}{2^{p-1}}+\frac{1}{2^{p-1}}=1+\frac{1}{2^{p-2}}
\end{align*}

Similarly for inequality (<ref>) we have:

\begin{align*}
(1-\frac{1}{2^{p}})^2=1+\frac{1}{2^{2p}}-\frac{1}{2^{p-1}}\geq 1-\frac{1}{2^{p-1}}>1-\frac{1}{2^{p-2}}
\end{align*}

For any $p\in \nat^{+}$ the following inequalities hold:

\begin{align}
\frac{2^{p}}{2^{p}-1}\leq 1+\frac{1}{2^{p-1}}\label{eqn:posprec1}\\
1-\frac{1}{2^{p-1}}\leq \frac{2^{p}}{2^{p}+1}\label{eqn:negprec1}
\end{align}

Since $p\in \nat^+$ we have $2^p\geq 2$. For inequality (<ref>) we have:

\begin{align*}
\frac{2^{p}}{2^{p}-1}=\frac{2^{p}-1+1}{2^{p}-1}=1+\frac{1}{2^{p}-1}\leq 1+\frac{1}{2^{p-1}}
\end{align*}

Similarly for inequality (<ref>) we can write:

\begin{align*}
\frac{2^{p}}{2^{p}+1}=\frac{2^{p}+1-1}{2^{p}+1}=1-\frac{1}{2^{p}+1}\geq 1-\frac{1}{2^{p-1}}
\end{align*}

For any $0<y<1$ the following inequality holds:

\begin{align}
\log_2(1+y)\leq \frac{y}{\ln(2)}
\end{align}

Based on Taylor's theorem, we can write the following expansion for $\log_2(1+y)$:

\begin{align*}
\log_2(1+y)=\sum_{i=1}^{N}(-1)^{i-1}\frac{y^i}{i\ln(2)}+\frac{(-1)^Ny^{N+1}}{(N+1)(1+c_y)^{N+1}\ln(2)}
\end{align*}

where $0<c_y<y$. For $N=1$ we obtain:

\begin{align*}
\log_2(1+y)=\frac{y}{\ln(2)}-\frac{y^2}{2(1+c_1)^2\ln(2)}
\end{align*}

The quantity $\frac{y^2}{2(1+c_1)^2\ln(2)}$ is positive and hence $\log_2(1+y)\leq \frac{y}{\ln(2)}$.

For $x\in \real\setminus\{0\}$ we have $|\frac{x}{(1+x^2)\arctan(x)}|<1$.

We analyze the derivative of $f(x)=|\frac{x}{(1+x^2)\arctan(x)}|=\frac{x}{(1+x^2)\arctan(x)}$ and find the intervals in which $f(x)$ is increasing/decreasing.

\begin{align*}
f'(x)=\frac{\arctan(x)(1-x^2)-x}{(1+x^2)^2(\arctan (x))^2}
\end{align*}

We analyze the derivative in the following cases:

* Suppose $0<x<1$. We use the Taylor expansion of $\arctan(x)$ to rewrite $\arctan(x)(1-x^2)-x$:

\begin{align*}
\arctan(x)(1-x^2)-x&=\taylorarctan{0}{\infty}{x}-x-x^2\arctan(x)\\
\end{align*}

The term $-x^2\arctan(x)$ is negative and $-x^{4i-1}(\frac{1}{4i-1}-\frac{x^2}{4i+1})$ is negative for $i\geq 1$ and $0<x<1$. Thus, the summation $\sum_{i=1}^{\infty}{-x^{4i-1}(\frac{1}{4i-1}-\frac{x^2}{4i+1})}-x^2\arctan(x)$ is negative and $f'(x)<0$. The function $f(x)$ is decreasing for $0<x<1$.

* Suppose $x\geq 1$. In this case, $\arctan(x)(1-x^2)-x<0$. Thus, $f'(x)<0$ and $f(x)$ is decreasing for $x\geq 1$.

* Suppose $x<0$. Since $f'(-x)=-f'(x)$, we can conclude from the first two cases that $f(x)$ is increasing when $x<0$.

The case analysis shows that $f(x)$ is bounded from above by $\lim_{x\to 0}f(x)=1$.

For any $y\in (-1,0) \cup (0,1)$ and $i\in \nat$ the following inequalities hold:

\begin{align*}
\frac{(-1)^{2i}}{(4i+1)!}y^{4i+1}+\frac{(-1)^{2i+1}}{(4i+3)!}y^{4i+3}>0 \qquad \text{if }y>0 \\
\frac{(-1)^{2i}}{(4i+1)!}y^{4i+1}+\frac{(-1)^{2i+1}}{(4i+3)!}y^{4i+3}<0 \qquad \text{if }y<0
\end{align*}

We can rewrite the summation as follows:

\begin{align*}
\frac{(-1)^{2i}}{(4i+1)!}y^{4i+1}+\frac{(-1)^{2i+1}}{(4i+3)!}y^{4i+3}=\frac{(-1)^{2i}}{(4i+1)!}y^{4i+1}(1-\frac{y^2}{(4i+2)(4i+3)})
\end{align*}

Since $1-\frac{y^2}{(4i+2)(4i+3)}\geq 1-\frac{1}{6}>0$ the sign of the summation is determined by the sign of $y^{4i+1}$.

For any $y\in (-1,0)\cup (0,1)$ and $i\in \nat$ the following inequality holds:

\begin{align*}
\frac{(-1)^{2i}}{(4i)!}y^{4i}+\frac{(-1)^{2i+1}}{(4i+2)!}y^{4i+2}>0
\end{align*}

We can rewrite the left hand side of the inequality as follows:

\begin{align*}
\frac{(-1)^{2i}}{(4i)!}y^{4i}+\frac{(-1)^{2i+1}}{(4i+2)!}y^{4i+2}=\frac{(-1)^{2i}}{(4i)!}y^{4i}(1-\frac{1}{(4i+1)(4i+2)}y^2)
\end{align*}

Since $1-\frac{1}{(4i+1)(4i+2)}y^2\geq 1-\frac{1}{2}>0$, the inequality holds.

For $x\in(-1,1)\setminus\{0\}$ we have $|x\cdot\cot(x)|<1$.

For the interval $x\in(-1,1)\setminus\{0\}$, we have $|x\cdot\cot(x)|=x\cdot\cot(x)$. We analyze the derivative of $f(x)=x\cdot\cot(x)$ in $(-1,1)\setminus\{0\}$.

\begin{align*}
\end{align*}

The point $x=0$ is a critical point for $f(x)$. Since $|\sin(2x)|\leq |2x|$, $f(x)$ is increasing in $(-1,0)$ and decreasing in $(0,1)$. Thus, $f(x)$ is bounded from above by $\lim_{x\to 0}f(x)=1$.

For $x\in(-1,1)$ we have $|x\cdot\tan(x)|<\tan(1)$.

For the interval $x\in(-1,1)$, we have $|x\cdot\tan(x)|=x\cdot\tan(x)$. To determine an upper-bound for $x\cdot\tan(x)$, we analyze the derivative of $f(x)=x\cdot\tan(x)$:

\begin{align*}
\end{align*}

The quantity $\sin(2x)+2x$ is negative in $(-1,0)$ and positive in $(0,1)$. Thus, $f(x)$ is decreasing in $(-1,0)$ and increasing in $(0,1)$. We conclude that $f(x)<\tan(1)$.

Let $x$ be a real number represented by $(m,n,p)$. Then $x^i$ can be represented by $(m^i,n^i,p-2\lceil \log_2^i \rceil)$.

We can apply Theorem <ref> to calculate $x^i=x^{\lceil \frac{i}{2} \rceil}\times x^{\lfloor \frac{i}{2} \rfloor}$. Let $k=\lceil \log_2^i \rceil$ and $P(i)$ denote the precision that we lose by calculating $x^i$. From Theorem <ref> we can write:

\begin{align*}
P(i)\leq P(\lceil \frac{i}{2} \rceil)+2\leq P(\lceil \frac{i}{2^k} \rceil)+2k=2k
\end{align*}

Thus, we lose $2k=2\lceil \log_2^i \rceil$ units of  precision by calculating $x^i$.