#  RLIBM-32: High Performance Correctly Rounded Math Libraries for 32-bit
Floating Point Representations
Rutgers Department of Computer Science Technical Report DCS-TR-754
Extended version of the PLDI 2021 paper

Jay P. Lim  0000-0002-7572-4017 Department of Computer ScienceRutgers
UniversityUnited States<EMAIL_ADDRESS>and  Santosh Nagarakatte
0000-0002-5048-8548 Department of Computer ScienceRutgers UniversityUnited
States<EMAIL_ADDRESS>
###### Abstract.

This paper proposes a set of techniques to develop correctly rounded math
libraries for 32-bit float and posit types. It enhances our RLibm approach
that frames the problem of generating correctly rounded libraries as a linear
programming problem in the context of 16-bit types to scale to 32-bit types.
Specifically, this paper proposes new algorithms to (1) generate polynomials
that produce correctly rounded outputs for all inputs using counterexample
guided polynomial generation, (2) generate efficient piecewise polynomials
with bit-pattern based domain splitting, and (3) deduce the amount of freedom
available to produce correct results when range reduction involves multiple
elementary functions. The resultant math library for the 32-bit float type is
faster than state-of-the-art math libraries while producing the correct output
for all inputs. We have also developed a set of correctly rounded elementary
functions for 32-bit posits.

correctly rounded math libraries, elementary functions, floating point,
posits, piecewise polynomials

††copyright: none††conference: Rutgers DCS Technical Report DCS-TR-754;
2021††journalyear: ;††ccs: Mathematics of computing Mathematical
software††ccs: Theory of computation Numeric approximation algorithms††ccs:
Mathematics of computing Linear programming

##  1\. Introduction

Math libraries provide implementations of elementary functions (_e.g._
$log(x)$, $e^{x}$, $cos(x)$) (Muller, 2005). They are crucial components in
various domains ranging from scientific computing to machine learning.
Designing math libraries is a challenging task because they are expected to
provide correct results for all inputs and also have high performance. These
elementary functions are typically approximated with some hardware supported
representation for performance.

Given a representation $\mathbb{T}$ with finite precision (_e.g._ , float),
the correctly rounded result of an elementary function $f$ for an input
$x\in\mathbb{T}$ is defined as the value of $f(x)$ computed with real numbers
and then rounded to a value in the representation $\mathbb{T}$. The IEEE-754
standard recommends the generation of correctly rounded results for elementary
functions. Seminal prior work on generating approximations for elementary
functions has resulted in numerous implementations that have reduced error
significantly (Lefèvre et al., 1998; Chevillard et al., 2010; Brisebarre et
al., 2006; Chevillard and Lauter, 2007; Chevillard et al., 2011; Kupriianova
and Lauter, 2014; Brunie et al., 2015; Jeannerod et al., 2011; Bui and Tahar,
1999; Ziv, 1991; Fousse et al., 2007). Further, numerous correctly rounded
libraries have also been developed (Daramy-Loirat et al., 2006; Fousse et al.,
2007). Unfortunately, they are not widely used due to performance
considerations. Moreover, widely used libraries do not produce correct results
for all inputs.

Mini-max approaches. Most prior approaches identify a polynomial that
minimizes the maximum error among all input points (_i.e._ , a mini-max
approach) compared to the real value of the elementary function using the
Weierstrass approximation theorem and the Chebyshev alternation theorem
(Trefethen, 2012). The Weierstrass approximation theorem states that if $f$ is
a continuous real-valued function on $[a,b]$ and $\epsilon>0$, there exists a
polynomial $P$ such that $|f(x)-P(x)|<\epsilon$ for all $x\in[a,b]$. The
Chebyshev alternation theorem provides the condition for such a polynomial: a
polynomial of degree $d$ that minimizes the maximum error will have at least
$d+2$ points where it has the absolute maximum error and the error alternates
in sign. Remez algorithm (Remes, 1934; Muller, 2005) is a procedure to
identify such mini-max polynomials. The maximum approximation error has to be
below the error threshold required to produce correct results for all inputs.

As approximating a polynomial in a small domain $[a,b]$ is much easier, the
input domain of the function is reduced using range reduction (Cody and Waite,
1980; Tang, 1990; Lim et al., 2020a). The approximated result is adjusted to
produce the result for the original input (_i.e._ , output compensation). Both
range reduction and polynomial evaluation in a representation with finite
precision will have some numerical errors. The combination of approximation
errors with the mini-max approach and numerical errors with polynomial
evaluation, range reduction, and output compensation can result in wrong
results.

Figure 1. Steps in our approach to generate correctly rounded libraries for
32-bit types ($\mathbb{T}$).

RLibm. Our RLibm approach (Lim et al., 2021, 2020a) generates polynomials that
approximate the correctly rounded result rather than the real value of the
elementary function. The generation of the polynomial considers errors in
polynomial approximation and numerical errors in polynomial evaluation, range
reduction, and output compensation to produce the correctly rounded output for
all inputs. The task of generating the polynomial is then structured as a
linear programming (LP) problem. The RLibm approach first computes the
correctly rounded result for each input in a target representation
$\mathbb{T}$ using an oracle (_e.g._ , the MPFR library (Fousse et al.,
2007)). Given the correctly rounded result for an input, it finds an interval
in double precision such that every value in the interval rounds to the
correctly rounded result, which is called the rounding interval. The rounding
intervals are further constrained to account for numerical errors during range
reduction and output compensation. Subsequently, it attempts to generate a
polynomial of degree $d$ using an LP solver, which when evaluated with an
input produces a result that lies within the rounding interval. Using the
RLibm approach, we have been successful in generating correctly rounded
libraries with 16-bit types such as bfloat16 and posit16.

Challenges in scaling to 32-bits. To extend our RLibm approach to 32-bit data
types, we have to address the following challenges. First, modern LP solvers
can handle a few thousand constraints. A naive use of the RLibm approach with
32-bit types will generate more than a billion constraints, which is beyond
the capabilities of current LP solvers. Second, it may not be feasible to
generate a single polynomial of a reasonable degree given the large number of
constraints. Third, LP solvers are sensitive to the condition number of the
system of constraints. LP solvers will not be able to solve an ill-conditioned
system of constraints. An effective range reduction is a strategy to address
it. Although there are excellent books on range reduction (Cody and Waite,
1980), these techniques need to be adapted to work with our RLibm approach.
Fourth, some range reduction strategies need multiple elementary functions
themselves (_e.g._ , $sinpi(x)$). Finally, we need to ensure that output
compensation does not experience pathological cancellation errors (_e.g._ ,
$cospi(x)$).

This paper. Our goal is to generate efficient implementations of elementary
functions that produce correctly rounded results for all inputs with 32-bit
types. This paper extends our RLibm approach to scale to 32-bit FP types to
address the challenges described above. We propose (1) sampling of inputs with
counterexample guided polynomial generation to handle the large input space,
(2) generation of piecewise polynomials for efficiency, (3) deduction of
rounding intervals when a range reduction technique uses multiple elementary
functions, and (4) modified range reduction techniques for some elementary
functions to address cancellation errors in output compensation. Figure 1
pictorially represents our approach to scale to 32-bit data types.

Counterexample guided polynomial generation. We sample inputs proportional to
the number of representable values in a given input domain $[a,b]$ with a
32-bit representation $\mathbb{T}$. To generate polynomials that produce the
correctly rounded result for every input, it is not necessary to consider
every input and its rounding interval. We primarily need to consider those
rounding intervals that are highly constrained. For each input in the sample,
we generate the oracle result using the MPFR library. We compute the rounding
interval in double precision (_i.e._ , set of values in the double type that
round to the oracle result). We generate LP constraints to create a polynomial
of degree $d$ such that it evaluates to a value in the rounding interval for
each input in the sample. If the initial sample generates a polynomial that
produces the correctly rounded output for all values in $[a,b]$, then the
process terminates. Otherwise, we add counterexamples to the sample and repeat
the process. The size of the sample is bounded by the number of constraints
that the LP solver can process.

Piecewise polynomials. When either the number of inputs in the sample exceeds
our LP constraint threshold or the LP solver is not able to generate a
polynomial, we split the input domain $[a,b]$ to $[a,b^{\prime})$ and
$[b^{\prime},b]$ to generate piecewise polynomials using the above process for
each input sub-domain. We choose the splitting point such that we can identify
the sub-domain quickly using a few bits of the input, which results in
efficient implementations. The ability to generate piecewise polynomials
ensures that our resultant polynomials are of a lower degree and provide
performance improvements when compared to state-of-the-art libraries.

Range reduction with multiple functions. We propose new algorithms to deduce
rounding intervals for a class of range reduction techniques that involve
multiple elementary functions. Range reduction reduces the input $x$ to
$x^{\prime}$. The creation of the polynomial happens with the reduced inputs.
The output of the polynomial $P(x^{\prime})$ should be adjusted to compute the
correctly rounded result for $x$, which is called output compensation. We have
to deduce the rounding intervals for the reduced input $x^{\prime}$ that
considers the numerical error in range reduction, polynomial evaluation, and
output compensation. We propose new techniques to create reduced rounding
intervals when range reduction uses multiple elementary functions (_e.g._ ,
$sinpi(x)$ in Section 2). These techniques allow us to perform range reduction
on functions that otherwise cause condition number issues with the LP
formulation (_i.e._ , $sinh(x)$ or $cosh(x)$). Further, we develop modified
range reduction techniques for some elementary functions to avoid cancellation
errors in output compensation (_e.g._ , $cospi(x)$ in Section 5).

The Rlibm-32 prototype. We have developed library generators and correctly
rounded libraries for multiple 32-bit data types: IEEE-754 float and posits.
Our elementary functions for floats are faster than existing libraries:
Intel’s libm, Glibc’s libm, CR-LIBM (Daramy-Loirat et al., 2006), and Metalibm
(Kupriianova and Lauter, 2014). Unlike existing libraries, our functions
produce correctly rounded results for all inputs. We have developed the first
correctly rounded implementations of functions for 32-bit posits.

##  2\. Overview of Our Approach with $sinpi(x)$

We provide an overview of our approach for generating piecewise polynomials
for $sinpi(x)$ (_i.e._ , $sin(\pi x)$) with a 32-bit float. The function
$sinpi(x)$ is defined for $x\in(-\infty,\infty)$. There are four billion
inputs with a 32-bit float. There are three kinds of special cases:

$sinpi(x)=\begin{cases}\pi x&\text{if }|x|<1.173\dots\times 10^{-7}\\\
0&\text{if }|x|>=2^{23}\\\ NaN&\text{if }x=NaN\text{ or
}x=\pm\infty\end{cases}$

For the first class of special cases, we compute $\pi x$ in double and round
the result to float, which produces the correctly rounded result for those
inputs.

###  2.1. Our Range Reduction for $sinpi(x)$

After considering special cases, there are close to 800 million float inputs
that need to be approximated with polynomials. Using RLibm’s approach directly
with an LP solver will fail. Next, we perform range reduction to reduce the
domain for polynomial approximation. The key idea is to use periodicity and
trigonometric identities of $sinpi(x)$. We transform input $x$ into
$x=2.0\times I+J$, where $I$ is an integer and $J\in[0,2)$. As a result of
periodicity, $sinpi(x)=sinpi(J)$. Next, we further split $J$ into $J=K+L$
where $K\in\\{0,1\\}$ is the integral part of $J$ and $L\in[0,1)$ is the
fractional part. Then, $sinpi(J)$ can be computed as,

$sinpi(J)=(-1)^{K}\times sinpi(L)$

Given that $sinpi$ between $[0.5,1)$ is a mirror image of values between
$[0,0.5)$, we further reduce as follows:

$L^{\prime}=\begin{cases}L&\text{if }L\leq 0.5\\\ 1.0-L&\text{if
}L>0.5\end{cases}$

From Sterbenz lemma (Sterbenz, 1974), the expression $1.0-L$ can be computed
exactly. Hence, $sinpi(L)=sinpi(L^{\prime})$. Even after reducing the input
$x$ to $L^{\prime}\in[0,0.5]$, there are around 184 million inputs with a
32-bit float in this reduced domain.

To enable easier polynomial approximation, we further reduce $L^{\prime}$ to a
value between $[0,\frac{1}{512}]$. We split $L^{\prime}$ as
$L^{\prime}=\frac{N}{512}+R$ where $N$ is an integer in the set
$\\{0,1,\dots,255\\}$ and $R$ is a fraction that lies in $[0,\frac{1}{512}]$.
There are $110$ million reduced inputs in $R$ ignoring special cases. Now,
$sinpi(L^{\prime})$ can be computed using the trigonometric identity
$sinpi(a+b)=sinpi(a)cospi(b)+cospi(a)sinpi(b)$ as follows,

$sinpi(L^{\prime})=sinpi(\frac{N}{512})cospi(R)+cospi(\frac{N}{512})sinpi(R)$

We precompute the values for $sinpi(\frac{N}{512})$ and $cospi(\frac{N}{512})$
in lookup tables (_i.e._ , 512 values in total). Finally, we approximate
$sinpi(R)$ and $cospi(R)$ for the reduced input domain
$R\in[0,\frac{1}{512}]$. _To approximate $sinpi(x)$ for the entire domain, the
range reduction requires us to approximate $sinpi$ and $cospi$ over the
reduced domain $R$_. We can compute the result for $sinpi(x)$ as follows,

$sinpi(x)=(-1)^{K}\times(sinpi(\frac{N}{512})cospi(R)+cospi(\frac{N}{512})sinpi(R))$

Figure 2. (a) A 32-bit float input $x_{1}=1.95\dots\times 10^{-3}$ and its
correctly rounded result of $sinpi(x_{1})$ (shown with a black circle). The
rounding interval in double is colored gray. (b) Input $x_{2}=2.14\dots\times
10^{-2}$, its correctly rounded result for $sinpi(x_{2})$, and the rounding
interval. (c) The reduced intervals for $sinpi(R)$ corresponding to $x_{1}$
and $x_{2}$ so that both $x_{1}$ and $x_{2}$ produce correctly rounded
results, respectively. Both $x_{1}$ and $x_{2}$ map to the same reduced input
$R$. The common interval for $R$ is highlighted with darker color. (d) To
approximate $sinpi(R)$, we create a piecewise polynomial with 32 sub-domains
in total. We use the 5-bits in the double representation of $R$ to identify
the sub-domain for the piecewise polynomial. (e) Our LP formulation for
generating a piecewise polynomial of degree 1 for the sub-domain with bit-
pattern ($10001$) with two reduced inputs in the sample. (f) The resulting
coefficients returned by the LP solver.

###  2.2. Generating Piecewise Polynomials for $sinpi(x)$

To produce correctly rounded results for $sinpi(x)$, our approach involves the
following steps. First, we identify the correctly rounded result and the
rounding interval for each input in the entire domain. Second, we identify the
reduced rounding interval after range reduction. Third, we split the reduced
domain into sub-domains to generate piecewise polynomials. Fourth, we perform
counterexample guided polynomial generation for each sub-domain. Finally, we
validate the generated piecewise polynomials for the entire input domain.

Step 1: Identifying the correctly rounded result and the rounding interval.
For each input $x$, we first identify the correctly rounded result of
$sinpi(x)$ using an oracle. Then, we identify an interval of values $[l,h]$ in
double where all values in the interval rounds to the correctly rounded
result. We call this interval the rounding interval. If our polynomial
approximation produces a value in the rounding interval, the rounded result is
the correct result. Consider the inputs:

$\displaystyle x_{1}$ $\displaystyle=1.95312686264514923095703125\times
10^{-3}$ $\displaystyle x_{2}$
$\displaystyle=2.148437686264514923095703125\times 10^{-2}$

We show the correctly rounded result of $sinpi(x)$ for these inputs with a
black circle in Figure 2(a) and Figure 2(b), respectively. It also shows the
rounding interval in gray.

Step 2: Identifying the reduced interval for input R. Range reduction
transforms input $x$ into $R$. To produce the result for $sinpi(x)$, we need
to compute both $sinpi(R)$ and $cospi(R)$ (_i.e._ , multiple elementary
functions). The result that we produce for $sinpi(R)$ and $cospi(R)$ should
allow us to produce the correct result for $sinpi(x)$ (_i.e._ , produce a
value within the rounding interval $[l,h]$ of input $x$).

To compute $sinpi(x)$, we will generate piecewise polynomials for $sinpi(R)$
and $cospi(R)$. Two inputs $x_{1}$ and $x_{2}$ (Figure 2(a) and 2(b)) map to
the same reduced input after range reduction,

$R=1.86264514923095703125\times 10^{-9}$

Now, we need to deduce an interval $[ls^{\prime},hs^{\prime}]$ for the output
of $sinpi(R)$ and an interval $[lc^{\prime},hc^{\prime}]$ for the output of
$cospi(R)$ such that the result of output compensation produces a value within
the rounding interval for $x$. We compute the correctly rounded value ($v$) of
$sinpi(R)$ in double using the oracle and set it as our initial guess for
$[ls^{\prime},hs^{\prime}]$ (_i.e._ , $[v,v]$ a singleton). Similarly, we
compute the interval $[lc^{\prime},hc^{\prime}]$ for $cospi(R)$. Now, we need
to check if these intervals are sufficient to produce the correct output for
the original input $x$. Section 3.2 provides our detailed algorithm. The key
idea is to simultaneously lower the lower bound for both $sinpi(R)$ and
$cospi(R)$ and check if output compensation produces the correct result for
all inputs. Similarly, we deduce the upper bound for both $sinpi(R)$ and
$cospi(R)$. The reduced interval for $sinpi(R)$ from the perspective of
$x_{1}$ is $[ls1^{\prime},hs1^{\prime}]$. Similarly, the reduced interval for
$sinpi(R)$ from the perspective of $x_{2}$ is $[ls2^{\prime},hs2^{\prime}]$.
These reduced intervals for $sinpi(R)$ corresponding to $x_{1}$ and $x_{2}$
are shown in Figure 2(c). They are not identical because our approach
considers the numerical error in both range reduction and output compensation.

Step 3: Splitting the reduced domain into sub-domains. Now that we have
reduced intervals for all reduced inputs, the next task is to generate
polynomials for $sinpi(R)$ and $cospi(R)$. We illustrate this process with
$sinpi(R)$. It is similar for $cospi(R)$. Even after range reduction, there
are approximately $110$ million unique reduced inputs for
$R\in[0,\frac{1}{512}]$. Using our counterexample guided polynomial generation
strategy (Step 4), we attempt to generate a polynomial for the entire reduced
domain. If we cannot generate a polynomial or the polynomial does not satisfy
the performance constraints, then we split the reduced input domain into
smaller sub-domains to generate piecewise polynomials. We iteratively split
the domain into smaller sub-domains until we can produce a polynomial that
produces the correct results for all inputs and satisfies the performance
criterion.

Let us say we want to generate 32 (_i.e._ $2^{5}$) piecewise polynomials for
the domain $[0,\frac{1}{512}]$. We use the bit-pattern of the reduced input
$R$ in double to identify the sub-domain. Although the domain of R is
$[0,\frac{1}{512}]$, the value of R in our reduced inputs ranges between
$[2^{-32},2^{-9}]$ along with $R=0$. There is a large gap of values between
the reduced input $0$ and $2^{-32}$. This is because we have already handled
special cases for the original input. Excluding the reduced input $R=0$, all
other reduced inputs in the double representation have the left-most six bits
identical. Hence, we use 5-bits after the six left-most bits to identify the
sub-domain for the piecewise polynomial. Figure 2(d) shows the reduced input
$R$, its double bit-pattern, and the 5-bits used to identify the sub-domain.

Step 4: Generating a polynomial for a sub-domain. The final step is to produce
a polynomial that approximates $sinpi(R)$ for a particular sub-domain. This
polynomial must produce a value within the reduced interval
$[ls^{\prime},hs^{\prime}]$ for each reduced input $R$ in the sub-domain. This
requirement can be encoded as a linear constraint for each reduced input,

$ls^{\prime}\leq P(R)\leq hs^{\prime}$

where $P(R)$ is a polynomial that approximates $sinpi(R)$.

We show the generation of the polynomial for sub-domain with bit-pattern
$10001$. First, we sample a portion of the reduced inputs (_e.g._ , 2 in
Figure 2(e)). Second, we encode the two reduced inputs and reduced intervals
as linear constraints to create a LP query (see Figure 2(e)). Third, we use a
LP solver to identify coefficients that satisfy the constraints. The generated
polynomial is shown in Figure 2(f). Fourth, we check if the generated
polynomial produces a value within the reduced interval for all inputs in the
sub-domain. In this case, there are two reduced inputs where the generated
polynomial does not produce a value within the reduced interval. Fifth, we add
both counterexamples (_i.e._ , reduced inputs) to the sample. Next, we create
a LP query using these four reduced inputs and intervals. Then, we check if
the generated polynomial satisfies all reduced inputs in the sub-domain
corresponding to bit-pattern $10001$.

After generating polynomials for all 32 sub-domains, we store the coefficients
of the piecewise polynomial in a table, which is indexed by the bit-pattern of
the reduced input that identifies the sub-domain. The approximation for the
elementary function $sinpi(x)$ is now ready. To produce the result for input
$x$, our library will perform range reduction on $x$, identify the reduced
input $R$, identify the sub-domain based on the bit-pattern of $R$, evaluate
the piecewise polynomial using the coefficients from the table, perform output
compensation, and round the result to a 32-bit float to produce the correctly
rounded result.

##  3\. Generating Piecewise Polynomials

Our goal is to generate polynomial approximations for elementary functions
$f(x)$ that produce the correctly rounded result for all inputs $x$ in 32-bit
target representations $\mathbb{T}$. Similar to our prior work on RLibm (Lim
et al., 2021, 2020a), we approximate the correctly rounded result rather than
the real value of $f(x)$. We extend it in three main directions. First, we
develop counterexample guided polynomial generation with sampling to make this
approach feasible with 32-bit types. Second, we design techniques to generate
piecewise polynomials, which provide performance improvements. Third, we
develop modified range reduction techniques for a class of elementary
functions and develop methods to deduce rounding intervals when range
reduction involves multiple functions.

Correctly rounded result. For a given input $x$ and elementary function $f$,
the output of our approximation is the correctly rounded result if it is equal
to the value of $f(x)$ computed with real numbers and then rounded to the
target representation. We use $RN_{\mathbb{T}}(f(x))$ to denote the rounding
function that rounds $f(x)$ computed with real numbers to target
representation $\mathbb{T}$. All internal computation such as range reduction,
polynomial evaluation, and output compensation is performed in representation
$\mathbb{H}$ where $\mathbb{H}$ has higher precision than $\mathbb{T}$. To
attain good performance, $\mathbb{H}$ is a representation that is supported in
hardware (_e.g._ , double).

Our approach. There are three main tasks in creating polynomial approximations
with our approach. First, we need to create a range reduction function, which
we denote as $RR_{\mathbb{H}}(x)$, that reduces input $x$ to a reduced input
$r$ in a smaller domain. Once we have the result of the elementary function
for the reduced input $r$ (let’s say $y^{\prime}=f(r)$), we need to develop an
output compensation function, which we denote as
$OC_{\mathbb{H}}(y^{\prime},x)$, to produce the result of $f(x)$ for input
$x$. Second, we need to generate polynomial approximations for each elementary
function $f_{i}(r)$ in the reduced domain (_e.g._ , there were two elementary
functions $sinpi$ and $cospi$ after range reduction in Section 2). We need to
generate polynomials $\Psi_{i}$ for each $f_{i}(r)$ in the reduced domain when
there are millions of reduced inputs in each reduced input domain. We have to
ensure that the polynomials generated for each $f_{i}(r)$ in the reduced
domain produce correctly rounded results for all inputs after output
compensation and polynomial evaluation is performed in $\mathbb{H}$. Third, we
may have to split the reduced input domain to generate piecewise polynomials
for each $f_{i}(r)$ to create efficient implementations.

1 Function  _CorrectPolys(_$f$ , $X$, $RR_{\mathbb{H}}$, $OC_{\mathbb{H}}$,
$d$_)_:

2 $Y\leftarrow\emptyset$

3 foreach  _$x\in X$_ do

4 $y\leftarrow RN_{\mathbb{T}}(f(x))$

5 $[l,h]\leftarrow$ RoundingInterval(_$y$ , $\mathbb{T}$, $\mathbb{H}$_)

6 $Y\leftarrow(x,[l,h])$

7

8 $\mathcal{L}\leftarrow$ ReducedIntervals(_$Y$ , $RR_{\mathbb{H}}$,
$OC_{\mathbb{H}}$_)

9 $Result\leftarrow\emptyset$

10

11 foreach  _$(f_{i},\mathcal{L}_{i})\in\mathcal{L}$_ do

12 if  _$\mathcal{L}_{i}\leftarrow\emptyset$_ then return _$\emptyset$_

13 $\Psi_{i}\leftarrow$ GenApproxFunc(_$\mathcal{L}_{i}$ , $d$_)

14 $Result\leftarrow(f_{i},\Psi_{i})\cup Result$

15

16 return _$Result$_

17

18

19Function  _RoundingInterval(_$y$ , $\mathbb{T}$, $\mathbb{H}$_)_:

20 $l\leftarrow min\\{v\in\mathbb{H}\mid v\leq y\text{ and
}RN_{\mathbb{T}}(v)=y\\}$

21 $h\leftarrow max\\{v\in\mathbb{H}\mid v\geq y\text{ and
}RN_{\mathbb{T}}(v)=y\\}$

22 return _$[l,h]$_

23

24

Algorithm 1 CorrectPolys computes piecewise polynomials of degree $d$ for each
elementary function $f_{i}$ used in output compensation, $OC_{\mathbb{H}}$, to
generate a math library for elementary function $f$. It produces the correctly
rounded result of $f(x)$ for each input $x\in X$. RoundingInterval computes
the rounding interval $[l,h]\subset\mathbb{H}$ of $y$ where all values in the
interval rounds to $y$ in $\mathbb{T}$. ReducedIntervals is shown in Algorithm
2 while GenApproxFunc is shown in Algorithm 3.

High-level sketch. Algorithm 1 provides a high-level sketch of our approach.
Given an elementary function $f(x)$ and a list of inputs $X$, we compute the
correctly rounded result $y$ in our target representation $\mathbb{T}$ (line
4) and compute the rounding interval of $y$ in $\mathbb{H}$ (lines 14-17). If
our approximation of $f(x)$ produces a value in the rounding interval, then
the result will round to $y$. Next, we compute the reduced input $r$ using
range reduction. The range reduction may require us to compute multiple
elementary functions $f_{i}$ to produce the result for $x$. Hence, we identify
the range of values that each function $f_{i}$ should produce such that the
result when used with output compensation produces a value in the rounding
interval of $y$ (line 7). We call this range of values for the reduced input
$r$ as the reduced interval (see Section 3.2). Finally, we approximate each
elementary function $f_{i}(r)$ used in output compensation with piecewise
polynomials of degree $d$ (line 11) with counterexample guided polynomial
generation and by using an LP solver. A single polynomial for each $f_{i}$ may
not be ideal for performance. To create efficient implementations, we
iteratively split the domain of the reduced input into multiple sub-domains
(see Section 3.3). Even such sub-domains for the reduced inputs can have
millions of reduced inputs. Hence, we create a sample of the reduced inputs,
generate constraints to ensure that the polynomial of degree $d$ produces a
value in the reduced interval for the reduced inputs in the sample, and query
the LP solver to solve for the coefficients. When the LP solver returns the
coefficients, we check whether the generated polynomial produces a value
within the reduced interval for all inputs in the sub-domain. We add any input
that violates the constraints to the sample and repeat this process. We call
this process as counterexample guided polynomial generation. At the end of
this process, our approach produces piecewise polynomials for each $f_{i}(r)$,
where the results of $f_{i}(r)$ when used with output compensation produces
the correctly rounded result for all inputs when rounded to $\mathbb{T}$.

###  3.1. Computing Rounding Intervals

Our approach approximates the correctly rounded result rather than the real
value. Hence, the first step is to identify the correctly rounded result using
an oracle and then identify all values in $\mathbb{H}$ that rounds to the
correct result in $\mathbb{T}$. As $\mathbb{H}$ has higher precision than
$\mathbb{T}$, there is a range of values in $\mathbb{H}$ that our approach can
produce and still round to the correctly rounded result in $\mathbb{T}$. We
call this range the rounding interval. Algorithm 1 illustrates our steps to
compute the rounding interval for each input $x\in X$ (lines 14-17). We
compute the oracle correctly rounded result, $y=f(x)$, using the MPFR math
library with a large number of precision bits. To compute the rounding
interval, we identify the smallest value $l\in\mathbb{H}$ that rounds to $y$
when rounded to $\mathbb{T}$ and the largest value $h\in\mathbb{H}$ that
rounds to $y$ when rounded to $\mathbb{T}$. This search procedure can be
efficiently implemented either using a binary search or by leveraging the
properties of $\mathbb{T}$ and $\mathbb{H}$. As long as our approach produces
a value in the rounding interval $[l,h]$ for input $x$, it will produce the
correctly rounded result.

1 Function  _ReducedIntervals(_$Y$ , $RR_{\mathbb{H}}$, $OC_{\mathbb{H}}$_)_:

2 if  _$OC_{\mathbb{H}}$ is not a monotonic function_ then return
_$\emptyset$_

3 $F\leftarrow\\{$list of functions used in $OC_{\mathbb{H}}\\}$

4 foreach  _$f_{i}\in F$_ do $\mathcal{L}_{i}\leftarrow\emptyset$

5 foreach  _$(x,[l,h])\in Y$_ do

6 $r\leftarrow RR_{\mathbb{H}}(x)$

7 $V\leftarrow\\{RN_{\mathbb{H}}(f_{i}(r))\mid f_{i}\in F\\}$

8 if  _$OC_{\mathbb{H}}(V,x)\notin[l,h]$_ then return _$\emptyset$_

9 //Set initial reduced range for each $f_{i}(r)$

10 $I^{\prime}\leftarrow\\{[v,v]\mid v\in V\\}$

11 //Decrease the lower bounds $l_{i}^{\prime}$ simulataneously

12 while  _true_ do

13 $A\leftarrow\\{$GetPrev(_$l_{i}^{\prime}$ , $\mathbb{H}$_)
$\mid[l_{i}^{\prime},h_{i}^{\prime}]\in I^{\prime}\\}$

14 if  _ $OC_{\mathbb{H}}(A,x)\notin[l,h]$ _ then break

15 $I^{\prime}\leftarrow\\{[$GetPrev(_$I_{i}^{\prime}$ ,
$\mathbb{H}$_)$,h_{i}^{\prime}]\mid[l_{i}^{\prime},h_{i}^{\prime}]\in
I^{\prime}\\}$

16 //Increase the upper bounds $h_{i}^{\prime}$ simulataneously

17 while  _true_ do

18 $B\leftarrow\\{$ GetNext(_$h_{i}^{\prime}$ , $\mathbb{H}$_)
$\mid[l_{i}^{\prime},h_{i}^{\prime}]\in I^{\prime}\\}$

19 if  _ $OC_{\mathbb{H}}(B,x)\notin[l,h]$ _ then break

20 $I^{\prime}\leftarrow\\{[I_{i}^{\prime},$ GetNext(_$h_{i}^{\prime}$ ,
$\mathbb{H}$_)$]\mid[l_{i}^{\prime},h_{i}^{\prime}]\in I^{\prime}\\}$

21 foreach  _$[l_{i}^{\prime},h_{i}^{\prime}]\in I^{\prime}$_ do

22
$\mathcal{L}_{i}\leftarrow\mathcal{L}_{i}\cup(r,[l_{i}^{\prime},h_{i}^{\prime}])$

23

24 return _$\\{(f_{i},\mathcal{L}_{i})\mid f_{i}\in F\\}$_

25

Algorithm 2 ReducedIntervals computes the reduced interval
$[l_{i}^{\prime},h_{i}^{\prime}]$ and the reduced input $r$ corresponding to
input $x$ for each function $f_{i}$ used with output compensation. If our
polynomial approximation for $f_{i}$ produces a value in
$[l_{i}^{\prime},h_{i}^{\prime}]$, then we can generate the correctly rounded
result for $x$. ReducedIntervals returns a list with
$(r,[l_{i}^{\prime},h_{i}^{\prime}])$ for each $f_{i}$. GetPrev(p,
$\mathbb{H}$) returns the preceding value of $p$ in $\mathbb{H}$ and
GetNext(p, $\mathbb{H}$) returns the succeeding value of $p$ in $\mathbb{H}$.

###  3.2. Computing Reduced Rounding Intervals

Range reduction is crucial for any technique that generates approximations for
elementary functions. It is particularly important with our approach for
32-bit types because the condition number of the LP problem increases
drastically if the input domain has both extremely large and small values.
Further, large inputs can cause overflows during the evaluation of a
polynomial with a large degree in $\mathbb{H}$.

After computing rounding intervals from Algorithm 1, we have a list of
constraints $(x,[l,h])$ that our approximation for $f(x)$ needs to satisfy for
each input $x$ to produce the correctly rounded result. The range reduction
and subsequent output compensation can require us to approximate multiple
elementary functions $f_{i}$. The next step is to identify reduced inputs to
$f_{i}$ and a range of values that $f_{i}$ should produce such that the result
of the output compensation produces a value in $[l,h]$ for each $x$. The input
to $f_{i}$ is the reduced input and the range of values that $f_{i}$ should
produce is the reduced interval.

Algorithm 2 shows the steps in deducing the reduced interval. For each
constraint $(x,[l,h])$, we can identify the reduced input $r$ by performing
range reduction on $x$ (line 6). However, computing the reduced interval is
challenging. We present an algorithm to deduce reduced intervals when output
compensation ($OC_{\mathbb{H}}$) is monotonic (either increasing or
decreasing), which is the case with all range reductions that we explore in
the paper.

To compute the reduced interval, we identify all functions $f_{i}$ used in
$OC_{\mathbb{H}}$ (line 4). Then, we compute the correctly rounded result
$v_{i}$ for each $f_{i}(r)$ in $\mathbb{H}$ using an oracle (line 7). If the
result of output compensation using $v_{i}$’s does not produce a value in the
rounding interval for $x$, then either the range reduction technique should be
redesigned or the precision of $\mathbb{H}$ should be increased.

Now, we have a candidate value (_i.e._ , $v_{i}$) for each $f_{i}(r)$ to
produce the correctly rounded result of $x$. We have to deduce the maximum
amount of freedom available for each $f_{i}(r)$. We initially set the reduced
intervals $[l_{i},h_{i}]$ for each $f_{i}$ to be $[v_{i},v_{i}]$ (line 10).
Next, we identify if we can decrease the lower bound of the intervals of
$f_{i}(r)$. For a given reduced input $r$ of input $x$, we check if using the
preceding values of $l_{i}$ in $\mathbb{H}$ for all $f_{i}$’s with output
compensation produces a value in the rounding interval $[l,h]$ of $x$. If it
does, then we widen the reduced interval by replacing each $l_{i}$ with the
preceding value. We repeat the process until the result of output compensation
using the preceding values no longer produces a value in $[l,h]$ (lines
12-15). This procedure to compute the lower bound can be efficiently
implemented by performing binary search between $v_{i}$ and the minimum
representable value.

Similarly, we identify if we can increase the upper bound of the interval for
each $f_{i}(r)$. For each upper bound $h_{i}$ of $f_{i}(r)$, we identify the
value that succeeds $h_{i}$ and check whether the result of output
compensation using the succeeding value produces a value in $[l,h]$. If it
does, then we widen the reduced interval by replacing each $h_{i}$ with the
succeeding value. We repeat the process until output compensation produces a
value outside the interval $[l,h]$ of input $x$ (lines 17-20). The upper bound
of the reduced interval can be efficiently computed by performing binary
search between $v_{i}$ and the maximum representable value. Finally, we store
the reduced constraints $(r,[l_{i}^{\prime},h_{i}^{\prime}])$ for each
function $f_{i}$ in a list $\mathcal{L}_{i}$.

Each $\mathcal{L}_{i}$ corresponding to $f_{i}$ contains reduced intervals
$(r,[l_{i}^{\prime},h_{i}^{\prime}])$ for the reduced input $r$ to produce a
correct result for input $x$. As multiple inputs can map to the same reduced
input $r$, there can be multiple reduced constraints
$(r,[l1_{i}^{\prime},h1_{i}^{\prime}])$ and
$(r,[l2_{i}^{\prime},h2_{i}^{\prime}])$ for the same reduced input $r$
corresponding to original inputs $x1$ and $x2$. The reduced intervals
$[l1_{i}^{\prime},h1_{i}^{\prime}]$ and $[l2_{i}^{\prime},h2_{i}^{\prime}]$
are not exactly identical to account for numerical errors in range reduction
and output compensation. Our polynomial approximation for $f_{i}$ must satisfy
the constraints $(r,[l1_{i}^{\prime},h1_{i}^{\prime}])$ to produce the
correctly rounded result for $x1$ and $(r,[l1_{i}^{\prime},h1_{i}^{\prime}])$
to produce the correctly rounded result for $x2$. Thus, we generate a single
combined interval by computing the common interval between them. If there is
no common interval between all reduced intervals corresponding to the same
reduced input, then it implies that there is no polynomial approximation for
$f_{i}$ that produces the correctly rounded results for all inputs $x$ in the
original domain. The library designer will have to redesign range reduction in
such cases.

1 Function  _GenApproxFunc(_$\mathcal{L}$ , $d$_)_:

2 $\mathcal{L}^{-}\leftarrow\\{(r,[l^{\prime},h^{\prime}])\in\mathcal{L}\mid
r<0\\}$

3 $\mathcal{L}^{+}\leftarrow\\{(r,[l^{\prime},h^{\prime}])\in\mathcal{L}\mid
r\geq 0\\}$

4 $\Psi^{-}\leftarrow$ GenApproxHelper(_$\mathcal{L}^{-},d$_)

5 $\Psi^{+}\leftarrow$ GenApproxHelper(_$\mathcal{L}^{+},d$_)

6 return _$\\{\Psi^{-},\Psi^{+}\\}$_

7

8Function  _GenApproxHelper(_$\mathcal{L}$ , $d$_)_:

9 $n\leftarrow 0$

10 while  _true_ do

11 $\Delta=$ SplitDomain(_$\mathcal{L}$ , $n$_)

12 $(status,\Psi)=$ GenPiecewise(_$\Delta$ , $d$_)

13 if  _$status=true$_ then return _$\Psi$_

14 $n\leftarrow n+1$

15

16

17Function  _GenPiecewise(_$\Delta$ , $d$_)_:

18 $\Psi\leftarrow\emptyset$

19 foreach  _$\Delta_{j}\in\Delta$_ do

20 $(status,\Psi_{j})\leftarrow$ GenPolynomial(_$\Delta_{j}$ , $d$_)

21 if  _status = false_ then return _$(false,\emptyset)$_

22 $\Psi\leftarrow\Psi\cup\Psi_{j}$

23

24 return _(true, $\Psi$)_

25

Algorithm 3 GenApproxFunc generates piecewise polynomials that produce a value
in the reduced interval for all reduced inputs in $\mathcal{L}$. It initially
attempts to produce a single polynomial for the entire reduced input domain.
If unsuccessful, then it splits the domain into multiple sub-domains.
SplitDomain (not defined in this algorithm) splits the reduced input domain
into sub-domains based on the bit-pattern of the reduced inputs in
$\mathbb{H}$. SplitDomain returns $\Delta$, which includes a set of reduced
constraints for each sub-domain $\Delta_{j}$. GenPiecewise generates a
polynomial for each sub-domain, which is shown in Algorithm 4.

###  3.3. Efficient Piecewise Polynomials

After the above steps, we have a list of reduced constraints
$(r,[l_{i}^{\prime},h_{i}^{\prime}])$ in $\mathcal{L}$ for each reduced input
$r$ and for each function $f_{i}$ that we need to approximate. The next step
in our approach is to generate polynomials that approximate $f_{i}$ and
satisfy the constraints in $\mathcal{L}_{i}$. Even after range reduction,
there can be hundreds of millions of reduced inputs. The counterexample guided
polynomial generation algorithm, which we describe in Section 3.4, can likely
generate a single polynomial in many cases. However, it will also have a large
degree and may not be efficient. To generate high performance math libraries,
we propose the generation of piecewise polynomials. Effectively splitting the
domain into smaller domains for the generation of piecewise polynomials is
essential to improve performance. Hence, we group the reduced input into sub-
domains based on the bit-patterns of the reduced input in $\mathbb{H}$.

Algorithm 3 describes our steps to generate piecewise polynomials. Range
reduction techniques for many elementary functions can create both positive
and negative reduced inputs (_e.g._ , $e^{x}$, $2^{x}$, $10^{x}$). The bit-
patterns for positive and negative reduced inputs in $\mathbb{H}$ will not
have common bits at the beginning (_e.g._ , the explicit sign bit
distinguishes positive and negative values in double). Hence, we separate the
reduced inputs (and their intervals) into two groups: $\mathcal{L}^{-}$ that
contains negative reduced inputs and $\mathcal{L}^{+}$ that contains non-
negative reduced inputs (lines 2-3). We create polynomial approximations for
each $\mathcal{L}^{-}$ and $\mathcal{L}^{+}$ (lines 4-5). This step also
allows us to subsequently group the reduced input into sub-domains in an
efficient manner.

If $\mathcal{L}$ contains only negative or positive reduced inputs, we try to
generate a single polynomial of degree $d$ that satisfies all reduced
constraints in $\mathcal{L}$ (line 11 and 17) using our counterexample guided
polynomial generation (see Section 3.4). If it cannot generate a polynomial of
degree $d$ that satisfies all constraints, then we split the reduced input
domain in $\mathcal{L}$ into multiple sub-domains (lines 9-13 in
GenApproxHelper). We iteratively split the domain of reduced inputs into
$2^{n}$ sub-domains based on the bit-pattern of $r$ in $\mathbb{H}$ (_i.e._ ,
SplitDomain call in line 10). To split the reduced input domain, we first
identify the smallest reduced input $R_{min}$ and the largest reduced input
$R_{max}$. Then, we compute the number of consecutive bits that are identical
in the bit-string representation of $R_{min}$ and $R_{max}$ in $\mathbb{H}$
starting from the most significant bit. We use the next $n$ bits to identify
the sub-domain for the piecewise polynomial. Subsequently, we group the
reduced inputs and reduced intervals based on the bit-pattern of the reduced
input into sub-domains ($\Delta$ returned by SplitDomain). We try to generate
a polynomial of degree $d$ that satisfies all reduced constraints in
$\Delta_{j}$ for all $\Delta_{j}$’s belonging to $f_{i}$ (lines 16-19). Using
bit-patterns of the reduced input in $\mathbb{H}$ allows us to efficiently
identify the sub-domain for the piecewise polynomial with two bitwise
operations (and and a shift). Once we generate a polynomial for each sub-
domain of every $f_{i}$, the coefficients of the polynomial are stored in a
table, which is indexed using the bit-pattern of the reduced input for each
$f_{i}$.

1 Function  _GenPolynomial(_$\Delta_{j}$ , $d$_)_:

2 $\mathcal{S}\leftarrow$ Sample(_$\Delta_{j}$_)

3 while  _true_ do

4 $\Psi_{j}\leftarrow$ GetCoeffsUsingLP(_$\mathcal{S}$ , $d$_)

5 if  _$\Psi_{j}=\emptyset$_ then return _$(false,\emptyset)$_

6 $(Done,\mathcal{S})\leftarrow$ Check(_$\Psi_{j}$ , $\Delta_{j}$,
$\mathcal{S}$_)

7 if  _$Done=true$_ then return _$(true,\Psi_{j})$_

8 if  _$|\mathcal{S}| >threshold$_ then return _$(false,\emptyset)$_

9

10

11

12Function  _Check(_$\Psi_{j}$ , $\Delta_{j}$, $\mathcal{S}$_)_:

13 $Done\leftarrow true$

14 foreach  _$(r,[l^{\prime},h^{\prime}])\in\Delta_{j}$_ do

15 if  _not $l^{\prime}\leq\Psi_{j}(r)\leq h^{\prime}$_ then

16 $\mathcal{S}\leftarrow\\{(r,[l^{\prime},h^{\prime}])\\}\cup\mathcal{S}$

17 $Done\leftarrow false$

18

19 return _$(Done,\mathcal{S})$_

Algorithm 4 GenPolynomial attempts to find a polynomial of degree $d$ that
satisfies the reduced input and interval constraints in $\Delta_{j}$ using our
counterexample guided sampling approach. If it is infeasible to find a
polynomial of degree $d$ or the size of the sample exceeds a threshold, then
it returns (false, $\emptyset$). GetCoeffsUsingLP generates the coefficients
of a polynomial that satisfies all constraints in $\mathcal{S}$ using the LP
solver. Check validates that the polynomial generated using the sample
satisfies all reduced input and interval constraints. We add counterexamples
(_i.e._ , all inputs where the polynomial does not satisfy the constraints) to
the sample and repeat the process.

###  3.4. Counterexample Driven Polynomial Generation

Once we have the reduced input and the reduced intervals, we structure the
problem of generating polynomials as a linear programming problem similar to
our prior work on RLibm (Lim et al., 2021, 2020a). Even after range reduction
and creation of sub-domains for the generation of piecewise polynomials, we
need to generate a polynomial approximation when there are several million
reduced inputs and reduced intervals in the context of 32-bit types. However,
they are beyond the capabilities of modern LP solvers, which can handle a few
thousand constraints. To address this issue, we propose counterexample guided
polynomial generation with sampling. The key insight is that we do not need to
add every reduced input and interval as a constraint in the LP formulation as
long as we identify and add the highly constrained intervals.

Our counterexample guided polynomial generation strategy takes as input the
set of reduced constraints $(r,[l^{\prime},h^{\prime}])$ corresponding to
reduced inputs that belong to a particular sub-domain. The goal is to generate
a polynomial of degree $d$ that produces a value in the reduced interval
$[l^{\prime},h^{\prime}]$ for each reduced input $r$. Each reduced input $r$
and the corresponding interval $[l^{\prime},h^{\prime}]$ specifies the
following linear constraint for a polynomial of degree $d$ that we want to
generate:

$l^{\prime}\leq c_{0}+c_{1}r+c_{2}r^{2}+\dots+c_{d}r^{d}\leq h^{\prime}$

The task of the polynomial generator is to find coefficients for the
polynomial.

To scale to 32-bit types, we sample a small fraction of the reduced input and
intervals. Algorithm 4 reports our counterexample guided polynomial generation
process. It takes two inputs: the degree of the polynomial and the set of
reduced inputs and intervals (_i.e._ , $\Delta_{j}$) for generating a
polynomial approximation for an elementary function $f_{i}$ on reduced inputs
for sub-domain $j$. We maintain the reduced inputs and their intervals in
increasing order. Then we uniformly sample the reduced inputs based on the
distribution of reduced inputs. If there are a large number of reduced inputs
in a particular region of the sub-domain, then our method has more samples
from that region. We also add highly constrained reduced inputs and intervals
(_i.e._ , the correctly rounded result and the lower bound/upper bound is less
than $\epsilon$, which is set by the math library designer) to the sample.

Then we express all constraints in the sample $(r,[l^{\prime},h^{\prime}])$
using a single system of linear inequalities and solve for the coefficients
using an LP solver (line 4). If there are $n$ points in the sample, the system
of linear inequalities is of the following form:

$\begin{bmatrix}l_{1}^{\prime}\\\ l_{2}^{\prime}\\\ \vdots\\\
l_{n}^{\prime}\end{bmatrix}\leq\begin{bmatrix}1&r_{1}&\dots&r_{1}^{d}\\\
1&r_{2}&\dots&r_{2}^{d}\\\ \vdots&\vdots&\vdots&\vdots\\\
1&r_{n}&\dots&r_{n}^{d}\\\ \end{bmatrix}\begin{bmatrix}c_{0}\\\ c_{1}\\\
\vdots\\\ c_{d}\end{bmatrix}\leq\begin{bmatrix}h_{1}^{\prime}\\\
h_{2}^{\prime}\\\ \vdots\\\ h_{n}^{\prime}\end{bmatrix}$

There are two issues with the polynomial generated using the sampled reduced
inputs that we need to address. First, as the LP solver returns coefficients
as real numbers, the coefficients are rounded to a value in $\mathbb{H}$. As a
result of rounding error, the result of polynomial evaluation for a particular
reduced input in the sample may not lie within its rounding interval. Second,
the polynomial generated using the sample may not satisfy the constraints for
the entire set of reduced inputs and their corresponding intervals.

We address the real coefficients issue with a search-and-refine procedure
similar to RLibm. When the LP solver returns real coefficients and we round it
to $\mathbb{H}$, we check whether evaluating the polynomial satisfies
constraints for every input in the sample. If it does not, then we select the
input and reduce its rounding interval (either replace the lower bound with
its succeeding value or replace the upper bound with the preceding value).
Then we repeat the above process until it generates a polynomial that either
satisfies all constraints in the sample when evaluated in $\mathbb{H}$ or
cannot find a polynomial of degree $d$. If we cannot find a polynomial that
satisfies all constraints in the sample, then we split the entire reduced
domain in $\mathcal{L}_{i}$ into even smaller sub-domains and repeat this
process.

If we successfully generate a polynomial $\Psi_{j}$ that satisfies all
constraints in the sample, then we check whether this polynomial satisfies all
constraints in $\Delta_{j}$ (line 10-15). If $\Psi_{j}$ satisfies all
constraints, then we return the polynomial (line 7). If there is any
constraint not satisfied by $\Psi_{j}$ in the entire set of reduced inputs,
then we add that reduced input and its interval to the sample (_i.e._ , adding
the counterexample in lines 12-13). We repeat the process of generating the
polynomial with the new sample. If the number of constraints in the sample
exceeds a threshold at any point, then we determine that we cannot generate a
polynomial for the sub-domain $\Delta_{j}$. Our function to generate the
coefficients for the polynomial (_i.e._ , GetCoeffsUsingLP) using an LP solver
generates a polynomial of a lower degree (than input degree $d$) if it is
possible to do so.

##  4\. Experimental Evaluation

Table 1. Generation of correctly rounded results for 32-bit floats with RLibm-32, Intel’s libm (float and double), glibc’s libm (float and double), CR-LIBM, and MetaLibm (float and double). ✓indicates that the library produces the correctly rounded result for all inputs. Otherwise, we use ✗. For each ✗, we show the number of inputs with wrong results. N/A indicates that the implementation is not available. | float
---
functions
| Using
---
RLibm-32
| Using
---
glibc float
| Using
---
glibc double
| Using
---
Intel float
| Using
---
Intel double
| Using
---
CR-LIBM
| Using
---
MetaLibm float
| Using
---
MetaLibm double
$\mathbf{ln(x)}$ | ✓ | ✗(4.2E5) | ✗(5) | ✗(1060) | ✗(5) | ✗(5) | N/A | N/A
$\mathbf{log2(x)}$ | ✓ | ✗(3.1E5) | ✓ | ✗(276) | ✓ | ✓ | N/A | N/A
$\mathbf{log10(x)}$ | ✓ | ✗(3.0E7) | ✗(1) | ✗(1.5E5) | ✗(1) | ✗(1) | N/A | N/A
$\mathbf{exp(x)}$ | ✓ | ✗(1.7E5) | ✓ | ✗(2.5E5) | ✓ | ✓ | ✗(5.1E8) | ✗(5.1E8)
$\mathbf{exp2(x)}$ | ✓ | ✗(1.7E5) | ✗(2) | ✗(7.2E5) | ✗(2) | N/A | ✗(6.5E7) | ✗(1026)
$\mathbf{exp10(x)}$ | ✓ | ✗(1.7E5) | ✓ | ✗(3.9E5) | ✓ | N/A | N/A | N/A
$\mathbf{sinh(x)}$ | ✓ | ✗(7.1E7) | ✗(2) | ✗(2.5E5) | ✗(2) | ✗(2) | N/A | N/A
$\mathbf{cosh(x)}$ | ✓ | ✗(1.8E7) | ✓ | ✗(1.4E5) | ✓ | ✓ | ✗(1.1E7) | ✓
$\mathbf{sinpi(x)}$ | ✓ | N/A | N/A | ✗(3.4E5) | ✓ | ✓ | N/A | N/A
$\mathbf{cospi(x)}$ | ✓ | N/A | N/A | ✗(3.8E5) | ✓ | ✓ | N/A | N/A

We provide details on our prototype, experimental methodology, and the results
of our experiments to check the correctness and performance of the generated
functions.

###  4.1. Experimental Setup and Methodology

Prototype. The RLibm-32 prototype generates correctly rounded elementary
functions for 32-bit floats and posit32, which is a 32-bit posit type
providing tapered precision (_i.e._ , more precision than float for values
near 1) (Gustafson, 2017). It contains ten correctly rounded elementary
functions for 32-bit floats and eight elementary functions for the posit32
type. To generate correctly rounded elementary functions with good
performance, the user can provide custom range reduction functions and specify
the degree or the structure of the polynomial (_i.e._ , odd or even). RLibm-32
uses the MPFR library (Fousse et al., 2007) with up to 400 precision bits to
compute the oracle for $f(x)$ and rounds it to the target representation,
which is good enough to compute the oracle result for double (Lefèvre and
Muller, 2001). RLibm-32 uses SoPlex (Gleixner et al., 2012), an exact rational
LP solver, for generating coefficients for the polynomials with a five minute
time limit. We use a threshold of fifty thousand reduced inputs and intervals
in the sample for counterexample guided polynomial generation. RLibm-32’s math
library performs range reduction, polynomial evaluation, and output
compensation using double precision. Polynomial evaluation uses the Horner’s
method (Borwein and Erdelyi, 1995). We designed novel extensions to range
reduction for many elementary functions, which is inspired by table-based
range reduction (Tang, 1989, 1990, 1991; Daramy-Loirat et al., 2006). The
appendix provides additional details about range reduction for each elementary
function. RLibm-32 is open source and publicly available (Lim and Nagarakatte,
2021).

Methodology. We test the elementary functions in RLibm-32 on two dimensions:
(1) ability to generate correct results and (2) performance in comparison to
state-of-the-art libraries. We compare RLibm-32’s functions with four
libraries: Intel’s libm, glibc’s libm, CR-LIBM (Daramy-Loirat et al., 2006),
and Metalibm (Kupriianova and Lauter, 2014). To use double precision
libraries, we convert the float input into double, use the double function,
and round the result back to float. Among these libraries, CR-LIBM has
correctly rounded functions for double precision. However, CR-LIBM does not
produce correctly rounded results for 32-bit floats due to double rounding.
There are no math libraries available for posit32. All posit32 values can be
exactly represented in double. Hence, we compare our posit32 library with
glibc and Intel’s double libm and CR-LIBM.

Experimental setup. We performed all our experiments on a 2.10GHz Intel Xeon
Gold 6230R machine with 187GB of RAM running Ubuntu 18.04. We disabled Intel
turbo boost and hyper-threading to minimize noise. We compiled RLibm-32’s math
library at the O3 optimization level. We used Intel’s libm from the oneAPI
Toolkit and glibc’s libm from glibc-2.33. We generated Metalibm
implementations with optimizations for AVX2 extensions enabled. Our test
harness that compares glibc’s libm, CR-LIBM, and Metalibm with RLibm-32 is
built using the gcc-10 compiler with -O3 -static -frounding-math -fsignaling-
nans flags. To use Intel’s libm, we have to use the Intel compiler. Hence, the
test harness that compares Intel libm with RLibm-32 is built using the icc
compiler with -O3 -no-ftz -fp-model strict -static to obtain as many correct
results as possible. Further, the size of the executable generated by
statically linking RLibm-32 is 2% smaller on average when compared to the
executable generated with Intel’s double libm.

Measuring performance. To compare performance, we measure the number of cycles
taken to compute the result for each input using hardware performance
counters. The total time taken is computed as the sum of the time taken by all
inputs (_i.e._ , all $2^{32}$ inputs for a 32-bit representation). We ran the
measurements for all inputs for each function six times. Then, we compute the
average time taken to compute each elementary function. As Intel’s compiler
performs vectorization by default at the O3 optimization level, our above
setup does not measure improvements due to vectorization. Hence, we created
another test harness that creates an array of 1024 floats (_i.e._ , $2^{10}$
inputs), populates it with different inputs, and measures the number of cycles
taken to compute the results of $2^{10}$ inputs using hardware performance
counters. We repeat this experiment $2^{22}$ times to compute the result and
measure the total time taken for all $2^{32}$ inputs.

###  4.2. Generation of Correctly Rounded Results

Table 1 reports the results of our experiments to check the correctness of
various elementary functions in RLibm-32 and other mainstream libraries.
RLibm-32 produces the correctly rounded results for all inputs for the ten
elementary functions for 32-bit floats. In contrast, elementary functions in
glibc, Intel, and MetaLibm’s float library do not produce the correct result
for all inputs. Multiple functions in glibc and MetaLibm’s float library
produce wrong results for several million inputs. Intel’s libm also produces
wrong results with several thousand inputs with the float version. When we use
double precision version of functions from glibc, Intel’s libm, and CR-LIBM,
it does not produce the correct result for $ln(x)$, $log10(x)$, $exp2(x)$, and
$sinh(x)$. These cases occur when the real value of $f(x)$ is extremely close
to the rounding boundary of a floating point value. Even with a smaller mini-
max approximation error in the double library compared to their float
versions, these libraries do not produce the correctly rounded result for all
inputs. CR-LIBM, which is a correctly rounded double library, produces wrong
results for float functions due to double rounding. We observed that functions
in MetaLibm do not produce correct results even when it internally uses Sollya
(Chevillard et al., 2010), which can be used to generate correctly rounded
implementations.

Table 2. Generation of correctly rounded results with posit32 functions for all inputs by RLibm-32, Intel and glibc’s double libraries, and CR-LIBM. ✓indicates that the library produces the correctly rounded result for all inputs and otherwise, we use ✗. | posit32
---
functions
| Using
---
RLibm-32
| Using
---
glibc double
| Using
---
Intel double
| Using
---
CR-LIBM
$\mathbf{ln(x)}$ | ✓ | ✗(22) | ✗(22) | ✗(22)
$\mathbf{log2(x)}$ | ✓ | ✗(19) | ✗(18) | ✗(18)
$\mathbf{log10(x)}$ | ✓ | ✗(26) | ✗(23) | ✗(23)
$\mathbf{Exp(x)}$ | ✓ | ✗(4.4E8) | ✗(4.4E8) | ✗(4.4E8)
$\mathbf{Exp2(x)}$ | ✓ | ✗(4.0E8) | ✗(4.0E8) | N/A
$\mathbf{Exp10(x)}$ | ✓ | ✗(5.2E8) | ✗(5.2E8) | N/A
$\mathbf{Sinh(x)}$ | ✓ | ✗(4.4E8) | ✗(4.4E8) | ✗(4.4E8)
$\mathbf{Cosh(x)}$ | ✓ | ✗(4.4E8) | ✗(4.4E8) | ✗(4.4E8)

Table 2 reports that RLibm-32 produces correctly rounded results with all
inputs for the eight posit32 functions. All posit32 values are representable
in double precision but they cannot be represented in 32-bit floats. Hence, we
use CR-LIBM, Intel and glibc’s double library to compare with RLibm-32. These
libraries for double precision do not produce correct results for all posit32
inputs. Unlike functions for 32-bit floats, they produce wrong results for
several million inputs especially for exponential and hyperbolic functions.
One of the key reasons for wrong results is the absence of overflows to
$\infty$ and underflows to 0 with the posit32 type. Instead, extremely large
values are rounded to the largest representable value. Similarly, extremely
small values are rounded to the smallest non-zero representable value in the
posit32 type.

Table 3. Details about the generated polynomials. For each elementary function, time taken to generate the polynomials in minutes, the size of the piecewise polynomial for approximating $f_{i}(r)$, the maximum degree of the polynomial, and the number of terms in the polynomial. $f(x)$ |  | Gen. Time
---
(Minutes)
| Reduced
---
Inputs
| # of Poly-
---
nomials
| Deg-
---
ree
| # of
---
Terms
float functions
$ln(x)$ | $218$ | 7.2E6 | $2^{10}$ | 3 | 3
$log2(x)$ | $251$ | 7.2E6 | $2^{8}$ | 3 | 3
$log10(x)$ | $429$ | 7.2E6 | $2^{8}$ | 3 | 3
$exp(x)$ | $117$ | 5.2E8 |  | $2^{7}$
---
$2^{7}$
| 4
---
4
| 5
---
5
$exp2(x)$ | $86$ | 3.0E8 |  | $2^{4}$
---
$2^{3}$
| 4
---
4
| 5
---
5
$exp10(x)$ | $169$ | 5.2E8 |  | $2^{6}$
---
$2^{7}$
| 4
---
3
| 5
---
4
$sinh(x)$ | $28$ | 1.5E8 | $2^{6}$ | 5 | 3
$cosh(x)$ | $24$ | 1.5E8 | $2^{6}$ | 4 | 3
$sinpi(x)$ | $30$ | 1.2E8 | 1 | 5 | 3
$cospi(x)$ | $19$ | 1.2E8 | 1 | 4 | 3
posit32 functions
$ln(x)$ | 264 | 1.1E8 | $2^{11}$ | 4 | 4
$log2(x)$ | 288 | 1.1E8 | $2^{8}$ | 4 | 4
$log10(x)$ | 685 | 1.1E8 | $2^{12}$ | 3 | 3
$exp(x)$ | 1089 | 3.5E9 |  | $2^{12}$
---
$2^{12}$
| $3$
---
$3$
| $4$
---
$4$
$exp2(x)$ | 814 | 7.9E8 |  | $2^{10}$
---
$2^{12}$
| $3$
---
$3$
| $4$
---
$4$
$exp10(x)$ | 1528 | 3.4E9 |  | $2^{13}$
---
$2^{13}$
| $3$
---
$3$
| $4$
---
$4$
$sinh(x)$ | 461 | 1.6E9 |  | $2^{14}$
---
$2^{14}$
| $5$
---
$4$
| $3$
---
$3$
$cosh(x)$ | 528 | 1.7E9 |  | $2^{14}$
---
$2^{12}$
| $3$
---
$6$
| $2$
---
$4$

Piecewise polynomials generated by RLibm-32. Table 3 provides details on the
piecewise polynomials generated by RLibm-32. Our goal is to get the best
possible performance within a given storage budget for piecewise polynomials
(_i.e._ , number of sub-domains when we split the range of reduced inputs).
Hence, we used the RLibm-32 to generate piecewise polynomials such that the
degree of each polynomial was less than or equal to 8 and the number of sub-
domains was less than or equal to $2^{14}$. The output compensation for
$sinh(x)$, $cosh(x)$, $sinpi(x)$, and $cospi(x)$ involves two elementary
functions. We generate two piecewise polynomials for each of those elementary
functions. There are both positive and negative reduced inputs for $exp(x)$,
$exp2(x)$, and $exp10(x)$. Hence, we created two piecewise polynomials: one
for the negative reduced inputs and another for positive reduced inputs.
Notably, we were able to generate a single polynomial of degree 5 and 4 that
satisfies all reduced constraints for $sinpi(r)$ and $cospi(r)$, respectively.
Both $sinpi(r)$ and $cospi(r)$ have close to 120 million reduced inputs. Our
counterexample driven polynomial generation with sampling was instrumental in
creating this efficient polynomial.

Time taken to generate RLibm-32 functions. Table 3 also reports the time taken
to generate the 32-bit float and the posit32 functions in RLibm-32. It ranges
from 19 minutes for $cospi(x)$ for the float type to approximately 25 hours
for $exp10(x)$ for the posit32 type. Majority of the total time total time is
spent in computing the oracle result and the rounding interval using the MPFR
library (_i.e._ , 86% of total time for 32-bit floats and 55% of total time
for the posit32 type). In contrast, counterexample guided polynomial
generation takes 14% and 45% of the total time for 32-bit floats and the
posit32 type, respectively. We noticed that it takes significantly longer to
generate posit32 functions. There are fewer special cases, which requires
longer oracle computation. Further, RLibm-32 generates larger piecewise
polynomials for posit32 functions to account for higher precision than a
32-bit float and saturating behavior with extremal values.

(a) Speedup of RLibm-32’s float functions over glibc libm

(b) Speedup of RLibm-32’s float functions over Intel libm

(c) Speedup of RLibm-32’s float functions over CR-LIBM

(d) Speedup of RLibm-32’s float functions over MetaLibm

Figure 3. (a) Speedup of RLibm-32’s float functions compared to glibc’s float
functions (left) and glibc’s double functions (right). (b) Speedup of
RLibm-32’s functions compared to Intel’s float functions (left) and Intel’s
double functions (right). (c) Speedup of RLibm-32’s functions compared to CR-
LIBM functions. (d) Speedup of RLibm-32’s functions compared to MetaLibm’s
float functions (left) and double functions (right) built with AVX2
optimizations.

###  4.3. Performance Evaluation of RLibm-32

Performance of float functions. Figure 3(a) presents the speedup of RLibm-32’s
float functions over glibc’s float functions (left bar in each cluster) and
double functions (right bar in each cluster). On average, RLibm-32’s float
functions have $1.1\times$ speedup over glibc’s float libm and $1.2\times$
speedup over glibc’s double libm. Figure 3(b) reports the speedup of
RLibm-32’s float functions over Intel’s float libm and double libm. RLibm-32’s
float functions have an average of $1.5\times$ speedup over Intel’s float
functions and $1.6\times$ speedup over Intel’s double functions. Figure 3(c)
reports that RLibm-32’s functions are on average $2\times$ faster than CR-LIBM
functions. Figure 3(d) reports the speedup of RLibm-32’s functions over
MetaLibm’s float and double functions. RLibm-32’s functions are on average
$2.5\times$ and $2.7\times$ faster than MetaLibm’s float and double functions,
respectively. RLibm-32’s functions are faster than all the corresponding
functions in Intel libm, CR-LIBM, and MetaLibm. RLibm-32’s functions are
faster than glibc’s functions except for $ln(x)$, $log_{2}(x)$, and
$log_{10}(x)$ for float and $ln(x)$ for double. However, glibc’s libm produces
a large number of wrong results for them. RLibm-32’s functions are not only
faster but also produce correctly rounded results for all inputs.

(a) Speedup of RLibm-32’s posit32 functions over glibc libm

(b) Speedup of RLibm-32’s posit32 functions over Intel libm

(c) Speedup of RLibm-32’s posit32 functions over CR-LIBM

Figure 4. (a) Speedup of RLibm-32’s posit32 functions compared to glibc’s
double functions. (b) Speedup of RLibm-32’s posit32 functions compared to
Intel’s double functions. (c) Speedup of RLibm-32’s posit32 functions compared
to CR-LIBM functions.

Performance of posit32 functions. The graphs in Figure 4(a), Figure 4(b), and
Figure 4(c) report the speedup of RLibm-32’s posit32 functions when compared
to math libraries created by re-purposing glibc’s, Intel’s, and CR-LIBM’s
double functions, respectively. On average, RLibm-32’s posit32 functions are
$1.1\times$, $1.1\times$, and $1.4\times$ faster than glibc’s libm, Intel’s
libm, and CR-LIBM, respectively. All three re-purposed math libraries produce
wrong results for some inputs. RLibm-32 provides the first correctly rounded
functions for the posit32 type.

Vectorization. Intel compiler uses vector instructions to improve performance
by default. In our experiments with vectorization using an array of 1024
inputs (see Section 4.1), RLibm-32 is on average 10% and 5% slower than
Intel’s float libm and double libm, respectively. However, Intel’s compiler
produces wrong results for several million inputs (without -no-ftz -fp-model
strict flags). In contrast, RLibm-32’s functions are almost as fast as
vectorized code while producing correct results for all inputs.

Figure 5. Performance speedup of $log2(x)$ and $log10(x)$ with an increase in
the number of sub-domains when compared to a single polynomial generated by
RLibm-32. All these polynomials produce the correctly rounded result for all
inputs. A circle represents a decrease in the degree of the piecewise
polynomial.

Performance impact of piecewise polynomials. To analyze the performance
benefits due to piecewise polynomials, we identified elementary functions for
which we could generate a single polynomial that produces correctly rounded
results for all inputs ($log2(x)$, $log10(x)$, $sinpi$, and $cospi$) . We
measured the change in performance with an increase in the number of sub-
domains ranging from $2^{0}$ (_i.e._ , a single polynomial) to $2^{12}$.
Figure 5 reports the performance of $log2(x)$ and $log_{10}(x)$ with an
increase in the number of sub-domains when compared to the performance of a
single polynomial. We validated that all these polynomials produce the correct
result for all inputs. Figure 5 does not report $sinpi$ and $cospi$ because
the single polynomial has the best performance. Initially, there is a small
decrease in performance by moving from a single polynomial to a piecewise
polynomial because the degree of the piecewise polynomial does not decrease
significantly to subsume the overhead of table lookup. On increasing the
number of sub-domains, we observed almost $1.2\times$ speedup with piecewise
polynomials having $2^{8}$ sub-domains. It requires 6KB for storing
coefficients of piecewise polynomials.

##  5\. Case Study with $cospi(x)$ for Float

We describe the case study with $cospi(x)=cos(\pi x)$ to illustrate the
importance of carefully designing range reduction to avoid cancellation errors
in output compensation. The elementary function $cospi(x)$ is defined for
$x\in(-\infty,\infty)$.

Special cases. There are three kinds of special cases:

$cospi(x)=\begin{cases}1.0&\text{if }|x|<7.771\times 10^{-5}\\\
(-1)^{(|x|\>mod\>2)}\times 1.0&\text{if }|x|\geq 2^{23}\\\ NaN&\text{if
}x=NaN\text{ or }x=\pm\infty\\\ \end{cases}$

All float values $\geq 2^{23}$ are integers. Hence, $cospi(x)=1.0$ for even
integers and $cospi(x)=-1.0$ for odd integer inputs.

Range reduction of $cospi(x)$. After excluding special cases, there are more
than 600 million float inputs that need to be approximated. Similar to range
reduction for $sinpi(x)$ (Section 2.1), we use periodicity and trigonometric
identities of $cospi(x)$ to reduce inputs to a smaller domain. We transform
input $x$ into $x=2.0\times I+J$ where $I$ is an integer and $J\in[0,2)$. Due
to periodicity, $cospi(x)=cospi(J)$. Next, we decompose $J$ into $J=K+L$ where
$K$ is the integral part of $J$ ($K\in\\{0,1\\}$) and $L\in[0,1)$ is the
fractional part. Then, $cospi(J)$ can be computed with,

$cospi(J)=(-1)^{K}cospi(L)$

To further reduce the range of $L$, we use the fact that $cospi(x)$ between
$[0.5,1)$ is a mirror image of $cospi(x)$ between $[0,0.5)$ with the opposite
sign. We decompose $L$ into $M$ and $L^{\prime}$ where

$M=\begin{cases}0&\text{if }\leq 0.5\\\ 1&\text{if }L>0.5\end{cases}\quad\quad
L^{\prime}=\begin{cases}L&\text{if }L\leq 0.5\\\ 1.0-L&\text{if
}L>0.5\end{cases}$

We have $cospi(L)=(-1)^{M}cospi(L^{\prime})$. After reducing the input $x$ to
$L^{\prime}\in[0,0.5]$, there are around 107 million inputs. Thus, we further
reduce $L^{\prime}$ to a value in $[0,\frac{1}{512}]$. We split
$L^{\prime}=\frac{N}{512}+Q$ where $N$ is an integer in the set
$\\{0,1,2,\dots,255\\}$ and $Q$ is a fractional value in $[0,\frac{1}{512}]$.
One possible method to compute $cospi(L^{\prime})$ is to use the trigonometric
identity $cospi(a+b)=cospi(a)cospi(b)-sinpi(a)sinpi(b)$,

$cospi\left(\frac{N}{512}+Q\right)=cospi\left(\frac{N}{512}\right)cospi(Q)-sinpi\left(\frac{N}{512}\right)sinpi(Q)$

The above formula is not monotonic and can have cancellation errors if $N\neq
0$ (if $N=0$, then $cospi(L^{\prime})=cospi(Q)$).

Creating monotonic output compensation. If $N\neq 0$, we transform $N$ and $Q$
to $N^{\prime}$ and $R$ such that $L^{\prime}=\frac{N^{\prime}}{512}-R$ to
create a monotonic output compensation function:

$N^{\prime}=\begin{cases}0&\text{if }N=0\\\ N+1&otherwise\end{cases}\quad\quad
R=\begin{cases}Q&\text{if }N=0\\\ \frac{1}{512}-Q&otherwise\end{cases}$

Then, we can compute $cospi(L^{\prime})=cospi(\frac{N^{\prime}}{512}-R)$ using
the trigonometric identity $cospi(a-b)=cospi(a)cospi(b)+sinpi(a)sinpi(b)$ as
follows,

$\displaystyle cospi(L^{\prime})=$
$\displaystyle\begin{cases}cospi(R)&\text{if }N=0\\\
cospi\left(\frac{N^{\prime}}{512}\right)cospi(R)+sinpi\left(\frac{N^{\prime}}{512}\right)sinpi(R)&\text{if
}N\neq 0\end{cases}$

This output compensation is monotonic and does not experience cancellation
error. The values of $N^{\prime}$ ranges from $0$ to $256$ and
$R\in[0,\frac{1}{512}]$. The computation $\frac{1}{512}-Q$ can be computed
exactly with float or double type for all values of $Q$ that corresponds to
$N\neq 0$. There are approximately 40 million values of $R$. We precompute the
values for $cospi\left(\frac{N^{\prime}}{512}\right)$ and
$sinpi\left(\frac{N^{\prime}}{512}\right)$ in lookup tables (_i.e._ , 514
values in total). We create polynomial approximations for $sinpi(R)$ and
$cospi(R)$ for the reduced input domain $R\in[0,\frac{1}{512}]$. Using
RLibm-32, we were able to generate a single $5^{th}$ degree odd polynomial for
$sinpi(R)$ and a single $4^{th}$ degree even polynomial for $cospi(R)$.
Finally, we can compute the result for $cospi(x)$ with the output compensation
function,

$cospi(x)=\begin{cases}S\times cospi(R)&\text{if }N=0\\\
S\times\left(cpn\times cospi(R)+spn\times sinpi(R)\right)&\text{if }N\neq
0\end{cases}$

where $S=(-1)^{K}\times(-1)^{M}$, $cpn=cospi(\frac{N^{\prime}}{512})$, and
$spn=sinpi(\frac{N^{\prime}}{512})$. These polynomials combined with the
output compensation functions produce correctly rounded results for all inputs
for $sinpi(x)$ and $cospi(x)$.

##  6\. Related Work

Multiple decades of seminal work has advanced the state-of-the-art on creating
approximations for FP representations (Jeannerod et al., 2011; Bui and Tahar,
1999; Ziv, 1991; Fousse et al., 2007; Muller, 2005; Trefethen, 2012; Remes,
1934; De Caro et al., 2017). Further, seminal research on range reduction has
made such approximation feasible (Tang, 1990, 1991, 1989; Story and Ping Tak
Peter Tang, 1999; Cody and Waite, 1980; Boldo et al., 2009). Simultaneously,
there are verification efforts to prove bounds for math libraries (Harrison,
2009, 1997a, 1997b; Sawada, 2002; Lee et al., 2017), identify numerical errors
with expressions that can be used in the implementation of math libraries
(Chowdhary et al., 2020; Benz et al., 2012; Fu and Su, 2019; Goubault, 2001;
Sanchez-Stern et al., 2018), and repair individual outputs of math libraries
(Panchekha et al., 2015; Yi et al., 2019; Zou et al., 2019).

Correctly rounded libraries. Numerous groups have developed correctly rounded
elementary functions (Jeannerod et al., 2011; Bui and Tahar, 1999). Some
correctly rounded libraries for FP are IBM LibUltim (Ziv, 1991), Sun
Microsystem’s LibMCR, CR-LIBM (Daramy-Loirat et al., 2006), MPFR math library
(Fousse et al., 2007), and RLIBM (Lim et al., 2021, 2020a). CR-LIBM is a
correctly rounded double library developed using Sollya (Chevillard et al.,
2010), which generates mini-max polynomials to approximate elementary
functions (Brisebarre et al., 2006; Brisebarre and Chevillard, 2007). Sollya
uses the modified Remez algorithm (Remes, 1934) using lattice basis reduction
and also computes the error bound of the polynomial (Chevillard and Lauter,
2007; Chevillard et al., 2011; Melquiond, 2019). Metalibm (Kupriianova and
Lauter, 2014; Brunie et al., 2015) builds on Sollya and generates efficient
mini-max polynomials with user-defined error bounds. It also uses domain
splitting and hardware specific optimizations (Kupriianova and Lauter, 2015).
Compared to mini-max approaches, our work approximates the correctly rounded
result of $f(x)$ and generates polynomials that already account for numerical
error in range reduction and output compensation. Hence, it generates
efficient and correctly rounded results for all inputs.

This paper extends our prior work on RLibm (Lim et al., 2021, 2020a) and John
Gustafson’s Minefield method (Gustafson, 2020), which advocate approximating
the correctly rounded value rather than real value of an elementary function.
Our prior work on RLibm also frames the problem of generating polynomials as
an LP problem. We have used RLibm to create correctly rounded functions for
16-bit types: bfloat16 and posit16. This paper extends RLibm to handle 32-bit
types with systematic counterexample guided polynomial generation, generation
of piecewise polynomials to improve performance, and new techniques to deduce
rounding intervals when range reduction involves multiple elementary
functions.

Posit libraries. SoftPosit-Math (Leong, 2019) and RLibm libraries provide
correctly rounded math functions for 16-bit posits. In our prior work, we have
produced approximations for a set of trigonometric functions using the CORDIC
method for posit32 (Lim et al., 2020b). However, it does not produce correct
results for all inputs. In this paper, we develop the first set of elementary
functions that produce correctly rounded results for all inputs for 32-bit
posits.

##  7\. Conclusion and Future Directions

Mainstream math libraries have been designed and improved by numerous
researchers spanning multiple decades. Yet, they fail to generate correct
results for all inputs. This paper advocates approximating the correctly
rounded value instead of the real value similar to our prior work on RLibm. It
extends RLibm to scale to 32-bit representations: (a) counterexample guided
polynomial generation with an LP solver to handle billions of inputs, (b)
generation of constraints to account for multiple elementary functions in
range reduction, and (c) generation of piecewise polynomials. The resulting
functions produce correct results for all inputs and are also faster than
existing libraries for 32-bit floats and posits.

Going forward, we plan to generate approximations for all commonly used
elementary functions with 32-bit types, which we believe can be accomplished
with our approach. However, it may require us to develop novel extensions to
range reduction. Further, it may be necessary to perform range reduction in
higher precision for some trigonometric functions such as sine and cosine that
use $\pi$. Beyond 32-bit types, we also plan to extend this approach to double
precision. Our approach can generate a polynomial that produces the correctly
rounded result for the sampled points in the double type. Validating the
correctness of the result produced by a polynomial generated using our
approach for all inputs in the double type is an open research problem. Our
long-term goal is to enable the standards of existing and new representations
to mandate correctly rounded results.

###### Acknowledgements.

We thank our shepherd Rahul Sharma and the PLDI reviewers for their feedback.
We thank John Gustafson for his inputs on the Minefield method and the posit
representation. This material is based upon work supported in part by the
Sponsor National Science Foundation http://dx.doi.org/10.13039/100000001 under
Grant No. Grant #1908798 and Grant No. Grant #1917897.

## References

  * (1)
  * Benz et al. (2012) Florian Benz, Andreas Hildebrandt, and Sebastian Hack. 2012\.  A Dynamic Program Analysis to Find Floating-point Accuracy Problems. In _Proceedings of the 33rd ACM SIGPLAN Conference on Programming Language Design and Implementation_ (Beijing, China) _(PLDI ’12)_. ACM, New York, NY, USA, 453–462.  https://doi.org/10.1145/2345156.2254118
  * Boldo et al. (2009) Sylvie Boldo, Marc Daumas, and Ren-Cang Li. 2009\.  Formally Verified Argument Reduction with a Fused Multiply-Add. In _IEEE Transactions on Computers_ , Vol. 58. 1139–1145.  https://doi.org/10.1109/TC.2008.216
  * Borwein and Erdelyi (1995) Peter Borwein and Tamas Erdelyi. 1995.  _Polynomials and Polynomial Inequalities_.  Springer New York.  https://doi.org/10.1007/978-1-4612-0793-1
  * Brisebarre and Chevillard (2007) Nicolas Brisebarre and Sylvvain Chevillard. 2007.  Efficient polynomial L-approximations. In _18th IEEE Symposium on Computer Arithmetic (ARITH ’07)_.  https://doi.org/10.1109/ARITH.2007.17
  * Brisebarre et al. (2006) Nicolas Brisebarre, Jean-Michel Muller, and Arnaud Tisserand. 2006.  Computing Machine-Efficient Polynomial Approximations. In _ACM ACM Transactions on Mathematical Software_ , Vol. 32. Association for Computing Machinery, New York, NY, USA, 236–256.  https://doi.org/10.1145/1141885.1141890
  * Brunie et al. (2015) Nicolas Brunie, Florent de Dinechin, Olga Kupriianova, and Christoph Lauter. 2015.  Code Generators for Mathematical Functions. In _2015 IEEE 22nd Symposium on Computer Arithmetic_. 66–73.  https://doi.org/10.1109/ARITH.2015.22
  * Bui and Tahar (1999) Hung Tien Bui and Sofiene Tahar. 1999.  Design and synthesis of an IEEE-754 exponential function. In _Engineering Solutions for the Next Millennium. 1999 IEEE Canadian Conference on Electrical and Computer Engineering_ , Vol. 1. 450–455 vol.1.  https://doi.org/10.1109/CCECE.1999.807240
  * Chevillard et al. (2011) Sylvain Chevillard, John Harrison, Mioara Joldes, and Christoph Lauter. 2011.  Efficient and accurate computation of upper bounds of approximation errors. In _Theoretical Computer Science_ , Vol. 412.  https://doi.org/10.1016/j.tcs.2010.11.052
  * Chevillard et al. (2010) Sylvain Chevillard, Mioara Joldes, and Christoph Lauter. 2010\.  Sollya: An Environment for the Development of Numerical Codes. In _Mathematical Software - ICMS 2010_ _(Lecture Notes in Computer Science, Vol. 6327)_. Springer, Heidelberg, Germany, 28–31.  https://doi.org/10.1007/978-3-642-15582-6_5
  * Chevillard and Lauter (2007) Sylvain Chevillard and Christopher Lauter. 2007.  A Certified Infinite Norm for the Implementation of Elementary Functions. In _Seventh International Conference on Quality Software (QSIC 2007)_. 153–160.  https://doi.org/10.1109/QSIC.2007.4385491
  * Chowdhary et al. (2020) Sangeeta Chowdhary, Jay P. Lim, and Santosh Nagarakatte. 2020\.  Debugging and Detecting Numerical Errors in Computation with Posits. In _41st ACM SIGPLAN Conference on Programming Language Design and Implementation_ _(PLDI’20)_.  https://doi.org/10.1145/3385412.3386004
  * Cody and Waite (1980) William J Cody and William M Waite. 1980.  _Software manual for the elementary functions_.  Prentice-Hall, Englewood Cliffs, NJ.  https://doi.org/10.1137/1024023
  * Daramy-Loirat et al. (2006) Catherine Daramy-Loirat, David Defour, Florent de Dinechin, Matthieu Gallet, Nicolas Gast, Christoph Lauter, and Jean-Michel Muller. 2006.  _CR-LIBM A library of correctly rounded elementary functions in double-precision_.  Research Report. Laboratoire de l’Informatique du Parallélisme.  https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804
  * De Caro et al. (2017) Davide De Caro, Ettore Napoli, Darjn Esposito, Gerardo Castellano, Nicola Petra, and Antonio G. M. Strollo. 2017.  Minimizing Coefficients Wordlength for Piecewise-Polynomial Hardware Function Evaluation With Exact or Faithful Rounding.  _IEEE Transactions on Circuits and Systems I: Regular Papers_ (2017).  https://doi.org/10.1109/TCSI.2016.2629850
  * Fousse et al. (2007) Laurent Fousse, Guillaume Hanrot, Vincent Lefèvre, Patrick Pélissier, and Paul Zimmermann. 2007\.  MPFR: A Multiple-precision Binary Floating-point Library with Correct Rounding.  _ACM Trans. Math. Software_ 33, 2, Article 13 (June 2007).  https://doi.org/10.1145/1236463.1236468
  * Fu and Su (2019) Zhoulai Fu and Zhendong Su. 2019.  Effective Floating-point Analysis via Weak-distance Minimization. In _Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation_ (Phoenix, AZ, USA) _(PLDI 2019)_. ACM, New York, NY, USA, 439–452.  https://doi.org/10.1145/3314221.3314632
  * Gleixner et al. (2012) Ambros M. Gleixner, Daniel E. Steffy, and Kati Wolter. 2012\.  Improving the Accuracy of Linear Programming Solvers with Iterative Refinement. In _Proceedings of the 37th International Symposium on Symbolic and Algebraic Computation_ (Grenoble, France) _(ISSAC ’12)_. Association for Computing Machinery, New York, NY, USA, 187–194.  https://doi.org/10.1145/2442829.2442858
  * Goubault (2001) Eric Goubault. 2001\.  Static Analyses of the Precision of Floating-Point Operations. In _Proceedings of the 8th International Symposium on Static Analysis_ _(SAS)_. Springer, 234–259.  https://doi.org/10.1007/3-540-47764-0_14
  * Gustafson (2017) John Gustafson. 2017\.  _Posit Arithmetic_.  https://posithub.org/docs/Posits4.pdf
  * Gustafson (2020) John Gustafson. 2020\.  _The Minefield Method: A Uniformly Fast Solution to the Table-Maker’s Dilemma_.  https://bit.ly/2ZP4kHj
  * Harrison (1997a) John Harrison. 1997a.  Floating point verification in HOL light: The exponential function. In _Algebraic Methodology and Software Technology_ , Michael Johnson (Ed.). Springer Berlin Heidelberg, Berlin, Heidelberg, 246–260.  https://doi.org/10.1007/BFb0000475
  * Harrison (1997b) John Harrison. 1997b.  Verifying the Accuracy of Polynomial Approximations in HOL. In _International Conference on Theorem Proving in Higher Order Logics_.  https://doi.org/10.1007/BFb0028391
  * Harrison (2009) John Harrison. 2009\.  HOL Light: An Overview. In _Proceedings of the 22nd International Conference on Theorem Proving in Higher Order Logics, TPHOLs 2009_ _(Lecture Notes in Computer Science, Vol. 5674)_ , Stefan Berghofer, Tobias Nipkow, Christian Urban, and Makarius Wenzel (Eds.). Springer-Verlag, Munich, Germany, 60–66.  https://doi.org/10.1007/978-3-642-03359-9_4
  * Jeannerod et al. (2011) Claude-Pierre Jeannerod, Hervé Knochel, Christophe Monat, and Guillaume Revy. 2011.  Computing Floating-Point Square Roots via Bivariate Polynomial Evaluation.  _IEEE Trans. Comput._ 60\.  https://doi.org/10.1109/TC.2010.152
  * Kupriianova and Lauter (2014) Olga Kupriianova and Christoph Lauter. 2014.  Metalibm: A Mathematical Functions Code Generator. In _4th International Congress on Mathematical Software_.  https://doi.org/10.1007/978-3-662-44199-2_106
  * Kupriianova and Lauter (2015) Olga Kupriianova and Christoph Lauter. 2015.  Replacing Branches by Polynomials in Vectorizable Elementary Functions. In _Scientific Computing, Computer Arithmetic, and Validated Numerics_ , Marco Nehmeier, Jürgen Wolff von Gudenberg, and Warwick Tucker (Eds.). Springer International Publishing, Cham, 14–22.  https://doi.org/10.1007/978-3-319-31769-4_2
  * Lee et al. (2017) Wonyeol Lee, Rahul Sharma, and Alex Aiken. 2017\.  On Automatically Proving the Correctness of Math.h Implementations.  _Proceedings of the ACM on Programming Languages_ 2, POPL, Article 47 (Dec. 2017), 32 pages.  https://doi.org/10.1145/3158135
  * Lefèvre and Muller (2001) Vincent Lefèvre and Jean-Michel Muller. 2001.  Worst Cases for Correct Rounding of the Elementary Functions in Double Precision. In _15th IEEE Symposium on Computer Arithmetic_ _(Arith ’01)_. 111–118.  https://doi.org/10.1109/ARITH.2001.930110
  * Lefèvre et al. (1998) Vincent Lefèvre, Jean-Michel Muller, and Arnaud Tisserand. 1998.  Toward correctly rounded transcendentals.  _IEEE Trans. Comput._ 47, 11 (1998), 1235–1243.  https://doi.org/10.1109/12.736435
  * Leong (2019) Cerlane Leong. 2019\.  _SoftPosit-Math_.  https://gitlab.com/cerlane/softposit-math
  * Lim et al. (2020a) Jay P. Lim, Mridul Aanjaneya, John Gustafson, and Santosh Nagarakatte. 2020a.  A Novel Approach to Generate Correctly Rounded Math Libraries for New Floating Point Representations.  arXiv:2007.05344  Rutgers Department of Computer Science Technical Report DCS-TR-753.
  * Lim et al. (2021) Jay P. Lim, Mridul Aanjaneya, John Gustafson, and Santosh Nagarakatte. 2021.  An Approach to Generate Correctly Rounded Math Libraries for New Floating Point Variants.  _Proceedings of the ACM on Programming Languages_ 6, POPL, Article 29 (Jan. 2021), 30 pages.  https://doi.org/10.1145/3434310
  * Lim and Nagarakatte (2021) Jay P. Lim and Santosh Nagarakatte. 2021.  _RLibm-32_.  https://github.com/rutgers-apl/rlibm-32
  * Lim et al. (2020b) Jay P. Lim, Matan Shachnai, and Santosh Nagarakatte. 2020b.  Approximating Trigonometric Functions for Posits Using the CORDIC Method. In _Proceedings of the 17th ACM International Conference on Computing Frontiers_ (Catania, Sicily, Italy) _(CF ’20)_. Association for Computing Machinery, New York, NY, USA, 19–28.  https://doi.org/10.1145/3387902.3392632
  * Melquiond (2019) Guillaume Melquiond. 2019\.  _Gappa_.  http://gappa.gforge.inria.fr
  * Muller (2005) Jean-Michel Muller. 2005\.  _Elementary Functions: Algorithms and Implementation_.  Birkhauser.  https://doi.org/10.1007/978-1-4899-7983-4
  * Panchekha et al. (2015) Pavel Panchekha, Alex Sanchez-Stern, James R. Wilcox, and Zachary Tatlock. 2015.  Automatically Improving Accuracy for Floating Point Expressions. In _Proceedings of the 36th ACM SIGPLAN Conference on Programming Language Design and Implementation_ , Vol. 50. Association for Computing Machinery, New York, NY, USA, 1–11.  https://doi.org/10.1145/2813885.2737959
  * Remes (1934) Eugene Remes. 1934\.  Sur un procédé convergent d’approximations successives pour déterminer les polynômes d’approximation.  _Comptes rendus de l’Académie des Sciences_ 198 (1934), 2063–2065.
  * Sanchez-Stern et al. (2018) Alex Sanchez-Stern, Pavel Panchekha, Sorin Lerner, and Zachary Tatlock. 2018.  Finding Root Causes of Floating Point Error. In _Proceedings of the 39th ACM SIGPLAN Conference on Programming Language Design and Implementation_ (Philadelphia, PA, USA) _(PLDI 2018)_. ACM, New York, NY, USA, 256–269.  https://doi.org/10.1145/3296979.3192411
  * Sawada (2002) Jun Sawada. 2002\.  Formal verification of divide and square root algorithms using series calculation. In _3rd International Workshop on the ACL2 Theorem Prover and its Applications_.
  * Sterbenz (1974) Pat H Sterbenz. 1974\.  _Floating-point computation_.  Prentice-Hall, Englewood Cliffs, NJ.
  * Story and Ping Tak Peter Tang (1999) Shane Story and Ping Tak Peter Tang. 1999.  New algorithms for improved transcendental functions on IA-64. In _Proceedings 14th IEEE Symposium on Computer Arithmetic_. 4–11.  https://doi.org/10.1109/ARITH.1999.762822
  * Tang (1989) Ping-Tak Peter Tang. 1989\.  Table-Driven Implementation of the Exponential Function in IEEE Floating-Point Arithmetic.  _ACM Trans. Math. Software_ 15, 2 (June 1989), 144–157.  https://doi.org/10.1145/63522.214389
  * Tang (1990) Ping-Tak Peter Tang. 1990\.  Table-Driven Implementation of the Logarithm Function in IEEE Floating-Point Arithmetic.  _ACM Trans. Math. Software_ 16, 4 (Dec. 1990), 378–400.  https://doi.org/10.1145/98267.98294
  * Tang (1991) P. T. P. Tang. 1991\.  Table-lookup algorithms for elementary functions and their error analysis. In _[1991] Proceedings 10th IEEE Symposium on Computer Arithmetic_. 232–236.  https://doi.org/10.1109/ARITH.1991.145565
  * Trefethen (2012) Lloyd N. Trefethen. 2012\.  _Approximation Theory and Approximation Practice (Other Titles in Applied Mathematics)_.  Society for Industrial and Applied Mathematics, USA.
  * Yi et al. (2019) Xin Yi, Liqian Chen, Xiaoguang Mao, and Tao Ji. 2019\.  Efficient Automated Repair of High Floating-Point Errors in Numerical Libraries.  _Proceedings of the ACM on Programming Languages_ 3, POPL, Article 56 (Jan. 2019), 29 pages.  https://doi.org/10.1145/3290369
  * Ziv (1991) Abraham Ziv. 1991\.  Fast Evaluation of Elementary Mathematical Functions with Correctly Rounded Last Bit.  _ACM Trans. Math. Software_ 17, 3 (Sept. 1991), 410–423.  https://doi.org/10.1145/114697.116813
  * Zou et al. (2019) Daming Zou, Muhan Zeng, Yingfei Xiong, Zhoulai Fu, Lu Zhang, and Zhendong Su. 2019\.  Detecting Floating-Point Errors via Atomic Conditions.  _Proceedings of the ACM on Programming Languages_ 4, POPL, Article 60 (Dec. 2019), 27 pages.  https://doi.org/10.1145/3371128

##  Appendix A Details on Range Reduction used in RLibm-32

In this section, we describe the special cases, the range reduction, and the
output compensation functions that we used to create the math library
functions in RLibm-32. A given family of elementary functions (_e.g._ ,
$ln(x)$, $log_{2}(x)$, and $log_{10}(x)$) uses similar approaches for range
reduction. Hence, we first describe the range reduction technique we used for
each family of elementary functions. In the subsequent sections, we describe
the special cases and specific details on range reduction and output
compensation that we used for float functions (Appendix B) and posit32
functions (Appendix C).

###  A.1. Logarithm Functions $log_{a}(x)$

We use a table-based approach (Tang, 1990) to perform range reduction for
$log_{a}(x)$ functions. The range reduction is generally applicable to any
value of $a$. To reduce the input $x$ to the reduced input in a smaller
domain, we first transform input $x$ into $x=m\times 2^{n}$ where $n$ is an
integer representing the exponent of $x$ and $m\in[0,1)$ is the significand.
Then, $log_{a}(x)$ can be computed with,

$log_{a}(x)=log_{a}(m\times 2^{n})=log_{a}(m)+nlog_{a}(2)$

Next, we further reduce the range of $m$ by transforming $m$ into $m=F+f$
where $F=1+\frac{J}{128}$, $J$ is a value in the set $\\{0,1,2,\dots,127\\}$,
and $f\in[0,\frac{1}{128})$ is the remaining value. Intuitively, F is the
value represented by the first 8 bits of the significand $m$ and $f$ is the
value represented by the rest of the bits. Then, $log_{a}(m)$ can be computed
with,

$\displaystyle log_{a}(m)$
$\displaystyle=log_{a}\left(F+f\right)=log_{a}\left(F\>\left(1+\frac{f}{F}\right)\right)$
$\displaystyle=log_{a}(F)+log_{a}\left(1+\frac{f}{F}\right)$

We can denote $r=\frac{f}{F}$. Then, $log_{a}(m)$ can be computed as,

$log_{a}(m)=log_{a}(F)+log_{a}\left(1+r\right)$

The reduced input $r$ is in the range of $[0,\frac{1}{128})$. The computation
$r=\frac{f}{F}$ can be efficiently performed by computing $f\times\frac{1}{F}$
if the value $\frac{1}{F}$ can be computed ahead of time. Thus, we pre-compute
the values for $log_{a}(F)$ (128 values for each $log_{a}$) and $\frac{1}{F}$
(128 values in total) in lookup tables. We approximate $log_{a}(1+r)$ for the
reduced input domain $r\in[0,\frac{1}{128})$. Finally, we can compute the
result of $log_{a}(x)$ with the output compensation function,

$log_{a}(x)=log_{a}(1+r)+log_{a}(F)+nlog_{a}(2)$

###  A.2. Exponential Functions $a^{x}$

We also use a table-based approach (Tang, 1989) to perform range reduction for
$a^{x}$ functions. The range reduction is applicable for different values of
$a$. To reduce the input $x$ we first transform $x$ into
$x=nlog_{a}(2)+\frac{j}{64}log_{a}(2)+r$ where $n$ is an integer, $j$ is a
value in a set $\\{0,1,2,\dots,63\\}$ and $|r|\leq\frac{log_{a}(2)}{64}$.
Intuitively, $n$ represents the integral part of the value
$\frac{x}{log_{a}(2)}$. The value $\frac{j}{64}$ represents the first 6
fractional bits of $\frac{x}{log_{a}(2)}$. Lastly, $r$ represents the
remaining value, $r=x-nlog_{a}(2)-\frac{j}{64}log_{a}(2)$.

The scaling by $log_{a}(2)$ allows us to create efficient output compensation
formula. The value of $a^{x}$ can be computed using the property
$a^{x+y}=a^{x}a^{y}$ and $a^{xlog_{a}(y)}=y^{x}$ as,

$\displaystyle a^{x}=a^{nlog_{a}(2)+\frac{j}{64}log_{a}(2)+r}$
$\displaystyle=a^{nlog_{a}(2)}\times a^{\frac{j}{64}log_{a}(2)}\times a^{r}$
$\displaystyle=2^{n}\times 2^{\frac{j}{64}}\times a^{r}$

Multiplication by $2^{n}$ can be computed efficiently using bit-wise
operations. We pre-compute and store the value of $2^{\frac{j}{64}}$ in a
table (_i.e._ 64 values in total) and we approximate $a^{r}$ for the input
domain of $r\in[-\frac{log_{a}(2)}{64},\frac{log_{a}(2)}{64}]$.

###  A.3. Sinpi(x)

The range reduction technique for $sinpi(x)$ first leverages the periodicity
of $sinpi(x)$ to reduce the input and then use a table-based range reduction
(Tang, 1991) to further reduce the input. First, we transform input $x$ into
$x=2i+j$ where $i$ is an integer and $j\in[0,2)$. Then, $sinpi(x)=sinpi(j)$.
Next, we decompose $j$ into $j=k+l$ where $k\in\\{0,1\\}$ is the integral part
of $j$ and $L\in[0,1)$ is the fractional part. Then, $sinpi(j)$ can be
computed with,

$sinpi(j)=(-1)^{k}\times sinpi(l)$

Third, we use the fact that $sinpi$ between $[0.5,1)$ is a mirror image of
values between $[0,0.5]$ and decompose $l$ into,

$l^{\prime}=\begin{cases}l&\text{if }l\leq 0.5\\\ 1.0-l&\text{if
}l>0.5\end{cases}$

Then, $sinpi(l)$ can be computed with $sinpi(l)=sinpi(l^{\prime})$.

We further reduce $l^{\prime}$ to a value between $[0,\frac{1}{512}]$ using
table-based range reduction. We split $l^{\prime}$ into
$l^{\prime}=\frac{n}{512}+r$ where $n$ is an integer in the set
$\\{0,1,\dots,255\\}$ and $r$ is a real value in $[0,\frac{1}{512}]$. The
value of $sinpi(l^{\prime})$ can be computed using the trigonometric identity
$sinpi(a+b)=sinpi(a)cospi(b)+cospi(a)+sinpi(b)$:

$sinpi(l^{\prime})=sinpi\left(\frac{n}{512}\right)cospi(r)+cospi\left(\frac{n}{512}\right)sinpi(r)$

We pre-compute and store the values of $sinpi\left(\frac{n}{512}\right)$ and
$cospi\left(\frac{n}{512}\right)$ in lookup tables (_i.e._ 512 values). We
approximate $sinpi(r)$ and $cospi(r)$ for the reduced input domain
$r\in[0,\frac{1}{512}]$. Finally, we can approximate $sinpi(x)$ using the
output compensation function,

$sinpi(x)=(-1)^{k}\times\left(sinpi\left(\frac{n}{512}\right)cospi(r)+cospi\left(\frac{n}{512}\right)sinpi(r)\right)$

###  A.4. Cospi(x)

Similar to $sinpi(x)$, we leverage the periodicity of $cospi(x)$ and table-
based range reduction (Tang, 1991) to reduce the input. Additionally, we apply
some modifications to create monotonic output compensation function. First, we
transform the input $x$ into $x=2i+j$ where $i$ is an integer and $j\in[0,2)$.
Due to periodicity, $cospi(x)=cospi(j)$. Second, we decompose $j$ into $j=k+l$
where $k\in\\{0,1\\}$ is the integral part of j and $l\in[0,1)$ is the
fractional part. Then, $cospi(j)$ can be computed with,

$cospi(j)=(-1)^{k}\times cospi(l)$

Third, we reduce $l$ using the fact that $cospi(x)$ between $[0.5,1)$ is a
mirror image of $cospi(x)$ between $[0,0.5]$ with the opposite sign. We
decompose $l$ into $m$ and $l^{\prime}$ where:

$m=\begin{cases}0&\text{if }l\leq 0.5\\\ 1&\text{if
}l>0.5\end{cases}\quad\quad l^{\prime}=\begin{cases}l&\text{if }l\leq 0.5\\\
1.0-l&\text{if }l>0.5\end{cases}$

Then, $cospi(l)$ can be computed with $cospi(l)=(-1)^{m}cospi(l^{\prime})$.

We further reduce $l^{\prime}$ to a value between $[0,\frac{1}{512}]$ using
table-based range reduction. We transform $l^{\prime}$ into $r$ and $n$ using,

$l^{\prime}=\begin{cases}r&\text{if }l^{\prime}<\frac{1}{512}\\\
\frac{n}{512}-r&otherwise\end{cases}$

where $n$ is an integer value in the set $\\{0,1,2,\dots,256\\}$ and $r$ is a
fractional value in $[0,\frac{1}{512}]$. Then, $cospi(l^{\prime})$ can be
computed with the trigonometric identity
$cospi(a-b)=cospi(a)copsi(b)+sinpi(a)sinpi(b)$,

$cospi(l^{\prime})=\begin{cases}cospi(r)&\text{if }l^{\prime}<\frac{1}{512}\\\
cpn\times cospi(r)+spn\times sinpi(r)&otherwise\end{cases}$

where $cpn=cospi\left(\frac{n}{512}\right)$ and
$spn=sinpi\left(\frac{n}{512}\right)$. This formula is monotonic for all
inputs $x$. We pre-compute and store the values of
$sinpi\left(\frac{n}{512}\right)$ and $cospi\left(\frac{n}{512}\right)$ in
lookup tables (_i.e._ 514 values). We approximate $sinpi(r)$ and $cospi(r)$
for the reduced input domain $r\in[0,\frac{1}{512}]$.

Finally, we compute the result of $cospi(x)$ with the output compensation
formula,

$cospi(x)=\begin{cases}s\times cospi(r)&\text{if }l^{\prime}<\frac{1}{512}\\\
s\times(cpn\times cospi(r)+spn\times sinpi(r))&otherwise\end{cases}$

where $s=(-1)^{k}\times(-1)^{m}$.

###  A.5. Sinh(x)

We use the range reduction technique similar to the table-based range
reduction technique used in CR_LIBM. First, the $sinh(x)$ function has a
property, $sinh(-x)=-sinh(x)$. Thus, the result of $sinh(x)$ for $x<0$ can be
derived by computing $-1\times sinh(|x|)$.

Next, we decompose $|x|$ into three parts:

$|x|=k{ln(2)}+\frac{j}{64}ln(2)+r$

where both $k$ and $j$ are integers, $k\geq 0$, $0\leq j<64$, and $r$ is a
value in $[0,\frac{ln(2)}{64})$. As we will show later, scaling by $ln(2)$
allows us to compute $sinh(k{ln(2)})$ and $cosh(k{ln(2)})$ efficiently with
minimal amount of error. If we denote $K=k{ln(2)}$ and $J=\frac{j}{64}ln(2)$,
then $sinh(|x|)$ can be computed using the hyperbolic identities,
$sinh(a+b)=sinh(a)cosh(b)+cosh(a)sinh(b)$ and
$cosh(a+b)=cosh(a)cosh(b)+sinh(a)sinh(b)$:

$\displaystyle sinh(|x|)$ $\displaystyle=SH\times cosh(r)+CH\times sinh(r)$
$\displaystyle SH$ $\displaystyle=sinh(K)\times cosh(J)+cosh(K)\times sinh(J)$
$\displaystyle CH$ $\displaystyle=cosh(K)\times cosh(J)+sinh(K)\times sinh(J)$

We store the values of $sinh(K)$, $cosh(K)$, $sinh(J)$, and $cosh(J)$ in
lookup tables. Storing $sinh(J)$ and $cosh(J)$ requires a total of 128 values.
The number of values to store for $sinh(K)$ and $cosh(K)$ depends on the
target representation. In the case of float, the values of $k$ does not exceed
$129$ because $|x|<130ln(2)$ for all non-special-case inputs. Thus, we store
$260$ values for $sinh(K)$ and $cosh(K)$.

Alternatively, we can choose to not store the values of $sinh(K)$ and
$cosh(K)$ by computing $sinh(K)$ and $cosh(K)$ manually using the following
properties,

$sinh(x)=\frac{e^{x}-e^{-x}}{2}\quad\quad cosh(x)=\frac{e^{x}+e^{-x}}{2}$

Then, $sinh(K)$ and $cosh(K)$ can be computed as follows,

$\displaystyle sinh(K)=sinh(kln(2))=2^{k-1}-2^{-k-1}$ $\displaystyle
cosh(K)=cosh(kln(2))=2^{k-1}+2^{-k-1}$

Although $sinh(K)$ and $cosh(K)$ cannot be exactly represented by double type
for all $K$, the correctly rounded value of $sinh(K)$ and $cosh(K)$ can be
computed efficiently using integer arithmetic, floating point subtraction, and
addition. We plan to incorporate this logic in our implementation of $sinh(x)$
to reduced the size of the lookup table.

We approximate $sinh(r)$ and $cosh(r)$ for the reduced inputs
$r\in[0,\frac{ln(2)}{64})$. Finally, $sinh(x)$ can be computed with the output
compensation function,

$sinh(x)=s\times\left(SH\times cosh(r)+CH\times sinh(r)\right)$

where $s$ is the sign of $x$.

###  A.6. Cosh(x)

The range reduction technique for $cosh(x)$ uses a similar technique as
$sinh(x)$ (described in Appendix A.5) to reduce the input $x$. Then it uses
the hyperbolic identities of $cosh(x)$ to perform output compensation. First,
the $cosh(x)$ function has a property, $cosh(-x)=cosh(x)$. Thus, the result of
$cosh(x)$ for $x<0$ can be derived by computing $cosh(|x|)$.

Next, we decompose $|x|$ into three parts similar to how we decompose the
input $|x|$ for $sinh(x)$:

$|x|=k{ln(2)}+\frac{j}{64}ln(2)+r$

Both $k$ and $j$ are integers, $k\geq 0$, $0\leq j<64$, and $r$ is a real
number value in $[0,\frac{ln(2)}{64})$. If we denote $K=k{ln(2)}$ and
$J=\frac{j}{64}ln(2)$, then $cosh(|x|)$ can be computed using the hyperbolic
identities, $sinh(a+b)=sinh(a)cosh(b)+cosh(a)sinh(b)$ and
$cosh(a+b)=cosh(a)cosh(b)+sinh(a)sinh(b)$:

$\displaystyle cosh(x)$ $\displaystyle=SH\times sinh(R)+CH\times cosh(R)$
$\displaystyle SH$ $\displaystyle=sinh(K)\times cosh(J)+cosh(K)\times sinh(J)$
$\displaystyle CH$ $\displaystyle=cosh(K)\times cosh(J)+sinh(K)\times sinh(J)$

We store the values of $sinh(K)$, $cosh(K)$, $sinh(J)$, and $cosh(J)$ in
lookup tables (_i.e._ 388 values). We approximate $sinh(r)$ and $cosh(r)$ for
the reduced inputs $r\in[0,\frac{ln(2)}{64})$. Finally, $cosh(x)$ can be
computed with the output compensation function,

$cosh(x)=SH\times sinh(r)+CH\times cosh(r)$

##  Appendix B Details on 32-bit Float Functions

In this section, we explain the 32-bit float functions in RLibm-32. We
describe the special cases and any specific details on the range reduction
technique used for each function.

###  B.1. $ln(x)$

The elementary function $ln(x)$ is defined over the input domain $(0,\infty)$.
There are four classes of special case inputs:

$ln(x)=\begin{cases}NaN&\text{if }x=NaN\\\ NaN&\text{if }x<0\\\
-\infty&\text{if }x=0\\\ \infty&\text{if }x=\infty\end{cases}$

We use the range reduction technique described in Appendix A.1 to decompose
$x$ into $n$, $F$, and $r$. In the case of $ln(x)$, the output compensation
can be mathematically computed as,

$ln(x)=ln(1+r)+ln(F)+n\times ln(2)$

To evaluate the output compensation function in double, we store the correctly
rounded value of $ln2=ln(2)$ in double. We order the operations in the
following way,

$ln(x)=(ln(1+r)+ln(F))+n\times ln2$

This reduces the amount of numerical error for the inputs $x$ where the
magnitude of $n$ is large. After range reduction, there were roughly 7.2
million reduced inputs. The reduced inputs are in the range of
$r\in[0,\frac{1}{128})$. We created an approximation of $ln(1+r)$ using
RLibm-32 to generate a piecewise polynomial with $2^{10}$ polynomials of
degree 3.

###  B.2. $log_{2}(x)$

The elementary function $log_{2}(x)$ is defined over the input domain
$(0,\infty)$. There are four classes of special case inputs:

$log_{2}(x)=\begin{cases}NaN&\text{if }x=NaN\\\ NaN&\text{if }x<0\\\
-\infty&\text{if }x=0\\\ \infty&\text{if }x=\infty\end{cases}$

We use the range reduction technique described in Appendix A.1 to decompose
$x$ into $n$, $F$, and $r$. In the case of $log_{2}(x)$, the output
compensation can be mathematically computed as,

$log_{2}(x)=log_{2}(1+r)+log_{2}(F)+n$

To evaluate the output compensation function in double, we order the
operations in the following way,

$log_{2}(x)=(log_{2}(1+r)+log_{2}(F))+n$

This reduces the amount of numerical error for the inputs $x$ where the
magnitude of $n$ is large. After range reduction, there were roughly 7.2
million reduced inputs. The reduced inputs are in the range of
$r\in[0,\frac{1}{128})$. We created an approximation of $log_{2}(x)$ using
RLibm-32 to generate a piecewise polynomial with $2^{8}$ polynomials of degree
3.

###  B.3. $log_{10}(x)$

The elementary function $log_{10}(x)$ is defined over the input domain
$(0,\infty)$. There are four classes of special case inputs:

$log_{10}(x)=\begin{cases}NaN&\text{if }x=NaN\\\ NaN&\text{if }x<0\\\
-\infty&\text{if }x=0\\\ \infty&\text{if }x=\infty\end{cases}$

We use the range reduction technique described in Appendix A.1 to decompose
$x$ into $n$, $F$, and $r$. In the case of $log_{10}(x)$, the output
compensation can be mathematically computed as,

$log_{10}(x)=log_{10}(1+r)+log_{10}(F)+n\times log_{10}(2)$

To generate efficient piecewise polynomial approximation, we store the value
of $log_{10}(2)$ in two double values, $t_{h}$ and $t_{l}$ such that
$t_{h}+t_{l}$ is a correctly rounded value of $log_{10}(2)$ with 106 precision
bits (each double value has 53 precision bits). The value $t_{h}$ stores the
higher 53 precision bits and $t_{l}$ stores the lower 53 precision bits. Then,
we evaluate the output compensation function in double with the following
order of operation:

$log_{10}(x)=((log_{10}(1+r)+n\times t_{l})+log_{10}(F))+n\times t_{h}$

This reduces the amount of numerical error for the inputs $x$ where the
magnitude of $n$ is large. After range reduction, there were roughly 7.2
million reduced inputs. The reduced inputs are in the range of
$r\in[0,\frac{1}{128})$. We created an approximation of $log_{10}(x)$ using
RLibm-32 to generate a piecewise polynomial with $2^{8}$ polynomials of degree
3.

###  B.4. $e^{x}$

The $e^{x}$ function is defined over the input domain $(-\infty,\infty)$.
There are 4 classes of special case inputs:

$e^{x}=\begin{cases}NaN&\text{if }x=NaN\\\ 0&\text{if }x\leq-103.97\dots\\\
1.0&\text{if }-2.98\dots\times 10^{-8}\leq x\leq 5.96\dots\times 10^{-8}\\\
\infty&\text{if }x\geq 88.72\end{cases}$

We use the range reduction technique described in Appendix A.2 to perform
range reduction and output compensation. After range reduction, there are
roughly 520 million reduced inputs in the domain
$r\in[-\frac{ln(2)}{64},\frac{ln(2)}{64}]$. We created an approximation
function for negative values of reduced inputs $r$ and an approximation
function for positive values of reduced inputs to efficiently split the
reduced input domain. Using RLibm-32, we generated a piecewise polynomial with
$2^{7}$ polynomials of degree 4 for negative reduced inputs and a piecewise
polynomial with $2^{7}$ polynomials of degree 4 for positive reduced inputs.

###  B.5. $2^{x}$

The $2^{x}$ function is defined over the input domain $(-\infty,\infty)$.
There are 4 classes of special case inputs:

$2^{x}=\begin{cases}NaN&\text{if }x=NaN\\\ 0&\text{if }x\leq-150\\\
1.0&\text{if }-4.29\dots\times 10^{-8}\leq x\leq 8.59\dots\times 10^{-8}\\\
\infty&\text{if }x\geq 128\end{cases}$

We use the range reduction technique described in Appendix A.2 to perform
range reduction and output compensation. After range reduction, there are
roughly 303 million reduced inputs in the domain
$r\in[-\frac{1}{64},\frac{1}{64}]$. We created an approximation function for
negative values of reduced inputs $r$ and an approximation function for
positive values of reduced inputs to efficiently split the reduced input
domain. Using RLibm-32, we generated a piecewise polynomial with $2^{4}$
polynomials of degree 4 for negative reduced inputs and a piecewise polynomial
with $2^{3}$ polynomials of degree 4 for positive reduced inputs.

###  B.6. $10^{x}$

The $10^{x}$ function is defined over the input domain $(-\infty,\infty)$.
There are 4 classes of special case inputs:

$10^{x}=\begin{cases}NaN&\text{if }x=NaN\\\ 0&\text{if }x\leq-45.15\dots\\\
1.0&\text{if }-1.29\dots\times 10^{-8}\leq x\leq 2.58\dots\times 10^{-8}\\\
\infty&\text{if }x\geq 38.53\dots\end{cases}$

We use the range reduction technique described in Appendix A.2 to perform
range reduction and output compensation. After range reduction, there are
roughly 521 million reduced inputs in the domain
$r\in[-\frac{log_{10}(2)}{64},\frac{log_{10}(2)}{64}]$. We created an
approximation function for negative values of reduced inputs $r$ and an
approximation function for positive values of reduced inputs to efficiently
split the reduced input domain. Using RLibm-32, we generated a piecewise
polynomial with $2^{6}$ polynomials of degree 4 for negative reduced inputs
and a piecewise polynomial with $2^{7}$ polynomials of degree 3 for positive
reduced inputs.

###  B.7. $Sinpi(x)$

The $sinpi(x)$ function is defined over the input domain $(-\infty,\infty)$.
There are three classes of special cases:

$sinpi(x)=\begin{cases}\pi x&\text{if }|x|<1.173\dots\times 10^{-7}\\\
0&\text{if }|x|\geq 2^{23}\\\ NaN&\text{if }x=NaN\text{ or
}x=\pm\infty\end{cases}$

We use the range reduction technique described in Appendix A.3 and evaluate
the range reduction and output compensation in double. After range reduction,
there were roughly 117 million reduced inputs in the domain
$r\in[0,\frac{1}{512}]$. The output compensation function for $sinpi(x)$ uses
both $sinpi(r)$ and $cospi(r)$. Similarly, the output compensation function
for $cospi(x)$ uses the approximation of both $sinpi(r)$ and $cospi(r)$. Thus,
we generated an approximation function for $sinpi(r)$ and an approximation
function for $cospi(r)$ that can be used to compute both $sinpi(x)$ and
$cospi(x)$. Using RLibm-32, we created a single polynomial of degree 5 for
$sinpi(r)$ and a single polynomial of degree 4 for $cospi(r)$.

###  B.8. $Cospi(x)$

The $cospi(x)$ function is defined over the input domain $(-\infty,\infty)$.
There are three classes of special cases:

$cospi(x)=\begin{cases}1.0&\text{if }|x|<7.771\times 10^{-5}\\\
(-1)^{(|x|\>mod\>2)}\times 1.0&\text{if }|x|\geq 2^{23}\\\ NaN&\text{if
}x=NaN\text{ or }x=\pm\infty\\\ \end{cases}$

We use the range reduction technique described in Appendix A.4. We evaluate
the range reduction and output compensation in double. After range reduction,
there were roughly 40 million reduced inputs in the domain
$r\in[0,\frac{1}{512}]$. The output compensation function for $cospi(x)$ uses
both $sinpi(r)$ and $cospi(r)$. These approximation functions are also used in
$sinpi(x)$. Thus, we generated approximation functions for $sinpi(r)$ and
$cospi(r)$ that can be used to compute the results for both $sinpi(x)$ and
$cospi(x)$. Using RLibm-32, we created a single polynomial of degree 5 for
$sinpi(r)$ and a single polynomial of degree 4 for $cospi(r)$.

###  B.9. $Sinh(x)$

The $sinh(x)$ function is defined over the input domain $(-\infty,\infty)$.
There are four classes of special cases:

$\text{sinh(x)}=\begin{cases}-\infty&\text{if }x<=-89.415\dots\\\ x&\text{if
}|x|\leq 4.436\dots\times 10^{-4}\\\ \infty&\text{if }x>=89.415\dots\\\
NaN&\text{if }x=NaN\end{cases}$

We use the range reduction technique described in Appendix A.5 and evaluate
the range reduction and output compensation functions in double. After range
reduction, there were roughly 147 million reduced inputs in the domain
$r\in[0,\frac{ln(2)}{64})$. Because the output compensation function for both
$sinh(x)$ and $cosh(x)$ uses approximations of $sinh(r)$ and $cosh(r)$, we
generated a piecewise polynomial for $sinh(r)$ and a piecewise polynomial for
$cosh(r)$ that can be used to compute both $sinh(x)$ and $cosh(x)$ correctly.
Using RLibm-32, we created a piecewise polynomial with $2^{6}$ polynomials of
degree 5 for $sinh(r)$ and a piecewise polynomial with $2^{6}$ polynomials of
degree 4 for $cosh(r)$.

###  B.10. $Cosh(x)$

The $cosh(x)$ functions is defined over the input domain $(-\infty,\infty)$.
There are three classes of special cases:

$\text{cosh(x)}=\begin{cases}1.0&\text{if }|x|\leq 3.452\dots\times 10^{-4}\\\
\infty&\text{if }|x|>=89.415\dots\\\ NaN&\text{if }x=NaN\end{cases}$

We use the range reduction technique described in Appendix A.6 and evaluate
the range reduction and output compensation functions in double. After range
reduction, there were roughly 151 million reduced inputs in the domain. We
approximate $sinh(r)$ and $cosh(r)$ for the reduced inputs in the domain
$r\in[0,\frac{ln(2)}{64})$. Since the output compensation function of
$sinh(x)$ also uses approximations of $sinh(r)$ and $cosh(r)$, we generate a
piecewise polynomial for $sinh(r)$ and a piecewise polynomial for $cosh(r)$
that can be used for the output compensation of both $sinh(x)$ and $cosh(x)$.
Using RLibm-32, we created a piecewise polynomial with $2^{6}$ polynomials of
degree 5 for $sinh(r)$ and a piecewise polynomial with $2^{6}$ polynomials of
degree 4 for $cosh(r)$.

##  Appendix C 32-bit Posit32 Functions

In this section, we explain the 32-bit posit (posit32) functions in RLibm-32.
We describe the special cases and any specific details on the range reduction
technique used for each function.

###  C.1. $ln(x)$

The elementary function $ln(x)$ is defined over the input domain $(0,\infty)$.
There are two classes of special case inputs of $ln(x)$ for posit32:

$ln(x)=\begin{cases}NaR&\text{if }x\leq 0\\\ NaR&\text{if }x=NaR\end{cases}$

We use the range reduction technique described in Appendix A.1 to decompose
$x$ into $n$, $F$, and $r$. In the case of $ln(x)$, the output compensation
can be mathematically computed as,

$ln(x)=ln(1+r)+ln(F)+n\times ln(2)$

To compute $n\times ln(2)$ accurately, we store the value of $ln(2)$ in two
double values, $t_{h}$ and $t_{l}$, such that $t_{h}+t_{l}$ stores the
correctly rounded value of $ln(2)$ with 106 precision bits (double type has 53
precision bits). The value $t_{h}$ stores the higher 53 precision bits and
$t_{l}$ stores the lower 53 precision bits. We evaluate the output
compensation function in double with the following order,

$ln(x)=((ln(1+r)+n\times t_{l})+ln(F))+n\times t_{h}$

This reduces the amount of numerical error for the inputs $x$ where the
magnitude of $n$ is large. After range reduction, there were roughly 115
million reduced inputs. The reduced inputs are in the range of
$r\in[0,\frac{1}{128})$. We created an approximation of $ln(1+r)$ using
RLibm-32 to generate a piecewise polynomial with $2^{10}$ polynomials of
degree 4.

###  C.2. $log_{2}(x)$

The elementary function $log_{2}(x)$ is defined over the input domain
$(0,\infty)$. There are two classes of special case inputs of $log_{2}(x)$ for
posit32:

$log_{2}(x)=\begin{cases}NaR&\text{if }x\leq 0\\\ NaR&\text{if
}x=NaR\end{cases}$

We use the range reduction technique described in Appendix A.1 to decompose
$x$ into $n$, $F$, and $r$. The output compensation function for $log_{2}(x)$
can be mathematically computed as,

$log_{2}(x)=log_{2}(1+r)+log_{2}(F)+n$

The output compensation function that we use for $log_{a}(x)$ experiences
cancellation error when $x=1-\epsilon$ for small values of $\epsilon$. In such
cases, $log_{2}(x)\approx 0$, $m=-1$, and $log_{2}(1+r)+log_{2}(F)\approx 1$.
The numerical error caused by cancellation error in the output compensation
function of $log_{2}(x)$ for posit32 poses a challenge in generating piecewise
polynomials of reasonable degree and size that satisfies all reduced input and
interval constraints.

Thus, we evaluate the output compensation function in double with the
following order of operations,

$log_{2}(x)=(log_{2}(1+r)+n)+log_{2}(F)$

This order subtracts two values with largest difference in magnitude, before
subtracting two values with similar magnitude. After range reduction, there
were roughly 115 million reduced inputs. The reduced inputs are in the range
of $r\in[0,\frac{1}{128})$. We created an approximation of $log_{2}(1+r)$
using RLibm-32 to generate a piecewise polynomial with $2^{8}$ polynomials of
degree 4. Comparatively, when we generated a piecewise polynomial for the
output compensation function that evaluates in the following order,

$log_{2}(x)=(log_{2}(1+r)+log_{2}(F))+n$

RLibm-32 generated a piecewise polynomial with $2^{11}$ polynomials of degree
4. In all other cases of $log_{a}(x)$ for float or posit32, adding
$log_{a}(1+r)+n\times log_{a}(2)$, then adding $log_{a}(F)$ at the end did not
produce piecewise polynomial with smaller number of polynomials.

###  C.3. $log_{10}(x)$

The elementary function $log_{10}(x)$ is defined over the input domain
$(0,\infty)$. There are two classes of special case inputs of $log_{10}(x)$
for posit32:

$log_{10}(x)=\begin{cases}NaR&\text{if }x\leq 0\\\ NaR&\text{if
}x=NaR\end{cases}$

We use the range reduction technique described in Appendix A.1 to decompose
$x$ into $n$, $F$, and $r$. The output compensation function of $log_{10}(x)$
can be mathematically computed as,

$log_{10}(x)=log_{10}(1+r)+log_{10}(F)+n\times log_{10}(2)$

To compute $n\times log_{10}(2)$ accurately, we store $log_{10}(2)$ in two
double values, $t_{h}$ and $t_{l}$. The sum $t_{h}+t_{l}$, if evaluated in
real numbers, is the correctly rounded value of $log_{10}(2)$ with 106
precision bits. The value $t_{h}$ stores the higher 53 precision bits and
$t_{l}$ stores the lower 53 precision bits. We evaluate the output
compensation function in double in the following order,

$log_{10}(x)=((log_{10}(1+r)+n\times t_{l})+log_{10}(F))+n\times t_{h}$

After range reduction, there were roughly 115 million reduced inputs. The
reduced inputs are in the range of $r\in[0,\frac{1}{128})$. We created an
approximation of $log_{10}(1+r)$ using RLibm-32 to generate a piecewise
polynomial with $2^{12}$ polynomials of degree 4.

###  C.4. $e^{x}$

The $e^{x}$ function is defined over the input domain $(-\infty,\infty)$.
There are 4 classes of special case inputs:

$e^{x}=\begin{cases}NaR&\text{if }x=NaR\\\ 2^{-120}&\text{if
}x\leq-81.7\dots\\\ 1.0&\text{if }-1.86\dots\times 10^{-9}\leq x\leq
3.72\dots\times 10^{-9}\\\ 2^{120}&\text{if }x\geq 81.7\dots\end{cases}$

We use a range reduction technique similar to the technique described in
Appendix A.2 to perform range reduction and output compensation. More
specifically, we split $x$ into $128$ segments instead of $64$ segments:

$x=nln(2)+\frac{j}{128}ln(2)+r$

where $n$ is an integer, $j$ is a value in a set $\\{0,1,2,\dots,127\\}$ and
$|r|\leq\frac{ln(2)}{128}$.

Then, the output compensation formula is adjusted accordingly:

$\displaystyle e^{x}=e^{nln(2)+\frac{j}{128}ln(2)+r}$
$\displaystyle=e^{nln(2)}\times e^{\frac{j}{64}ln(2)}\times e^{r}$
$\displaystyle=2^{n}\times 2^{\frac{j}{128}}\times e^{r}$

We pre-compute and store the value of $2^{\frac{j}{128}}$ in a table (_i.e._
128 values in total) and approximate $e^{r}$ for the input domain of
$r\in[-\frac{ln(2)}{128},\frac{ln(2)}{128}]$.

After range reduction, there are roughly 3.5 billion reduced inputs in the
domain $r\in[-\frac{ln(2)}{128},\frac{ln(2)}{128}]$. We created an
approximation function for negative values of reduced inputs $r$ and an
approximation function for positive values of reduced inputs to efficiently
split the reduced input domain. Using RLibm-32, we generated a piecewise
polynomial with $2^{12}$ polynomials of degree 3 for negative reduced inputs
and a piecewise polynomial with $2^{12}$ polynomials of degree 3 for positive
reduced inputs.

###  C.5. $2^{x}$

The $2^{x}$ function is defined over the input domain $(-\infty,\infty)$.
There are 4 classes of special case inputs:

$2^{x}=\begin{cases}NaR&\text{if }x=NaR\\\ 2^{-120}&\text{if
}x\leq-118\dots\\\ 1.0&\text{if }-2.68\dots\times 10^{-9}\leq x\leq
5.37\dots\times 10^{-9}\\\ 2^{120}&\text{if }x\geq 118\dots\end{cases}$

We use the range reduction technique described in Appendix A.2 to perform
range reduction and output compensation. After range reduction, there are
roughly 790 million reduced inputs in the domain
$r\in[-\frac{1}{64},\frac{1}{64}]$. We created an approximation function for
negative values of reduced inputs $r$ and an approximation function for
positive values of reduced inputs to efficiently split the reduced input
domain. Using RLibm-32, we generated a piecewise polynomial with $2^{10}$
polynomials of degree 3 for negative reduced inputs and a piecewise polynomial
with $2^{12}$ polynomials of degree 3 for positive reduced inputs.

###  C.6. $10^{x}$

The $10^{x}$ function is defined over the input domain $(-\infty,\infty)$.
There are 4 classes of special case inputs:

$10^{x}=\begin{cases}NaR&\text{if }x=NaR\\\ 2^{-120}&\text{if
}x\leq-35.5\dots\\\ 1.0&\text{if }-8.08\dots\times 10^{-10}\leq x\leq
1.61\dots\times 10^{-9}\\\ 2^{120}&\text{if }x\geq 35.5\dots\end{cases}$

We use the range reduction technique described in Appendix A.2 to perform
range reduction and output compensation. After range reduction, there are
roughly 3.4 billion reduced inputs in the domain
$r\in[-\frac{log_{10}(2)}{64},\frac{log_{10}(2)}{64}]$. We created an
approximation function for negative values of reduced inputs $r$ and an
approximation function for positive values of reduced inputs to efficiently
split the reduced input domain. Using RLibm-32, we generated a piecewise
polynomial with $2^{13}$ polynomials of degree 3 for negative reduced inputs
and a piecewise polynomial with $2^{13}$ polynomials of degree 3 for positive
reduced inputs.

###  C.7. $Sinh(x)$

The $sinh(x)$ function is defined over the input domain $(-\infty,\infty)$.
There are four classes of special cases:

$\text{sinh(x)}=\begin{cases}-2^{120}&\text{if }x<=-82.4\dots\\\ x&\text{if
}|x|\leq 2.79\dots\times 10^{-4}\\\ 2^{120}&\text{if }x>=82.4\dots\\\
NaR&\text{if }x=NaR\end{cases}$

We use the range reduction technique described in Appendix A.5 and evaluate
the range reduction and output compensation functions in double. After range
reduction, there were roughly 1.6 billion reduced inputs in the domain
$r\in[0,\frac{ln(2)}{64})$. Because the output compensation function for
$sinh(x)$ uses approximations of $sinh(r)$ and $cosh(r)$, we generated a
piecewise polynomial for $sinh(r)$ and a piecewise polynomial for $cosh(r)$.
Using RLibm-32, we created a piecewise polynomial with $2^{14}$ polynomials of
degree 5 for $sinh(r)$ and a piecewise polynomial with $2^{14}$ polynomials of
degree 4 for $cosh(r)$.

###  C.8. $Cosh(x)$

The $cosh(x)$ functions is defined over the input domain $(-\infty,\infty)$.
There are three classes of special cases:

$\text{cosh(x)}=\begin{cases}1.0&\text{if }|x|\leq 8.63\dots\times 10^{-5}\\\
2^{120}&\text{if }|x|>=82.4\dots\\\ NaR&\text{if }x=NaR\end{cases}$

We use the range reduction technique described in Appendix A.6 and evaluate
the range reduction and output compensation functions in double. After range
reduction, there were roughly 1.7 billion reduced inputs in the domain. We
approximate $sinh(r)$ and $cosh(r)$ for the reduced inputs in the domain
$r\in[0,\frac{ln(2)}{64})$. Since the output compensation function of
$cosh(x)$ uses approximations of $sinh(r)$ and $cosh(r)$, we generate a
piecewise polynomial for $sinh(r)$ and a piecewise polynomial for $cosh(r)$.
Using RLibm-32, we created a piecewise polynomial with $2^{14}$ polynomials of
degree 3 for $sinh(r)$ and a piecewise polynomial with $2^{12}$ polynomials of
degree 6 for $cosh(r)$.