RLIBM-ALL: A Novel Polynomial Approximation Method to Produce Correctly Rounded Results for Multiple Representations and Rounding Modes

Mainstream math libraries for floating point (FP) do not produce correctly rounded results for all inputs. In contrast, CR-LIBM and RLIBM provide correctly rounded implementations for a specific FP representation with one rounding mode. Using such libraries for a representation with a new rounding mode or with different precision will result in wrong results due to double rounding. This paper proposes a novel method to generate a single polynomial approximation that produces correctly rounded results for all inputs for multiple rounding modes and multiple precision configurations. To generate a correctly rounded library for $n$-bits, our key idea is to generate such a polynomial approximation for a representation with $n+2$-bits using the \emph{round-to-odd} mode. We prove that the resulting polynomial approximation will produce correctly rounded results for all five rounding modes in the standard and for multiple representations with $k$-bits such that $|E| +1<k \leq n$, where $|E|$ is the number of exponent bits in the representation. Building on our prior work in the RLIBM project, we also approximate the correctly rounded result when we generate the library with $n+2$-bits using the round-to-odd mode. We also generate polynomial approximations by structuring it as a linear programming problem but propose enhancements to polynomial generation to handle the round-to-odd mode. Our prototype is the first 32-bit float library that produces correctly rounded results with all rounding modes in the IEEE standard for all inputs with a single polynomial approximation. It also produces correctly rounded results for any FP configuration ranging from 10-bits to 32-bits while also being faster than mainstream libraries.

INTRODUCTION
The floating point (FP) representation is widely used to approximate real numbers. The two main attributes of the FP representation are its dynamic range (i.e., the range of values that can be represented) and precision (i.e., the accuracy of each value represented). As some real numbers cannot be accurately represented in the FP representation, they need to be rounded to the nearest result according to the rounding mode. Further, FP performance is important in various domains ranging from scientific computing to machine learning. Hence, modern accelerators, processors, and systems have explored new variants of the FP representation that vary the dynamic range and the precision. Intel's Bfloat16 [Kalamkar et al. 2019] and FlexPoint [Köster et al. 2017], Nvidia's TensorFloat32 [NVIDIA 2020], Microsoft's MSFP [Rouhani et al. 2020], Facebook's Log Number System [Johnson 2018], and Posits [Gustafson and Yonemoto 2017] are examples of such recent FP variants. All these representations need math libraries that provide approximations for various elementary functions (e.g., ( ), , . . . ) [Muller 2005].
An output of an elementary function is a correctly rounded result if it matches the result that is computed with infinite precision and then rounded to the target FP representation. Correctly rounded elementary functions improve the portability and reproducibility of applications. Unfortunately, mainstream FP libraries do not produce correct results for all inputs. Libraries like CR-LIBM [Daramy et al. 2003;Daramy-Loirat et al. 2006] and RLibm [Lim et al. 2020Lim and Nagarakatte 2021a,b] provide correctly rounded implementations for some FP representations. CR-LIBM provides correctly rounded elementary functions for the double type with a given rounding The next task is to generate polynomial approximations using the round-to-odd mode for a representation with + 2-bits (i.e., T +2 ). We extend our prior work in the RLibm project to generate polynomial approximations with the round-to-odd mode. Specifically, we approximate the correctly rounded result of an elementary function ( ) with the round-to-odd mode rather than the real value of ( ). Subsequently, we identify an interval of real values around the correctly rounded result with the round-to-odd mode such that any value in the interval rounds to the correctly rounded result, which we call the odd interval. We show that any real value in the odd interval that rounds to the round-to-odd result in T +2 will subsequently round to the correctly rounded result for any representation T .
One challenge in generating polynomial approximations using the odd intervals is the presence of singleton values with the round-to-odd mode. When the real value is exactly representable in T +2 and the bit-string of that value in T +2 is even when interpreted as an unsigned integer, the odd interval for that value will be a singleton element. We use the mathematical properties of the elementary function to identify such singletons. Concretely, we identify the inputs of an elementary function that has rational outputs (because the round-to-odd result is exactly representable in T +2 ). We subsequently develop efficient table-lookups to produce round-to-odd results for such inputs (Section 4.3).
Once we identify the odd interval for every input and the singleton odd intervals among them, we can subsequently use the RLibm approach to generate polynomial approximations. Each nonsingleton interval imposes a constraint on the output of the generated polynomial for a given input. Similar to our prior work in the RLibm project [Lim 2021], we structure the problem of generating a polynomial approximation that satisfies the odd interval for each input as a linear programming problem. We use efficient range reduction and output compensation functions while accounting for numerical errors with them. To account for numerical errors, we further constrain the odd intervals given to the LP formulation. We also employ counterexample guided polynomial generation and generate piecewise polynomials for efficiency.
Our prototype, RLibm-All, is open-source and publicly available [Lim and Nagarakatte 2021c]. It contains the polynomial generator for ten elementary functions and includes efficient polynomial approximations for them. The resulting polynomial generated by our prototype for the FP representation with 34-bits (i.e., T 34 ) using the round-to-odd mode produces correctly rounded results for all FP representations ranging from 10-bits (i.e., T 10 ) to 32-bits (i.e., T 32 ) and for all rounding modes. It includes bfloat16 and tensorfloat32 representations. RLibm-All is the first math library that produces correctly rounded results for 32-bit floats for all rounding modes with a single polynomial approximation. Our implementations are faster than mainstream libraries for 32-bit floats while producing correctly rounded results for all inputs and for all rounding modes.

BACKGROUND
We provide background on the FP representation, the process of rounding, and the various rounding modes in the IEEE-754 standard. As we extend our RLibm approach, we also provide a brief background on generating polynomial approximations with it. When a real value R is not exactly representable in T, then R is rounded to one of the two adjacent values sm , l 2 T depending on the rounding mode used. We show the range of R that rounds to sm (blue box) and l (green box) if we use (a) rne mode when the bit-string of sm is even, (b) rne mode when the bit-string of l is even, (c) rna mode when x < 0, (d) rna mode when x > 0, (e) rnz mode when x < 0, (f) rnz mode when x > 0, (g) rnp mode, and (h) rnn mode.
sign bit is 0, then the value is positive. Otherwise, it is negative. The value represented by the FP bit-string can be classied into three classes: normal values, denormal values, and special values. The value represented by the FP bit-string is a normal value if the bit-string E is neither all ones nor all zeros (i.e., 0 < E < 2 |E | 1). The normal value represented by this bit-string is where bias is 2 |E | 1 1. If the exponent eld E is all zeros (i.e., E = 0), then FP value is a denormal value. Denormal values are used to represent values close to zero. The number represented by this denormal value is ( F 2 |F | ) ⇥ 2 1 bias . When the exponent eld is all ones (i.e., E = 2 |E | 1), the FP bit-strings represent special values. If F = 0, then the bit-string represents ±1 depending on the sign and in all other cases, it represents not-a-number (NaN).

Rounding a Real Number to the FP representation
Any FP representation can represent a nite number of real values. Hence, many real values ( R ) cannot be exactly represented. It is rounded to either the largest FP value smaller than R ( sm ) or the smallest FP value larger than R ( l ).
The rounding mode, which we represent by rm, species whether R rounds to sm or l . We denote the operation of rounding R using a rounding mode rm to a value in the representation FP with RN F n, |E | ,rm ( R ). The IEEE-754 standard species ve dierent rounding modes: round to Proc. ACM Program. Lang., Vol. 1, No. POPL, Article 1. Publication date: January 2022. (a) Round-to-nearest-tie-goes-to-even when v sm is even (b) Round-to-nearest-tie-goes-to-even when v  2. When a real value R is not exactly representable in T, then R is rounded to one o values sm , l 2 T depending on the rounding mode used. We show the range of R th (blue box) and l (green box) if we use (a) rne mode when the bit-string of sm is even, ( the bit-string of l is even, (c) rna mode when x < 0, (d) rna mode when x > 0, (e) rnz mo rnz mode when x > 0, (g) rnp mode, and (h) rnn mode.
sign bit is 0, then the value is positive. Otherwise, it is negative. The value repre bit-string can be classied into three classes: normal values, denormal values, an The value represented by the FP bit-string is a normal value if the bit-string ones nor all zeros (i.e., 0 < E < 2 |E | 1). The normal value represented by is not exactly representable in T, then R is rounded to one of the two adjacent ng on the rounding mode used. We show the range of R that rounds to sm ) if we use (a) rne mode when the bit-string of sm is even, (b) rne mode when ) rna mode when x < 0, (d) rna mode when x > 0, (e) rnz mode when x < 0, (f) mode, and (h) rnn mode. e is positive. Otherwise, it is negative. The value represented by the FP into three classes: normal values, denormal values, and special values. by the FP bit-string is a normal value if the bit-string E is neither all < E < 2 |E | 1). The normal value represented by this bit-string is bias is 2 |E | 1 1. If the exponent eld E is all zeros (i.e., E = 0), then lue. Denormal values are used to represent values close to zero. The s denormal value is ( F 2 |F | ) ⇥ 2 1 bias . When the exponent eld is all ones t-strings represent special values. If F = 0, then the bit-string represents and in all other cases, it represents not-a-number (NaN). e the 16-bit half type (F 16,5 ), the 32-bit oat type (F 32,8 ), and the 64-bit d these types, recent extensions have increased the dynamic range and/or ontext of machine learning. The new types include boat16 (F 16,8 ) 8 ) [NVIDIA 2020].
umber to the FP representation represent a nite number of real values. Hence, many real values ( R ) ted. It is rounded to either the largest FP value smaller than R ( sm ) or r than R ( l ).
R } ich we represent by rm, species whether R rounds to sm or l . We unding R using a rounding mode rm to a value in the representation he IEEE-754 standard species ve dierent rounding modes: round to 1, No. POPL, Article 1. Publication date: January 2022. Anon.
v sm (even) v lg (odd) tie-point (a) Round-to-nearest-tie-goes-to-even when v sm is even v sm (odd) All values in this interval rounds to v lg (even) tie-point All values in this interval rounds to (b) Round-to-nearest-tie-goes-to-even when v lg is even When a real value R is not exactly representable in T, then R is rounded to one of the two adjacent values sm , l 2 T depending on the rounding mode used. We show the range of R that rounds to sm (blue box) and l (green box) if we use (a) rne mode when the bit-string of sm is even, (b) rne mode when the bit-string of l is even, (c) rna mode when x < 0, (d) rna mode when x > 0, (e) rnz mode when x < 0, (f) rnz mode when x > 0, (g) rnp mode, and (h) rnn mode.
sign bit is 0, then the value is positive. Otherwise, it is negative. The value represented by the FP bit-string can be classied into three classes: normal values, denormal values, and special values. The value represented by the FP bit-string is a normal value if the bit-string E is neither all ones nor all zeros (i.e., 0 < E < 2 |E | 1). The normal value represented by this bit-string is where bias is 2 |E | 1 1. If the exponent eld E is all zeros (i.e., E = 0), then FP value is a denormal value. Denormal values are used to represent values close to zero. The number represented by this denormal value is ( F 2 |F | ) ⇥ 2 1 bias . When the exponent eld is all ones (i.e., E = 2 |E | 1), the FP bit-strings represent special values. If F = 0, then the bit-string represents ±1 depending on the sign and in all other cases, it represents not-a-number (NaN).

Rounding a Real Number to the FP representation
Any FP representation can represent a nite number of real values. Hence, many real values ( R ) cannot be exactly represented. It is rounded to either the largest FP value smaller than R ( sm ) or the smallest FP value larger than R ( l ).
The rounding mode, which we represent by rm, species whether R rounds to sm or l . We denote the operation of rounding R using a rounding mode rm to a value in the representation

The Floating Point Representation
The floating point representation, which is specified by the IEEE-754 standard, is parameterized by the total number of bits and the number of bits for the exponent | |, which we represent as F , | | . The total number of bits and the number of bits for the exponent determine the dynamic range and precision of the representation. The FP bit-string consists of a sign bit, | | bits to represent the exponent, and − 1 − | | bits to represent the mantissa ( ). Figure 1 shows the bit-string for a standard 32-bit float and the custom 5-bit and 4-bit FP representations. If the sign bit is 0, then the value is positive. Otherwise, it is negative. The value represented by the FP bit-string can be classified into three classes: normal values, denormal values, and special values.
The value represented by the FP bit-string is a normal value if the bit-string is neither all ones nor all zeros (i.e., 0 < < 2 | | − 1). The normal value represented by this bit-string is (1 + 2 | | ) × 2 − , where bias is 2 | |−1 − 1. If the exponent field is all zeros (i.e., = 0), then the FP value is a denormal value. Denormal values are used to represent values close to zero. The real number represented by this denormal value is ( 2 | | ) × 2 1− . When the exponent field is all ones (i.e., = 2 | | − 1), the FP bit-strings represent special values. If = 0, then the bit-string represents ±∞ depending on the sign and in all other cases, it represents not-a-number (NaN).
The default FP types are the 16-bit half type (F 16,5 ), the 32-bit float type (F 32,8 ), and the 64-bit double type (F 64,11 ). Beyond these types, recent extensions have increased the dynamic range and/or precision especially in the context of machine learning. The new types include bfloat16 (F 16,8 ) [Tagliavini et al. 2018] and tensorfloat32 (F 19,8 ) [NVIDIA 2020].

Rounding a Real Number to the FP Representation
Any FP representation can represent a finite number of real values. Hence, many real values ( R ) cannot be exactly represented. It is rounded to either the largest FP value smaller than R ( ) or the smallest FP value larger than R ( ).
rounding mode, which we represent by , specifies whether R rounds to or . We denote the operation of rounding R using a rounding mode to a value in the representation The IEEE-754 standard specifies five different rounding modes: roundto-nearest-ties-to-even ( ), round-to-nearest-ties-to-away ( ), round-towards-zero ( ), roundtowards-positive-infinity ( ), and round-towards-negative-infinity ( ). The standard mandates correct rounding for primitives operations (i.e., +, −, * , /).
The round-to-nearest-ties-to-even ( ) mode. This rounding mode rounds R to a FP value that is closer to R among and . If | − R | < | − R |, then R rounds to . If | − R | > | − R |, then R rounds to . If R is exactly in the middle of and (i.e., R = + 2 ), R is rounded to a value whose bit-string is even when interpreted as an unsigned integer. The mode is the most commonly used rounding mode. Figure 2(a) and (b) illustrate rounding with the mode depending on whether is even or odd, respectively.
The round-to-nearest-ties-to-away ( ) mode. This rounding mode also rounds R to a FP value that is closer to R among and . If | − R | < | − R |, then R rounds to . If | − R | > | − R |, then R rounds to . When R is exactly in the middle of and , R is rounded to a value that is farther away from 0. Specifically, R rounds to Figure 2(c) and (d) illustrate rounding with mode depending on whether R < 0 or R > 0, respectively.
The round-towards-zero ( ) mode. In this mode, R is rounded to a value that is closer to 0. Here, R is rounded to if R > 0 and R is rounded to if R < 0. The mode is equivalent to truncating the fraction bits of R that cannot fit within the mantissa bits of the representation. Figure 2(e) and (f) illustrate the mode depending on whether R < 0 or R > 0, respectively.
The round-towards-positive-infinity ( ) mode. This mode always rounds R to the larger value , which is the value that is closer to +∞. This mode is also known as rounding up. Figure 2(g) pictorially shows the mode.
The round-towards-negative-infinity ( ) mode. The round-towards-negative-infinity ( ) mode always rounds R to the smaller value (i.e., a value that is closer to −∞). This mode is also known as rounding down. Figure 2(h) demonstrates the mode.

A Systematic Method for Rounding
We describe a systematic procedure for rounding a real number, which will be useful later for understanding our proofs. As we described above, we need to identify the two values and in the F , | | representation that are adjacent to R and then decide between or . We will identify four pieces of information ( , − , , ) from the real value R that will be sufficient to correctly round according to the various rounding modes. We call them rounding components. The first component, represents the sign (-1 or 1) and identifies whether R is positive or negative. The smaller of or in magnitude is represented by − . The components and encode information about whether R is in the middle, closer to , or closer to . To identify the rounding components, we represent R in the FP representation with an infinite number of mantissa bits while having the same number of exponents bits as F , | | . We call this representation extended infinite precision representation (i.e., F ∞, | | ). Effectively, this extended precision representation is similar to the F , | | representation but has a large number of bits for the mantissa. When R is larger than the dynamic range of F , | | , we represent it with the largest representable value in F ∞, | | (i.e., the exponents bits correspond to the largest normal value in F , | | and all the mantissa bits are ones). We cannot use ∞ to represent | R | because we need to make a clear distinction between a real number and ∞. The round-towards-zero mode never rounds a real value to ∞. Similarly, the round-towards-positive-infinity mode does not round negative real values to −∞ and the round-towards-negative-infinity mode does not round positive real values to ∞.
Here, 1 is the sign bit and bits 2 . . . | |+1 represent the exponent bits. The rest of the bits starting from | |+2 . . . represent the mantissa.
Identify rounding components. To identify and that R can round to, we identify two positive values − and + adjacent to | R | (i.e., the magnitude of R ) in F , | | . Here, − represents the largest value that is smaller than or equal to | R | and + represents the smallest value larger than | R |. To identify − , we truncate | R | to bits, − = 0 2 3 4 . . . −1 Note that the sign bit is 0 because we are just considering the magnitude for − . We call − the truncated value, which is a rounding component. Then, the succeeding value of − in F , | | is + , which is obtained by adding 1 to − . We maintain the invariant: − ≤ | R | < + . In the context of rounding R to F , | | , − and + satisfy the following property, Once we identify , − , and + , we can compute and as follows.
Rounding bit. To determine the rounding decision for the and mode, we must determine whether R is closer to × − , closer to × + , or exactly in the middle of the two values. We extract the ( + 1) ℎ -bit from our extended precision representation of | R | , which we call as the rounding bit (i.e., ). The rounding bit describes whether | R | is closer to − than + . If the rounding bit is 0, then | R | is closer to − (i.e., − ≤ | R | < − + + 2 ). If the rounding bit is 1, then | R | is at the middle or close to + . Figure 3 illustrates the range of real values where the rounding bit is 0 or 1.
Sticky bit. While the rounding bit tells us whether | R | is closer to − , it does not tell us whether | R | is exactly equal to − or is exactly in the middle of − and + (i.e., | R | = − + + 2 ). When we look at the bit-string | R | , | R | is equal to − when the ( + 1) ℎ -bit (i.e., ) is 0 and the remaining bits from the ( + 2) ℎ -bit are all zeros. If = 0 and any bit afterwards is 1, then | R | is not equal to − . Similarly, | R | is exactly in the middle of − and + when the ( + 1) ℎ -bit is 1 and the remaining bits from the ( + 2) ℎ -bit are all 0's in | R | . In both these cases, we need to determine if all the 757:7 . We illustrate rounding when R is positive or negative and when − is even or odd. When a interval is colored green, all real values in the interval round to the FP value colored green. Similarly, all real values in the interval colored blue will round to the FP value colored blue.
bits starting from +2 are zeros. We define the sticky bit as the bitwise OR of all bits starting from the ( + 2) ℎ -bit in the extended precision representation.
Using the rounding components ( , − , , ), we can identify the relationship between | R | and the nearest FP values for any rounding mode in the standard.
We can compute + from the rounding components. Figure 3 pictorially shows the rounding bit and sticky bit for various real values between − and + .
Rounding to various modes with the rounding components. Figure 4 shows rounding R using the rounding components for the mode. Similarly, Figure 5 illustrates rounding with the rounding components for the other four rounding modes.

The RLibm Approach
We provide a brief background on our RLibm project as we build on top of it in this paper. In the RLibm project, we make a case for approximating the correctly rounded result rather than the real value of an elementary function [Lim et al. 2020Lim and Nagarakatte 2021a,b]. When we approximate the correctly rounded result, there is an interval of real values around the correctly rounded result for each input such that producing any value in the interval produces the correct result. This interval is further constrained to account for numerical errors that can occur with polynomial evaluation, range reduction, and output compensation. This interval can be used to generate polynomial approximations. It represents the maximum amount of freedom available to produce the correct result. Figure 6 illustrates our RLibm approach.
The RLibm approach consists of four steps. The first step is to use an oracle to compute the correctly rounded result of an elementary function ( ) for each input ∈ T, where T is the target representation. The second step is to identify an interval [ , ℎ] around the correctly rounded result such that any value in [ , ℎ] rounds to the correctly rounded result in T, which is known as the rounding interval. Since polynomial evaluation, range reduction, and output compensation happen in representation with higher precision H, the rounding intervals are also in H. The third step is to employ range reduction to transform input to ′ . The generated polynomial will approximate the result for ′ . Subsequently, we use output compensation to produce the final correctly rounded output for . Both range reduction and output compensation happen in H and can experience numerical errors. These numerical errors should not affect the generation of correctly rounded results. Hence, it is necessary to deduce intervals for the reduced domain so that the polynomial evaluation over the reduced input produces the correct results for the original inputs. Given and its rounding interval [ , ℎ], reduced input ′ is computed with range reduction. The next task before polynomial generation is to identify the reduced rounding interval for ( ′ ) that when used with output compensation produces the correctly rounded result. We use the inverse of the output compensation function to identify the reduced interval [ ′ , ℎ ′ ].
The last step is to synthesize a polynomial of a degree using an arbitrary precision linear programming (LP) solver that satisfies the constraints (i.e., ′ ≤ ( ′ ) ≤ ℎ ′ ) when given a set of inputs ′ . Approximating the correctly rounded result with the RLibm approach provides more freedom in generating polynomials. Hence, the resulting RLibm libraries are more efficient compared to mainstream libraries.
rounds to rounding interval of v 2 Fig. 6. Illustration of the RLibm approach. (a) The values 1 , 2 , and 3 are representable values in representation T. The real value of ( ) for a given input cannot be exactly represented in T and it is rounded to 2 . The RLibm approach identifies the rounding interval of 2 (shown in gray box). (b) Polynomial generation using the rounding interval (i.e., [ , ℎ]) for each input with an LP formulation.

AN ILLUSTRATIVE EXAMPLE OF OUR APPROACH
We describe our entire approach with an end-to-end example for creating a polynomial approximation for ( ) that produces correctly rounded results for a 5-bit FP representation with 2 exponent bits (FP5) and a 4-bit FP with 2 exponent bits (FP4) for all standard rounding modes (i.e., , , , , and ). Figure 1(b) and Figure 1(c) show the bit-string of FP5 and FP4, respectively. Although we illustrate our approach with FP5 and FP4 for ease of exposition, it is beneficial in practice to create table-lookups for FP5 and FP4 because there are only 32 and 16 distinct bit-patterns, respectively.
The ( ) function is defined over the input domain (0, ∞). The result of ( ) is NaN when < 0 or when x is NaN. The result is ∞ when the input is ∞ and −∞ when the input is 0. There are only 11 non-special case inputs, which range from 0.25 to 3.5 in FP5. Similarly, there are only 5 non-special case inputs in FP4. Now, our goal is to generate a single polynomial approximation that produces correctly rounded results for both FP5 and FP4 with all five rounding modes.
To accomplish this goal, we will generate a polynomial approximation that produces correctly rounded results for a 7-bit FP representation (FP7) with the round-to-odd mode. Here, FP7 has exactly the same number of exponent bits as FP5 and FP4 (i.e., 2 exponent bits). Effectively, FP7 has 2 additional fraction bits when compared to FP5. Every value that is representable in FP5 and FP4 is also representable in FP7. While rounding with the round-to-odd mode, if the polynomial approximation produces a value that is exactly representable in FP7, then it is unchanged. Otherwise, the result of the polynomial approximation (which is implemented in double precision) is rounded to the nearest FP7 value whose bit-string is odd (i.e., the last bit is a 1). When this FP7 round-to-odd result is rounded to a value in FP5 or FP4 according to any of the five standard rounding modes, it produces the correct result for them.
Why does a correctly rounded result with the round-to-odd mode for FP7 work with FP5/FP4? As the number of exponent bits is identical in FP7, FP5, and FP4, every value that is representable in FP5 and FP4 is also representable in FP7. Let us consider an input 1.5. We want to produce correctly rounded results of (1.5) for all the rounding modes with FP5 and FP4. The first row of Figure 7 shows the real result (i.e., a star) and the correctly rounded result of FP5 with the mode. If we want to generate the correctly rounded result for (1.5) in FP5 with the mode using polynomial approximations, there is an interval of real values around the correctly rounded result such that producing any value in that interval produces the correct result (shaded in gray). The subsequent rows show the correctly rounded result and the rounding intervals for other rounding modes of FP5 and FP4 for the same input 1.5. When compared to FP5, the rounding interval for FP4 will be larger because the distance between adjacent points is larger. Intuitively, a single polynomial can produce the correctly rounded result of (1.5) for both FP5 and FP4 with all rounding modes if it produces a value that lies in the common interval among all these modes and precision configurations (i.e., an intersection of the rounding intervals). The correctly rounded result of (1.5) for FP5 and FP4 with some subset of the rounding modes and their rounding intervals (gray box). The gray star represents the real value of (1.5). Values that are representable in FP7, FP5, and FP4 are shown with rhombus, circle, and square, respectively. Solid shapes represents the correctly rounded result for the chosen representation and rounding mode. The last row shows the odd interval to produce the correctly rounded result of (1.5) in FP7 with the round-to-odd mode. The odd interval is a subset of the intersection of the rounding intervals of these configurations.
We show that computing the correctly rounded result of (1.5) with FP7 using the round-to-odd mode and identifying the interval around this result in FP7 is an effective way to compute the common interval described above. The last row of Figure 7 shows the correctly rounded result in FP7 with the round-to-odd mode and the interval to produce that value. The interval for the round-to-odd result in FP7 is smaller than the common interval among FP5 and FP4 with all the rounding modes because it works for many other representations beyond FP5 and FP4.
In FP7, there are three additional values between the two adjacent FP5 values. Hence, the roundto-odd result with FP7 preserves enough information to produce the correctly rounded result with FP5 and FP4 with any rounding mode. Our proofs in Section 5 show that this is a generic result for any representation with -bits.
Generating polynomial approximations. The first step is to identify rounding intervals for producing the correctly rounded result of FP7 with the round-to-odd mode, which we call the odd interval. In our approach, polynomial evaluation happens with double precision. Hence, we identify an interval of values in double precision such that any value in that interval rounds to the correctly rounded result in FP7 with the round-to-odd mode. When the correctly rounded result in FP7 with the round-to-odd mode is even (i.e., the bit-string is even when interpreted as an unsigned integer), the rounding interval is a singleton. For example, the odd interval for (1.0) is a singleton because the correctly rounded result is 0. If the correctly rounded result with the round-to-odd mode is not even, then we can identify the odd interval as follows. We identify the preceding value (l) and the succeeding (h) value corresponding to the correctly rounded result in FP7. Then, the open interval ( , ℎ) is the odd interval. Figure 8(b) shows the odd interval for each input (shaded in blue). Next, we need to create a polynomial approximation that produces a value in the odd interval for each input.
Creating polynomial approximations with singleton odd intervals is challenging because there is no freedom for the polynomial generator. For the ( ) function with FP7, there is only one input (i.e., 1.0) whose odd interval is a singleton. We treat it as a special case. In general, we use mathematical properties of the elementary function for larger data types to effectively handle such 757:11 1:14 Anon. that solves for the coecients of P(x). We create a system of linear inequalities that specify the constraints of P(x) (Figure 10(a)) and use an LP solver to solve for the coecients of P(x) that satises all the constraints. If the LP solver produces a solution, then the resulting polynomial P(x) created using the coecients (Figure 10(b)) satises all the constraints specied by the generic intervals. Figure 10(c) pictorially shows that the polynomial we generate does indeed produce a value in the generic interval for each input. Thus, the result of the polynomial will round to the correctly rounded result of f (x) when rounded to FP5 or FP4 representation using any rounding mode.
Proc. ACM Program. Lang., Vol. 1, No. POPL, Article 1. Publication date: January 2022. that solves for the coecients of P(x). We create a system of linear inequalities that specify the constraints of P(x) (Figure 10(a)) and use an LP solver to solve for the coecients of P(x) that satises all the constraints. If the LP solver produces a solution, then the resulting polynomial P(x) created using the coecients (Figure 10(b)) satises all the constraints specied by the generic intervals. Figure 10(c) pictorially shows that the polynomial we generate does indeed produce a value in the generic interval for each input. Thus, the result of the polynomial will round to the correctly rounded result of f (x) when rounded to FP5 or FP4 representation using any rounding mode.  singleton odd intervals (Section 4.3). Figure 8(b) shows the remaining inputs and their non-singleton odd intervals. The next step is to generate a polynomial that produces a value in the odd interval for all inputs. We show the constraints imposed by the odd interval on the output of the polynomial in Figure 8(a). As there only 10 non-special case inputs, we encode them as a system of linear inequalities similar to our prior work in the RLibm project and use an LP solver to solve for the coefficients of a 4 ℎ degree polynomial ( ) (see Figure 8(a)). For larger representations, we employ sophisticated range reduction, counterexample guided polynomial generation, and generate piecewise polynomials [Lim and Nagarakatte 2021a]. Figure 8(b) pictorially shows the generated polynomial, which produces a value in the odd interval for each input. This polynomial will produce the correctly rounded result of ( ) when the result is rounded to FP5 and FP4 with any of the five rounding modes in the standard.

OUR APPROACH TO GENERATE A GENERIC POLYNOMIAL APPROXIMATION
Our goal is to generate a single polynomial approximation of an elementary function that produces correctly rounded results for all inputs for multiple precision and rounding configurations. Let T be a -bit FP representation (i.e., F , | | ). Let T be a representation where T has no more precision bits compared to T with the same number of exponent bits. Specifically, T = F , | | where | | + 1 < ≤ . Note that all values exactly representable in T are also exactly representable in T (i.e., T ⊆ T ). We define to be a rounding mode in the standard (i.e., ∈ { , , , , }). Our goal is to generate a polynomial approximation H ( ), which is implemented in representation H, of an elementary function ( ) that produces correctly rounded results for all inputs for any representation T and any rounding mode . Specifically, rounding the result of H ( ) to any representation T with the rounding mode must result in the same value as computing ( ) in real numbers and rounding the result to T with the rounding mode, for all inputs in T .
T , ( H ( )) = T , ( ( )) Main insight. To generate correct results for T , our key insight is to create a polynomial approximation that produces the correctly rounded results for T +2 with the round-to-odd mode. We prove that it produces the correctly rounded result for any representation T with all standard rounding modes when we round the round-to-odd result to the target representation (Section 5). Intuitively, our approach works because the round-to-odd result with T +2 maintains sufficient information w 0 w 1 w 2 w 3 w 4 All values in this region rounds to w 1 All values in this region rounds to w 3 (even) (even) (even) (odd) (odd) Fig. 9. The round-to-odd ( ) rounding mode. We show the rounding of R with the mode. Here, 0 , 1 , 2 , 3 , and 4 are values representable in representation T. If R is exactly representable in T, then R rounds to that value. Otherwise, R rounds to the nearest value in T that is odd. about the real value of an elementary function ( ) that is required for correct rounding for all representations T with any standard rounding mode.
To generate a polynomial approximation for T +2 with the round-to-odd mode, we use our RLibm approach. We approximate the correctly rounded result rather than the real value. We extend the RLibm approach to handle the round-to-odd mode. Specifically, we need to generate an interval of values around the correctly rounded round-to-odd result for each input, which we call the odd interval. If the generated polynomial produces a value in the odd interval for a particular input, then it produces the correct result for all representations T and for all rounding modes. One unique challenge that we address is the presence of singleton odd intervals, which happens when the correctly rounded result in T +2 is even. To scale to 32-bit floats (i.e., T +2 = 34-bit float), we employ range reduction, counterexample guided polynomial generation, and generate piecewise polynomials. Finally, we use a linear programming solver to solve for the coefficients of a polynomial given a system of linear constraints generated from the odd intervals.

Generating the Correctly Rounded Result for T +2 with the Round-to-Odd Mode
As we make a case for creating polynomial approximations for T +2 with the round-to-odd mode, we formally define it and describe rounding with the round-to-odd mode using the rounding components. We describe the properties of the correctly rounded result with the round-to-odd mode and provide intuition on why rounding the round-to-odd result in T +2 to any representation T produces the correct result.
The round-to-odd ( ) is a non-standard rounding mode that has been previously used to avoid double rounding issues while converting a binary FP number to a decimal FP number [Goldberg 1991] and while performing primitive operations Melquiond 2005, 2008]. Given a real value R , the round-to-odd mode rounds R as follows. If R is exactly representable as value in the target representation, then R rounds to . Otherwise, R rounds to the nearest odd value in the target representation. Figure 9 illustrates the round-to-odd rounding mode. Using the rounding components ( , − , , ) from Section 2.2, the round-to-odd mode can be defined as follows: where + is the adjacent value to − in T. Our contribution is to use the round-to-odd mode to generate correctly rounded elementary functions for multiple representations and multiple types. Specifically, we prove the following theorem in Section 5, which forms the foundation for our approach.  the exact result of f (x) for a given input x 2 T n , then Tn+2,rno (f (x))) = RN Tk,rm (f (x)) sult of RN Tn+2,rno (f (x)). Then, if our approximation ll inputs x 2 T n , then A(x) rounds to the correctly sentations T k with any rounding mode rm. Since any re also exactly representable in T n , it is sufficient for es y rno for all x 2 T n . Technically, the interval that , is a generic interval for x. However, generating A(x) (where the only option is to produce y rno ) for all x 2 T n to find a range of values near y rno that is guaranteed to lt of f (x) for all T k and rm.
Odd Interval of v rno ing double rounding by first rounding v R to T n+2 using rno = RN Tn+2,rno (v R ) no to T k using a rounding mode rm produces v rm , the v R directly to T k using the same rounding mode rm, The odd interval of v rno is a range e v in the interval rounds to v rno when rounded to T n+2 148 m produces v rm , the rounding mode rm, rval of v rno is a range hen rounded to T n+2 v 1 (odd) v 0 (even) w 0 w 1 w 2 w 3 w 4 mode rm. If we substitute v R with the exact result of f (x) for a given input x 2 T n , then it is guaranteed that, RN Tk,rm (RN Tn+2,rno (f (x))) = RN Tk,rm (f (x)) Let us denote y rno to be the result of RN Tn+2,rno (f (x)). Then, if our approximation function A(x) produces y rno for all inputs x 2 T n , then A(x) rounds to the correctly rounded result of f (x) for all representations T k with any rounding mode rm. Since any exactly representable value in T k are also exactly representable in T n , it is sufficient for us to make sure that A(x) produces y rno for all x 2 T n . Technically, the interval that contains a single value, [y rno , y rno ], is a generic interval for x. However, generating A(x) that produces a value in [y rno , y rno ] (where the only option is to produce y rno ) for all x 2 T n is infeasible. Thus, our next step is to find a range of values near y rno that is guaranteed to round to the correctly rounded result of f (x) for all T k and rm.

Formal Reasoning On the Odd Interval of v rno
Theorem 6.12 tells us that performing double rounding by first rounding v R to T n+2 using and then subsequently rounding v rno to T k using a rounding mode rm produces v rm , the value that we obtain when rounding v R directly to T k using the same rounding mode rm, unding v rno to T k using a rounding mode rm produces v rm , the n rounding v R directly to T k using the same rounding mode rm, v rm = RN Tk,rm (v R ) RN Tk,rm (v rno ) = v rm define the odd interval of v rno . The odd interval of v rno is a range re any value v in the interval rounds to v rno when rounded to T n+2 148 v 1 (odd) v 0 (even) w 0 w 1 w 2 w 3 w 4 mode rm. If we substitute v R with the exact result of f (x) for a given input x 2 T n , then it is guaranteed that, Let us denote y rno to be the result of RN Tn+2,rno (f (x)). Then, if our approximation function A(x) produces y rno for all inputs x 2 T n , then A(x) rounds to the correctly rounded result of f (x) for all representations T k with any rounding mode rm. Since any exactly representable value in T k are also exactly representable in T n , it is sufficient for us to make sure that A(x) produces y rno for all x 2 T n . Technically, the interval that contains a single value, [y rno , y rno ], is a generic interval for x. However, generating A(x) that produces a value in [y rno , y rno ] (where the only option is to produce y rno ) for all x 2 T n is infeasible. Thus, our next step is to find a range of values near y rno that is guaranteed to round to the correctly rounded result of f (x) for all T k and rm.

Formal Reasoning On the Odd Interval of v rno
Theorem 6.12 tells us that performing double rounding by first rounding v R to T n+2 using and then subsequently rounding v rno to T k using a rounding mode rm produces v rm , the value that we obtain when rounding v R directly to T k using the same rounding mode rm, Now, let us formally define the odd interval of v rno . The odd interval of v rno is a range of values [l rno , h rno ] where any value v in the interval rounds to v rno when rounded to T n+2 148 and then subsequently rounding v rno to T k using a rounding mode rm produces v rm , the value that we obtain when rounding v R directly to T k using the same rounding mode rm, Now, let us formally define the odd interval of v rno . The odd interval of v rno is a range of values [l rno , h rno ] where any value v in the interval rounds to v rno when rounded to T n+2 148 v 1 (odd) v 0 (even) (a) (b) (c) Fig. 10. An example to show that the round-to-odd result in T +2 maintains sufficient information to produce correctly rounded results for T when the round-to-odd result is double rounded to T with the mode. Real value is represented with a red star. Here, 0 , 1 , 2 , 3 , and 4 are representable values in T +2 . Values 0 and 1 are adjacent values representable in T . As all values in T are representable in T +2 , 0 = 0 and 4 = 1 . Solid arrow represents directly rounding the real value to T . The dotted arrows represent the process of double rounding from the real value to the round-to-odd result in T +2 and subsequently to T . (a) When the real value is in the interval between 0 and 2 . (b) When the real value is exactly equal to 2 , which is the midpoint of 0 and 1 . (c) When the real value is between 2 and 4 .
We propose an efficient procedure to create polynomial approximations H ( ) of an elementary function ( ) that produces values in the odd interval of the correctly rounded result in T +2 . Using Theorem 1, rounding any value in the odd interval (i.e., H ( )) to T using a rounding mode produces the correctly rounded result of ( ) in T using the same rounding mode .
An example to show why the round-to-odd result avoids double rounding errors. We provide intuition on how rounding with the round-to-odd mode avoids double rounding errors in Figure 10. Any value that is representable in T is also representable in T +2 . Further, there are three additional values ( 1 , 2 , 3 ) in T +2 between 0 and 4 . Here, 0 and 4 are also representable in T . In the round-to-odd mode, any real value between 0 and 2 rounds to 1 . Similarly, any real value between 2 and 4 rounds to 3 . If the real value is exactly equal to 0 , then the round-to-odd mode with T +2 also rounds to 0 (similarly 2 and 4 with T +2 ). Figure 10 illustrates the task of rounding the real value directly to T with the mode (solid arrow) and the result produced from double rounding the result from T +2 to T using the mode.
In the context of rounding a real value directly to T with rounding components, the last bit of the round-to-odd result in T +2 captures the sticky bit. Similarly, the penultimate bit of the round-to-odd result in T +2 captures the rounding bit. In summary, the round-to-odd result in T +2 maintains sufficient information about the real value so that when the round-to-odd result is (double) rounded to T with any rounding mode, it produces the correctly rounded result for T .

Polynomials for Correctly Rounded Results with the Round-to-Odd Mode in T +2
Our strategy is to create a generic polynomial approximation that produces correctly rounded results for T +2 using the round-to-odd mode. Next, we describe our approach to generate such a polynomial approximation. Algorithm 1 provides a high-level sketch of this process. Given an elementary function ( ) and a list of inputs in the T representation (i.e., ⊆ T ), the first step is to compute the correctly rounded result in representation T +2 using the round-to-odd mode (i.e., for each input ∈ ). Figure 11 shows our algorithm to compute the round-to-odd result for each input using the real value from the oracle. Subsequently, we compute the odd interval of each result such that any real value in the odd interval rounds to . Figure 11 also provides our algorithm to compute the odd interval. The odd intervals of some inputs can be a singleton (i.e., only one value in the odd interval), which we handle separately (Section 4.3). Once we have a set of non-singleton odd intervals for all inputs, we use our prior work in the RLibm project [Lim and Nagarakatte 2021a] to generate piecewise polynomials. Algorithm 1: A sketch of our approach to generate piecewise polynomials of degree for elementary function ( ) in the representation T +2 using the round-to-odd mode. The resulting polynomial when used with range reduction ( H ) and output compensation ( H ) produces correctly rounded results for all inputs ∈ with all representations T for all standard rounding modes. CalcResultsInRO computes the round-to-odd result using an oracle (see Figure 11). CalcOddIntervals computes the set of odd intervals ( ) and set ( ) of singleton odd intervals (see Figure 11). Once we have the odd intervals and singletons, we use RLibm's polynomial generation procedure (RLibmPolyGen) to obtain the generic polynomial.  Fig. 11. CalcResultsInRO computes the correctly rounded result of ( ) in T +2 using the round-to-odd rounding mode for each input ∈ . CalcOddIntervals computes the odd interval for each input based on the correctly rounded result in T +2 . The list is the set of inputs that have a singleton as the odd interval. The list contains inputs and the corresponding odd intervals. GetPrecVal( , T) returns the value preceding in the representation T. GetSuccVal( , T) returns the value succeeding in the representation T.
At the end of this process, we will have two main components that together can produce correctly rounded results for ( ). First, our approach produces a set that contains inputs whose odd interval is a singleton. For the resulting math libraries to be efficient, we need a fast method to check these inputs and compute results for them either using table lookups or using function-specific mathematical properties (Section 4.3). Second, our approach produces piecewise polynomials that when used with output compensation produces correct results for all inputs when rounded to any T with all standard rounding modes.
Computing the round-to-odd result from a real value. The first step in our approach is to identify the correctly rounded result for input . Figure 11 provides the steps to compute the round-to-odd result in T +2 given a real value of ( ). We compute the real value = ( ) for each input using an oracle (e.g., MPFR library). Then, we obtain the rounding components ( , − , , ) as described in Section 2.2. When the real value is exactly representable ( =0 and = 0) or when − is odd, then the round-to-odd result is − . Otherwise, the round-to-odd result in T +2 is the value succeeding − in T +2 .
Deducing the odd interval of an input. Once we determine the correctly rounded result of ( ) in representation T +2 using the round-to-odd mode, the next step is compute the interval of values in representation H, which is used for polynomial evaluation and range reduction, such that producing any value in the interval rounds to , which we call as the odd interval. The function CalcOddIntervals in Figure 11 describes the steps to compute the odd interval. If the correct rounded result in T +2 is even, then the odd interval is a singleton. In such cases, the only value that rounds to with the round-to-odd mode is itself. Generating polynomial approximations with singletons is challenging because they limit the amount of freedom available to the polynomial generator. Hence, we identify such inputs and handle them separately. If is odd, then all values in H that are strictly greater than the preceding value of in T +2 and strictly less than the succeeding value of in T +2 forms the odd interval. Any value in this odd interval rounds to in T +2 with the round-to-odd mode. We deduce the odd interval for each input. In Figure 11, represents the set of non-singleton odd intervals for all inputs, which is given to the polynomial generator.
Piecewise polynomial generation using the odd intervals. The next step is to generate piecewise polynomials that produce a value in the odd interval for all inputs. Each input and odd interval pair (i.e., ( , [ , ℎ]) ∈ specifies the constraints on the polynomial approximation H ( ) for each input . Using the RLibm methodology, we create an LP problem with these constraints to deduce the coefficients of a polynomial with degree . Similarly, we generate piecewise polynomials and use counterexample guided polynomial generation to facilitate the entire process. We also make sure that the generated polynomial produces a value in the odd interval after range reduction and output compensation.
Implementation of the polynomial approximation for T . At the end of polynomial generation, we will have a a set of inputs whose odd interval is a singleton and a polynomial approximation of ( ) that produces the correct round-to-odd result in T +2 for all inputs. We implement the polynomial approximation as follows. Given an input , we first check whether the input 's odd interval is a singleton. If so, we either use the precomputed round-to-odd result with table lookups or efficiently compute the round-to-odd result using function-specific properties. Otherwise, we use perform range reduction and use Horner's method for polynomial evaluation to compute the round-to-odd result in T +2 for the input. Finally, we round the round-to-odd result in T +2 to T using the user specified rounding mode to return the final result. We guarantee that our implementation produces the correctly rounded result of ( ) for any representation T with all the standard rounding modes for all inputs .

Computing Round-to-Odd Results for Inputs with Singleton Odd Intervals
One of the challenging issues for polynomial generation with odd intervals is the presence of singletons, which happens when the correctly rounded result with the round-to-odd mode is even (i.e., the round-to-odd result in T +2 matches the real value). We want to identify such inputs efficiently. As both T and T +2 are finite precision representations, all values in T and T +2 are rational values. If the real value matches the round-to-odd result in T +2 exactly, then it is a rational value. Hence, our task of identifying inputs with singleton odd intervals corresponds to the problem of identifying rational inputs that produce rational outputs for various elementary functions, which is well studied [Aigner and Ziegler 2009;Baker 1975;Cohn 1974;Niven 1956]. We first use the mathematical properties of the elementary function ( ) to identify all rational inputs such that ( ) is a rational value. Then, we check if these inputs 1 and the corresponding result ( 1 ) are exactly representable in T and T +2 , respectively. If so, such inputs are of interest. Then, we need to develop a quick way to identify those inputs and compute the round-to-odd results for them without using a multi-way branch.
We now describe the specific mathematical properties of elementary functions that we use to identify inputs whose odd interval is a singleton and the mechanism that we use to efficiently compute round-to-odd results for them.
Functions and ( ). From the Lindemann-Weierstrass theorem [Baker 1975], if the input is a non-zero rational value, then cannot be a rational value. Hence, the only value that will have a singleton odd interval with is = 0. Similarly, ( ) will produce a rational output only when = 1, which also follows from the Lindemann-Weierstrass theorem [Baker 1975]. We have a single branch to check this input and return the pre-computed correctly rounded result.
Functions 2 and 10 . The function 2 can produce a rational result only when is an integer and the value of 2 is less than the dynamic range of the T +2 representation. When T is a 32bit float, T +2 can represent all values of 2 for between −151 ≤ ≤ 127. Thus, any integer input between −151 and 127 (279 inputs in total) can produce a singleton odd interval. Hence, our implementation checks whether is an integer within a certain bound (i.e., −151 ≤ ≤ 127) and directly computes the result, 2 , using bit-wise operations.
Similar to 2 , 10 produces a rational value when is a positive integer. In contrast to 2 , 10 grows much faster and there are a few inputs for which 10 is exactly representable in T +2 . For a 32-bit float (T ), there are only 12 inputs ranging from 0 to 11 that are exactly representable in a 34-bit float (T +2 ). We use a precomputed table to store the correct results for these 12 inputs and use a switch statement for it.

Functions
2 ( ) and 10 ( ). The 2 ( ) function produces a rational result when is a power of 2 (i.e., = 2 and is an integer). We use bitwise operations to check if the input is a power of two. In contrast to 2 ( ), 10 ( ) produces a rational result when is a positive power of 10 (i.e., = 10 and is a positive integer). This difference between 2 ( ) and 10 ( ) is due to the fact that T cannot exactly represent negative powers of 10. When we are generating a polynomial to approximate the round-to-odd result with a 34-bit float (i.e., T +2 ), there are 11 inputs that can produce singletons. We create table-lookups for them.
The hyperbolic functions, ℎ( ) and ℎ( ). If the input is a non-zero rational value, then = ℎ( ) or = ℎ( ) cannot be a rational value using the Lindemann-Weierstrass theorem [Niven 1956]. Hence, the only input whose odd interval is a singleton is 0, for which we use a branch condition.
The ( ) function. The function ( ) is equal to ( ). By Niven's theorem [Niven 1956], the only rational values of between 0 ≤ ≤ 1 2 where ( ) is also a rational value are when = 0, = 1 6 , and = 1 2 . Among these three inputs, 1 6 is not exactly representable in T . Given that is a periodic function, there are only three cases of inputs in ∈ T where the result of ( ) is representable in T +2 when we extend the domain of x to the set of all inputs: We need to implement the floating point modulo operation efficiently using integer operations. Consider the case where T is a 32-bit float. All inputs ∈ T greater than or equal to 2 23 are integers where ( ) is always 0. Next, if < 2 23 , then we need to identify whether is either an integer, a multiple of 0.5, or a multiple of 1.5. To determine this condition, we compute 2 with a 32-bit float and then cast the result (i.e., 2 ) to a 32-bit integer to obtain the value . This operation of casting the value 2 to an integer truncates the value of 2 to the integral part of 2 . Now if we cast back to a 32-bit float value and the resulting float value is exactly equal to 2 , then is an integer, which implies that is either an integer, or a multiple of 0.5, or a multiple of 1.5. Finally, we compute the result of ( ) based on as shown below: The ( ) function. Similarly, ( ) = ( ) produces a rational value representable in T +2 in the following cases.
if is an even integer −1 if is an odd integer 0 if ( ) ≡ 0.5 These checks can be performed efficiently using a similar strategy illustrated for the ( ) function. In summary, handling singleton odd intervals efficiently is important for performance when we generate polynomials for correctly rounded results in T +2 with the round-to-odd mode.

PROOF THAT THE ROUND-TO-ODD RESULT WITH T +2 PRODUCES CORRECT
RESULTS FOR T We provide a proof of Theorem 1 in this section. We prove that the round-to-odd result in T +2 produced by our polynomial approximation when rounded to T with any of the standard rounding modes produces the correctly rounded result for T .

Unique Properties of the Round-to-Odd Result
We prove the unique properties of the round-to-odd result, which we subsequently use to prove Theorem 1. When we use R to represent the real value, we refer to it in the extended infinite precision representation. Lemma 1. The round-to-odd result in T +2 preserves the sign of R .
Proof. The value zero is representable in T +2 . The only value that rounds to zero is zero itself. Hence, all positive real values will round to a positive value in the round-to-odd mode. Similarly all negative real values will round to a negative value in the round-to-odd mode. □.
Proof. The result is created using the rounding components ( , − , , ). The round-to-odd mode preserves the sign of R in . Without loss of generality, we assume R is positive for the rest of the proof. Further, − is the truncated value of R (see Section 2.2). Hence, all the ( + 2)-bits of − and R are identical. After rounding with the round-to-odd mode, is either equal to − or the succeeding value of − in T +2 . We prove that the ( + 1)-bits of and R are identical by looking at the possible values of − and its relation to . First case, when − is odd. Then, = − . Hence, all the ( + 2)-bits of and R are identical. Second case, when − is even. Hence, the last bit of − is 0. Now, there are two sub-cases. (1) If = 0 and = 0, then = − . Hence, all the ( + 2)-bits of and R are identical. (2) If ≠ 0 or ≠ 0, then is equal the succeeding value of − . The only bit that changes between − and its succeeding value is the ( + 2) ℎ -bit. Hence, the first ( + 1)-bits of and R are identical. □.
Lemma 3. The ( + 2) ℎ -bit of is equal to the bitwise OR of all the bits of R starting from the ( + 2) ℎ -bit.
Proof. We prove this lemma using a strategy similar to Lemma 2. Intuitively, this lemma states that the last bit of is 0 if and only all bits starting from the ( + 2) ℎ -bit of R is 0. As the round-to-odd mode preserves sign, we assume R is positive for the rest of the proof without loss of generality. Let us say the rounding components for are ( , − , , ). In the round-to-odd mode with T +2 , will be equal to either − or a succeeding value of − in T +2 . Now, we look at the cases where − is odd and even to complete the proof.
In the first case, − is odd. Then, = − . The ( +2) ℎ -bit of − is 1. As − is a truncated value of R , the ( + 2) ℎ -bit of R is 1. Hence, the bitwise-OR of all bits of R starting from the ( + 2) ℎ -bit is 1, which is equal to the ( + 2) ℎ -bit of in T +2 . In the second case, − is even. There are two cases depending on the values of and . In the first sub-case, = 0 and = 0, then = − . The ( + 2) ℎ -bit of − is 0. So is the ( + 2) ℎ -bit of R . From the definition of rounding components for T +2 , is the value of the bit at position ( + 3) in R and is the bitwise-OR of bits starting from ( + 4) ℎ -bit in R . Hence, all the bits of R starting from the ( + 2) ℎ -bit are 0, which matches the ( + 2) ℎ -bit of . The next sub-case is when ≠ 0 or ≠ 0. In this case, is equal to the succeeding value of − in T +2 , which is odd. Hence, the ( + 2) ℎ -bit of is 1. Both and are not zeros, one of the bits starting from ( + 2) ℎ -bit in R is 1. Hence, the bitwise-OR of all bits starting from the ( + 2) ℎ -bit is 1, which matches the ( + 2) ℎ -bit of . □.
Proof. This lemma directly follows from the definition of rounding components in Section 2.2. Intuitively, this lemma states that identifying the correctly rounded result of the real value in T only depends on the rounding components and the rounding mode . □.

Proof that Double Rounding the Round-to-Odd Result Produces Correct Results
for all T We now sketch the proof of Theorem 1. We show that rounding a real value R to the FP representation T +2 = F +2, | | using the round-to-odd mode to produce and then subsequently rounding the result ( ) to T using a rounding mode produces the same value as rounding R directly to T using the same rounding mode , as long as | | + 1 < ≤ . More formally, we prove that Fig. 12. Rounding components while rounding R and to T . We show the bit-string of R in extended infinite precision representation. Note is a value in T +2 .
for rounding to T , then we prove that T , ( R ) = T , ( ) for all rounding modes . Hence, our strategy is to show that the rounding components for rounding R and to T are identical. As the the round-to-odd mode preserves the sign, we will consider R to be positive in the rest of the proof.
The representation of R in extended infinite precision representation ( R ) is as follows: . . Let us say ( 1 , − 1 , 1 , 1 ) are the rounding components for rounding R to T . Then, − 1 is the truncated value in T . While rounding to T , the rounding bit, 1 , is the ( + 1) ℎ -bit of R . The sticky bit, 1 is the bitwise-OR of all bits starting from the ( + 2) ℎ -bit of R .
Hence, all the rounding components for rounding R to T directly and rounding to T are identical. Hence, T , ( R ) = T , ( ) from Lemma 4. □ Theorem 1 directly follows from Theorem 2. Theorem 2 states that , which is produced by rounding a real value R to T +2 using the round-to-odd mode, rounds to the same value as if R is directly rounded to T using the same rounding mode . If we substitute R with the exact result of the elementary function ( ) for a given input ∈ T , then T , ( T +2 , ( ( ))) = T , ( ( )) Further, by definition, all values in the odd interval of in T +2 round to . Hence, any value in the odd interval rounds to the correctly rounded result for representations T using any rounding mode ∈ { , , , , }.

EXPERIMENTAL EVALUATION
We describe our prototype, experimental methodology, and the results of our experiments to check the correctness and performance of the generated polynomial approximations.
Prototype. The prototype, RLibm-All, is open-source and publicly available [Lim and Nagarakatte 2021c]. RLibm-All is a generator and a collection of correctly rounded implementations of polynomial approximation for multiple representations and rounding modes. RLibm-All contains ten functions that produce the correctly rounded result of ( ) for the 34-bit FP representation (i.e., T +2 ) with 8 bits of exponent (FP34) with the round-to-odd mode. As FP34 is not supported in hardware, RLibm-All maintains the FP34 result in double precision. RLibm-All's functions produce the correct result for all -bit FP representations with 8-bits for the exponent with all the rounding modes in the IEEE standard where 9 < ≤ 32. This includes 32-bit float, bfloat16, and tensorfloat32. RLibm-All uses the MPFR library [Fousse et al. 2007] with up to 1000 precision bits to compute the oracle value of ( ). RLibm-All uses SoPlex [Gleixner et al. 2012], an exact rational LP solver, to generate the coefficients of the polynomials with a time limit of five minutes. We limit the size of the LP formulation to contain up to fifty thousand reduced input and interval constraints. We use the range reduction and output compensation functions from the RLibm prototype. RLibm-All performs range reduction, polynomial evaluation, and output compensation using the double precision. The polynomial evaluation uses the Horner's method [Borwein and Erdelyi 1995] for efficiency.
Experimental methodology and setup. We compare RLibm-All's functions with Intel's libm, glibc's libm, CR-LIBM [Daramy-Loirat et al. 2006], and RLibm-32. Among these, CR-LIBM provides four implementations for each elementary function that produces the correctly rounded results in double precision with the , , , and mode, respectively. CR-LIBM does not provide implementations for the mode. RLibm-32 provides correctly rounded functions for a 32-bit float with the mode. To produce the result in a target representation T that is not natively supported by these libraries, we first convert the input in T to the representation supported by the library, use the elementary function, and round the result back to T. We perform our experiments on a 2.10GHz Intel Xeon Gold 6230R machine with 192GB of RAM running Ubuntu 18.04. We disabled Intel turbo boost and hyper-threading to minimize noise. All our libraries are compiled with O3 optimizations. We use Intel's libm from the oneAPI Toolkit and glibc's libm from glibc-2.33. The test harness for comparing glibc's libm, CR-LIBM, and RLibm-32 is built using the gcc-10 compiler with -O3 -static -frounding-math -fsignaling-nans flags. Because Intel's libm is only supported in the Intel's compiler, we built a test harness that compares Intel's libm against RLibm-All using the icc compiler with -O3 -static -no-ftz -fp-model strict flags to obtain as many accurate results as possible. To compare performance, we measure the number of cycles taken to compute the result for each input using rdtscp. We then measured the total time taken to compute the elementary function as the sum of the time taken by all inputs. Table 1 provides details on the properties of the polynomials generated by RLibm-All. Our attempt was to generate piecewise polynomials with degree less than or equal to 8. We also restricted the Table 1. Details about the generated polynomials. For each function, we show the time taken to generate the polynomial in minutes, the size of the piecewise polynomial, the maximum degree, the number of terms, and whether the generated polynomial produces correct results in FP34 using the round-to-odd mode for all inputs. number of sub-domains for the piecewise polynomials to 2 15 . The output compensation function for ℎ( ), ℎ( ), ( ), and cospi( ) uses two elementary functions and we generate two piecewise polynomials for each function. The , 2 , and 10 functions have both negative and positive reduced inputs. Hence, we create two piecewise polynomials: one for the negative reduced inputs and the other for the positive reduced inputs. As RLibm-All generates piecewise polynomials for FP34 with the round-to-odd mode, the number of sub-domains used in the resulting piecewise polynomials are bigger than RLibm-32. However, the degrees of the polynomial for each sub-domain was similar to RLibm-32. The amount of time taken to generate the piecewise polynomials ranged from approximately 2 hours to 9 hours. About 79% of the total time on average is spent in computing the oracle result using the MPFR library. In contrast, computing the intervals and generating the polynomials using the LP solver takes 15% and 5% of the total time on average, respectively.

Does RLibm-All Produce Correct Results?
We experimentally show that RLibm-All produces correctly rounded results for all rounding modes for multiple representations. We built a harness that checks if RLibm-All's functions produce correctly rounded results for all inputs with 161 different FP representations where the number of exponent bits ranged from 2 to 8 and the number of mantissa bits ranged from 1 to 23 bits (i.e., 23 * 7 = 161) for all five rounding modes. RLibm-All produces the correct results for all these representations with all standard rounding modes with all inputs. RLibm-All is the first efficient library that provides correctly rounded results for all rounding modes for a 32-bit float.
Correct results with all rounding modes for a 32-bit float. Table 2 reports the results of our experiments to check whether existing libraries produce correct results for a 32-bit float type. While RLibm-All produces the correct float results with all five rounding modes, mainstream libraries (glibc and Intel's libm) do not produce correctly rounded results for all rounding modes for all inputs for many of the elementary functions. When CR-LIBM's implementations, which produces correctly rounded double results with the mode for all inputs, is used to produce the results for a 32-bit float, it does not produce correctly rounded results for all inputs with several functions due to double rounding. CR-LIBM's , , and implementations produce correctly rounded float results for the , , and mode respectively. Double rounding with these three rounding modes using CR-LIBM's implementations do not generate wrong results.
RLibm-32 produces correctly rounded results for all inputs for the mode with a 32-bit float. However, it does not produce correct results for other rounding modes. In contrast, RLibm-All's Table 2. Ability to generate correct results for a 32-bit float for all inputs with each of the five standard rounding modes with RLibm-All, glibc's double libm, Intel's double libm, CR-LIBM, and RLibm-32. CR-LIBM provides separate correctly rounded functions for each rounding mode except the mode, which we use to check its correctness. ✓indicates that the library produces the correct result using a given rounding mode for all inputs. Otherwise, we use ✗. Table 3. Ability to generate correct results with tensorfloat32 for all inputs with various rounding modes. ✓indicates that the library produces the correct tensorfloat32 result using a given rounding mode for all inputs. Otherwise, we use ✗. Table 4. Generation of correctly rounded results for bfloat16 using the five standard IEEE-754 rounding modes , , , , and . We show the results with the elementary functions from RLibm-All, glibc's double libm, Intel's double libm, RLibm-32, and CR-LIBM. Then, we convert the output to bfloat16 values. ✓indicates that the library produces the correctly rounded bfloat16 result using a given rounding mode for all inputs. Otherwise, we use ✗. We use the separate approximation provided for each of the four rounding modes with CR-LIBM.  2020, 2021; Lim and Nagarakatte 2021a; Ziv 1991]. They are restricted to a specific representation and a rounding mode.
CR-LIBM [Daramy et al. 2003;Lefèvre et al. 1998] is a correctly rounded collection of elementary functions for double precision. It was developed using Sollya [Chevillard et al. 2010], which is a tool and a library for developing FP code. Sollya can generate polynomials of degree with coefficients in a representation used for the implementation (H) that has the minimum infinity norm [Brisebarre and Chevillard 2007]. Sollya uses a modified Remez algorithm to produce polynomials. It also computes and proves the error bound on the polynomial evaluation using interval arithmetic [Chevillard et al. 2011;Chevillard and Lauter 2007]. Metalibm [Brunie et al. 2015;Kupriianova and Lauter 2014] is a customization infrastructure also built using Sollya. MetaLibm is able to automatically identify range reduction and domain splitting techniques for some transcendental functions. It has been used to trade-off correctness and performance while approximating elementary functions for float and double precision types.
A modified Remez algorithm has also been used to generate polynomials that minimizes the infinity norm compared to an ideal elementary function [Arzelier et al. 2019]. It can be useful for generating correctly rounded results for a specific precision and a rounding mode when range reduction is not necessary. This paper builds on our prior work in the RLibm project [Lim 2021;Lim et al. 2020Lim and Nagarakatte 2021a,d] that creates polynomials using the correctly rounded value rather than the real value of the elementary function. Like the RLibm project, we structure the problem of generating polynomials as an LP problem. We also use RLibm's range reduction strategies. The RLibm project has generated correctly rounded libraries with the commonly used mode for multiple types: bfloat16, posit16, 32-bit float, and posit32. However, it is necessary to create individual polynomial approximation for each representation with each rounding mode with the RLibm project to avoid double rounding errors. In contrast, this paper shows that by generating polynomial approximations for T +2 with the round-to-odd mode, we can create a single polynomial approximation that works for multiple representations T with multiple rounding modes.

CONCLUSION
This paper proposes a novel method to generate a single polynomial approximation that produces correctly rounded results for multiple representations and rounding modes. The key idea is to create a polynomial approximation that produces the correctly rounded result for T +2 with the round-to-odd mode when the goal is to generate correct results for T , where ≤ , with all rounding modes. We address the issue of singletons while generating polynomials that approximate the correctly rounded result with the round-to-odd mode. We provide the first correctly rounded implementations of elementary functions for multiple representations. We believe that our results make a strong case for mandating correctly rounded results at least with any representation that has fewer than or equal to 32-bits.