Round to Odd
When writing numerical programs, rounding is a subtle but important aspect that can significantly affect the accuracy and stability of computations. In...
Added:
Last updated:
This blog post was split into two parts. This blog post covers rounding in detail and discusses some theory; this blog post focuses on design principles for number libraries.
Number libraries simulate number systems — floating-point, fixed-point, posits, and more — beyond those offered by standard hardware or language runtimes. Well known examples of these libraries include MPFR for arbitrary-precision floating-point numbers, GMP for arbitrary-precision integers and rationals, SoftFloat for emulating IEEE 754 floating-point numbers in software, and Universal for supporting formats across the machine learning landscape. They are essential tools for analyzing numerical error in multi-precision numerical algorithms, exploring different implementation trade-offs, and verifying the correctness of numerical software and hardware.
This flexibility over number systems comes at a cost: number libraries are complex, requiring expert knowledge to build and maintain, and significant effort to provide useful features and ensure correctness. Maintainers of these libraries face significant challenges as their tools serve as trusted reference implementations but must also be efficient for simulation and verification tasks. Due to intense demand from machine learning, there has been a proliferation of new number formats and rounding behaviors. Each new combination of operation, number format, rounding mode, overflow behavior, and special value handling requires careful implementation and testing, multiplying the complexity and maintenance burden. Developers must choose: stick with a smaller set of features or invest significant effort to meet user demand.
I maintain number libraries that fall into the latter category: they aim to support a wide variety of number formats and rounding behaviors. I want to reflect on two design principles that have helped mitigate some of these challenges.
First,
the round-to-odd rounding mode [1] allows decoupling
of arithmetic operations from rounding.
Applying this insight to number libraries,
a number library may consist of two independent parts:
a core arithmetic engine performing round-to-odd operations,
and a core rounding library that safely re-rounds under
the desired number format or rounding mode.
As a result,
each operation provided by the number library
is the composition of an operation in the arithmetic engine
and the round method of a rounding context instance from
the core rounding library.
This trick was pioneered by Bill Zorn [2]
in his number library found in the Titanic repo [3].
Second,
the core rounding library should be organized
into rounding contexts which encapsulate number format,
rounding mode, and other rounding information;
providing a single round method.
These rounding contexts are often composable:
a rounding context for an IEEE 754 floating-point number
can reuse the same logic as a rounding context for
a p-bit floating-point number,
with added checks for overflow.
Thus,
the rounding library can be broken down
further into composable components that can
be assembled to implement rounding for a
variety of number formats.
Put together, these two design principles enable significant modularity and code reuse within number libraries, resulting in several benefits. They dramatically improve maintainability: the library is smaller; it consists of smaller, composable components rather than a monolithic implementation where each combination of operation and rounding mode requires separate code. They ensure extensibility: adding features is easier and less error-prone; adding a new operation requires implementing it only once in the arithmetic engine, while adding a new format or rounding mode requires implementing it only once in the rounding library. They improve correctness: testing can be done compositionally; components can be extensively tested in isolation, and their composition inherits the correctness guarantees.
I will explore these two design principles in greater detail and illustrate their benefits through examples. As someone working in this space, I hope this blog post provides some insight into the challenges in this domain and possible solutions.
The first design principle is to separate rounding from arithmetic operations using round-to-odd arithmetic [1]. My blog post covers this topic in detail. I’ll summarize the key ideas here and illustrate how they apply to number libraries.
For a real number function \(f\), we say that an implementation of \(f\), say \(f^{*}\), is correctly-rounded when it produces the infinitely precise result of \(f\) rounded to the target number format. This rounding operation is usually parameterized by a rounding mode that specifies which representable value to choose when the result is not exactly representable in the target number format; the IEEE 754 standard specifies several such rounding modes. For now, we’ll only consider floating-point numbers with a fixed precision of \(p\) digits, so I’ll denote the rounding operation in the style of Boldo and Melquiond [1] as \(\mathrm{rnd}^{p}_{rm}\). Using this notation, the correctly-rounded implementation of \(f\) in precision \(p\) under rounding mode \(rm\) can be expressed as:
\[f^{*} = \mathrm{rnd}^{p}_{rm} \circ f.\]Boldo and Melquiond [1] describe a non-standard rounding mode called round to odd and prove that the rounding mode permits safe re-rounding. Rounding with round to odd at precision \(p + k\) (for \(k \geq 2\)) followed by re-rounding under a standard rounding mode at precision \(p\), yields the same result as rounding directly at precision \(p\) under the desired rounding mode.
Theorem 1. Let \(p, k \geq 2\) be integers; and \(rm\) be a standard rounding mode. Then,
\[\mathrm{rnd}^{p}_{rm} = \mathrm{rnd}^{p}_{rm} \circ \mathrm{rnd}^{p+k}_{\text{RTO}}.\]Applying this result to the definition of a correctly-rounded implementation, we can derive that \(f^{*}\) is the composition of (i) a round-to-odd implementation of \(f\) at higher precision, followed by (ii) re-rounding under the desired rounding mode:
\[f^{*} = \mathrm{rnd}^{p}_{rm} \circ f_{\text{RTO}}^{*}.\]The result is not novel: one successful application is found in the RLibm project [4] which automatically generates efficient, correctly-rounded elementary functions by generating a polynomial approximation with additional bits of precision using round-to-odd arithmetic that will be correctly rounded when re-rounded under the desired rounding mode.
For number libraries, the result has significant implications. It suggests that arithmetic and rounding can be decoupled: a number library can consist of two independent components: an arithmetic engine that implements each mathematical operation using round-to-odd arithmetic, and a rounding library that implements rounding operations for various number formats and rounding modes. The only interaction between the two components is that the rounding library must provide the precision \(p + k\) required for safe re-rounding.
To illustrate this approach in practice, consider a number library providing an implementation of multiplication.
module Engine:
def rto_mul(x, y, p):
...
module Round:
def round(x, p, rm):
...
def rto_prec(p):
...
def mul(x, y, p, rm):
rto_p = Round.rto_prec(p)
result = Engine.rto_mul(x, y, rto_p)
return Round.round(result, p, rm)
The Engine module implements arithmetic operations
that produce round-to-odd results with at least precision p.
The Round module handles rounding operations
and precision calculations:
round rounds x to precision p using the specified rounding mode,
while rto_prec calculates the precision needed for safe re-rounding.
The mul function provided by the number library
composes these functions by
performing round-to-odd multiplication,
and re-rounding to the desired precision and rounding mode.
To implement rto_mul,
we can use existing libraries like MPFR
that have been extensively tested.
MPFR provides a narrower interface
that performs floating-point arithmetic at
a specified precision and rounding mode.
Although MPFR does not directly support round-to-odd,
we can implement it as described by Boldo and Melquiond [1]:
we use MPFR’s implementation at precision p - 1 with
round towards zero, followed by an additional step to adapt
the result to be round-to-odd at precision p.
def rto_mul(x, y, p):
r = MPFR.mul(x, y, p - 1, MPFR.RTZ)
return rto_fixup(r)
Note by Theorem 1,
we can request higher precision than p for safe re-rounding
at precision p - 2.
Since the arithmetic engine and rounding library are separate,
the exact implementation of rto_mul can be changed
without affecting the correctness of any mul implementation,
as long as it produces correct round-to-odd results.
For example,
assume that floating-point numbers in our library
have the following structure:
class Float: # (-1)^sign * c * 2^exp
sign: bool # sign
exp: int # exponent
c: int # significand (c >= 0)
We can implement rto_mul manually as follows:
def rto_mul(x, y, p):
s = x.sign != y.sign
exp = x.exp + y.exp
c = x.c * y.c
return Float(s, exp, c)
Notice that this implementation does not
perform any rounding, so p is unused.
Since the result is infinitely precise,
it can be re-rounded under any precision.
Extending this implementation to handle special values
(\(+\infty\), \(-\infty\), NaN, etc.)
just requires careful case analysis.
For the Round module,
we need to implement two functions:
rto_prec(p) computes the precision required
for safe re-rounding to precision p; and
round(x, p, rm) rounds x to precision p
using the specified rounding mode rm.
The implementation of rto_prec is straightforward:
it returns p + k for a constant k >= 2.
Implementing the round function
requires care as it’s the core method of the rounding library.
The exact implementation is too verbose
to include here, but I’ll outline its basic structure.
For this example,
we’ll ignore special values like NaN and infinity.
def round(x, p, rm):
n = x.e - p # compute where to round off digits
hi, lo = x.split(n) # split into significant and leftover digits
rbits = lo.round_bits(n) # summarize leftover digits for rounding decision
increment = decide_increment(hi, rbits, rm) # round away from zero based on rounding mode?
if increment: # adjust hi if needed
hi = hi.increment()
return hi
The round function first determines
where to round off digits based on the
normalized exponent x.e of the argument and
the target precision p.
Any digit above this point is significant
while digits below must be rounded off.
A method split divides the number based on this point,
and the lower part is summarized into rounding bits.
Since hi represents the round-towards-zero result,
we decide whether to round away from zero to the next representable
value based on the rounding bits and the specified rounding mode.
The correctly-rounded result is hi.
A number library built in this style
can be extended to support new operations
and rounding modes easily.
For a new operation \(g\),
we must only implement \(g^{*}_{\text{RTO}}\)
in the arithmetic engine,
either from existing libraries or manually.
The function round in the rounding library
implements the function \(\mathrm{rnd}^{p}_{rm}\).
Their composition is then
the correctly-rounded implementation \(g^{*}\) of \(g\).
For a new rounding mode,
only the round function needs to be updated;
the arithmetic engine remains unchanged.
Compare this approach to a monolithic implementation.
Each new operation requires implementing
operation-specific rounding logic;
each new rounding mode requires updating
every operation.
To briefly summarize, separating arithmetic from rounding achieves the following benefits:
Maintainability: The arithmetic engine implements each operation only once, while the rounding library implements rounding logic separately.
Extensibility: Any new mathematical operation can be composed with the existing rounding library. Similarly, any new rounding logic can be composed with the existing arithmetic engine.
Correctness: Once a mathematical operation is verified in the arithmetic engine, its correctness applies to any composition with the rounding library. Similarly, once a rounding context is verified, its correctness applies to any mathematical operation. Thus, testing can be done modularly and the results reused.
While separating rounding from arithmetic allows developers to split the number library into two independent components, the rounding library must support a large number of number formats and rounding behaviors. To manage this complexity, we can organize the rounding library into rounding contexts.
A rounding context encapsulates all information
required to round a number correctly:
the number format (precision, exponent range, etc.),
the rounding mode,
overflow behavior (saturating, wrapping, exceptions, etc.),
and special value handling (NaN, infinity, etc.).
Rather than having a single round function
that consumes all rounding information as
a exhaustive list of parameters,
we convert Round into an interface;
each rounding context will implement the interface.
interface Round:
def rto_prec():
...
def round(x):
...
def round_core(x, p, rm): # default method
...
def mul(x, y, ctx):
rto_p = ctx.rto_prec()
result = Engine.rto_mul(x, y, rto_p)
return ctx.round(result)
The two interface methods are similar to before:
rto_prec computes the precision required to safely re-round
a number under the rounding context; andround rounds x under the rounding context.The previous round implementation is now round_core.
The mul function now accepts
a rounding context instance rather than explicit rounding parameters:
its implementation must be adapted slightly.
One strategy for implementing rounding contexts is
to organize each implementation of Round
based on families of number formats.
For example,
we can implement p-digit floating-point numbers,
as a family of rounding contexts.
class MPFloat(Round):
p: int # p >= 1
rm: RoundingMode
def rto_prec(self):
return self.p + 2
def round(self, x):
return self.round_core(x, self.p, self.rm)
Notice that we completely
reuse the round_core implementation for the rounding logic.
Critically,
we add support for new number formats
by simply implementing a new class that
implements the Round interface.
For example,
consider supporting the IEEE 754 floating-point numbers
parameterized by es, the size of the exponent field,
and nbits, the total number of bits of the representation.
While IEEE 754 numbers have restrictions
on exponent ranges, its core rounding behavior
is similar to p-digit floating-point numbers:
we can reuse the MPFloat rounding logic.
The exact implementation is too verbose,
but an outline of the implementation is as follows:
class IEEEFloat(Round):
es: int # exponent size (es >= 2)
nbits: int # total number of bits (nbits >= es + 2)
rm: RoundingMode
def emin(self): # minimum (normalized) exponent
return 1 - (1 << (es - 1))
def rto_prec(self):
p = self.nbits - self.es
return MPFloat(p).rto_prec() # need at least this much precision
def round(self, x):
max_p = self.nbits - self.es # maximum allowable precision
e_diff = x.e - self.emin() # e_diff < 0 if subnormal
p = min(max_p, max_p + e_diff + 1) # adjust precision for subnormals
r = MPFloat(p, rm).round(x) # re-use rounding logic
# handle overflow based on rounding mode
...
return r
The implementation of round is more complicated.
First,
we must deal with subnormal numbers,
i.e., numbers with magnitude below \(2^{emin}\)
which have reduced precision:
the min(max_p, ...) expression
computes the effective precision.
We adjust the precision accordingly,
construct the appropriate MPFloat context,
and reuse its round method to round
without exponent bounds.
Finally,
we handle overflow based on the specified rounding mode.
What about fixed-point formats?
To support fixed-point numbers,
we first alter round_core to accept a parameter n
which sets n directly:
def round_core(x, p, n, rm): # added n parameter
...
The caller must specify either p or n.
The arithmetic engine must also be altered
to support a stopping point n rather than
precision p when performing fixed-point operations.
For example,
if we want to round with n=-1, keeping only
integer digits, then the arithmetic engine must produce
a round-to-odd result with enough arbitrary precision
to preserve all integer digits plus extra digits
for safe re-rounding.
Like before,
one of p or n must be specified.
module Engine:
def rto_mul(x, y, p, n):
...
def mul(x, ctx):
p, n = ctx.rto_params()
result = Engine.rto_mul(x, y, p, n)
return ctx.round(result)
Consider implementing a rounding context
for an arbitrary-precision fixed-point number
that must round any digit less significant than
the n+1-th digit:
class MPFixed(Round):
n: int # first insignificant digit (drop digits <= n)
rm: RoundingMode
def rto_params(self):
return None, self.n - 2 # 2 additional bits for safe re-rounding
def round(self, x):
return self.round_core(x, None, self.n, self.rm)
I could continue to list implementations of
rounding contexts for other number formats,
but I believe the pattern is clear:
to support new number formats, we implement
a rounding context class that can often reuse
existing rounding logic.
Machine integers and fixed-point formats
can compose MPFixed and apply the appropriate
overflow behavior.
Floating-point formats like those
described in the OCP MX standard [5] might
benefit from a similar approach to
IEEE 754 floating-point numbers.
Posit numbers [6] are floating-point numbers
with tapered precision;
one approach might be to use MPFloat
with variable precision based on the magnitude of the number.
Organizing the rounding library into rounding contexts achieves the following benefits:
Maintainability: Rounding contexts encapsulate rounding logic and often compose together existing rounding contexts to implement its rounding behavior.
Extensibility: Implementing a new number format only requires implementing a new rounding context class, often reusing existing rounding logic; any instance of the new rounding context can be used with all existing mathematical operations.
Correctness: Rounding can be decomposed into many smaller, often reusable components; verifying each component in isolation increases confidence in the correctness of each rounding context implementation.
To demonstrate the benefits of these design principles, I applied them to the number library underlying the runtime of the FPy language [7]. FPy is an embedded Python DSL for specifying numerical algorithms, with explicit control over rounding via rounding contexts, including first-class rounding context values. The number library supports many families of number formats from IEEE 754 floating-point numbers, OCP MX floating-point numbers, fixed-point numbers, and more. The core arithmetic engine uses MPFR [8] to implement round-to-odd arithmetic operations.
Due to its modular design, the FPy number library remains compact, composing core components to support a wide variety of number formats and operations.
The FPy number library consists of four major components:
Number module defines FPy’s number representation;Rounding module implements the round_core;Arithmetic module implements round-to-odd arithmetic using MPFR;Contexts module implements various rounding contexts.| Component | LOC |
|---|---|
| Number | 1725 |
| Rounding | 400 |
| Arithmetic | 1500 |
| Contexts | 4350 |
| Total | 7975 |
The total code size is approximately 8,000 lines of code (LOC)
with 4750 lines dedicated to rounding alone.
Since FPy’s arithmetic engine uses MPFR,
the arithmetic engine is relatively small at 1500 lines of code.
Rounding contexts in FPy implement more
than the essential round method
with additional methods for encoding and decoding bit patterns,
constructing special values,
and more.
The largest rounding context is almost 850 lines of code,
while the smallest is only 60 lines of code.
An exhaustive list of rounding contexts is below:
| Rounding Context | Parameters | Description |
|---|---|---|
MPFloat |
p | p-digit floating-point number |
MPSFloat |
p, emin | MPFloat with minimum exponent |
MPBFloat |
p, emin, max | MPSFloat with maximum value |
IEEEFloat |
es, nbits | IEEE 754 floating-point number |
EFloat |
es, nbits, I, O, E | generalized IEEE 754 format [9] |
MPFixed |
n | unbounded fixed-point number |
MPBFixed |
n, max | fixed-point with maximum value |
Fixed |
scale, nbits | nbits fixed-point with scale \(2^{scale}\) |
SMFixed |
scale, nbits | sign-magnitude nbits fixed-point |
ExpFloat |
nbits | nbits exponential floating-point number |
The MPFloat and MPFixed rounding contexts
call the core rounding logic directly,
while other rounding contexts
compose existing rounding contexts
to implement their rounding behavior.
Every rounding context may be used
with any arithmetic operation
provided by the arithmetic engine.
This compositional design
ensures the codebase remains relatively small,
comprising of modular components
rather than monolithic implementations
where small changes require large modifications
or modifications across many parts of the codebase.
To demonstrate extensibility of FPy,
I will showcase the implementation of
a correctly-rounded implementation of 1/x^2
and the ExpFloat number format.
1/x^2A correctly-rounded implementation 1/x^2
provides additional accuracy compared
to composing 1/x with x^2, separately.
To implement this operation in FPy,
we only need to implement its round-to-odd implementation.
To illustrate one such implementation,
I implemented a digit recurrence algorithm
that iteratively computes the significant digits of 1/x^2.
The implementation adapts the classic reciprocal digit-recurrence algorithm. Assuming that
\[x = {(-1)}^s * 1.m * 2^{e} = {(-1)}^s * c * 2^{exp},\]we know the result is of the form
\(q * 2^{-2e}\) where \(q = 1/(1.m)^2\).
The algorithm first computes the square of the argument:
adding the exponents and squaring the significand.
Then,
it computes the expected exponent, both \(e\) and \(exp\),
and checks for the special case where the (fractional)
significand is exactly 1.0, in which case the result
is just \(2^{-2e}\).
Otherwise,
it performs digit-recurrence to compute all
but the last significant digit of the result.
Finally,
the last digit is determined by round-to-odd:
if we have a non-zero remainder,
we need to round up to the result with
odd significand.
def rto_recip_sqr(x, p):
assert x != 0
# square the argument
x = x * x
e = -x.e # result normalized exponent
exp = e - p + 1 # result unnormalized exponent
m = x.c # argument significand (in 1.M)
one = 1 << (x.p - 1) # representation of 1.0 (fixed-point)
if m == one:
# special case: m = 1 => q = 1.0
q = 1 << (p - 1)
else:
# general case: m > 1 => q \in (0.5, 1.0)
# step 1. digit recurrence algorithm for 1/m
# trick: skip first iteration since we always extract 0
q = 0 # quotient
r = one << 1 # remainder (constant fold first iter)
for _ in range(1, p): # compute p - 1 bits
q <<= 1
if r >= m:
q |= 1
r -= m
r <<= 1
# step 2. generate last digit by inexactness
q <<= 1
if r > 0:
q |= 1
# step 3. adjust exponent so that q \in [1.0, 2.0)
exp -= 1
# result
return Float(s, exp, q)
The implementation is clearly naive,
but it verifiably satisfies the round-to-odd contract
of the arithmetic engine.
Composing this implementation
with any rounding context in FPy
yields a correctly-rounded implementation of 1/x^2
under that number format and rounding mode.
For example,
computing with single-precision floating-point:
the computed result is 0.101321185
compared to 0.101321183 when computed
separately as 1/x and x^2.
ExpFloatThe ExpFloat number format
encodes exponential numbers, i.e., 2^{e}
where e is an integer exponent,
that is, positive, power-of-two numbers.
Unlike standard floating-point numbers,
ExpFloat numbers cannot be negative,
zero, infinity, or NaN.
In the OCP MX standard [5],
the “E8M0” format is an 8-bit exponential number
encoding the possible values of the exponent field
of a single-precision IEEE 754 floating-point number.
The ExpFloat rounding context
is parameterized by nbits, the total number of bits
in the representation.
To implement ExpFloat in FPy,
we compose with the MPFloat rounding context implementation:
class ExpFloat(Round):
nbits: int # total number of bits (nbits >= 2)
rm: RoundingMode # rounding mode
... # overflow/underflow behavior
def mp_ctx(self): # corresponding MPFloat context
return MPFloat(1, rm)
def rto_params(self):
return self.mp_ctx().rto_params()
def round(self, x):
# user-defined behavior for negative, zero, Inf, or NaN
if x.is_nar() or x <= 0:
...
r = self.mp_ctx().round(x) # re-use rounding logic
# handle overflow/underflow based on rounding mode
...
return r
The implementation of ExpFloat is highly compact:
it reuses the MPFloat rounding logic entirely,
wrapping it with custom behavior for invalid values,
overflow, and underflow behavior.
In FPy,
the actual implementation of the round method
of the ExpFloat context comes out to 50 logical lines
of code (110 in total), with most lines dedicated
to overflow and underflow handling based on
the specified rounding mode;
the full context implementation with encoding logic,
value constructors, and predicates comes out to 500 lines.
Rather than having to test each operator implementation in its entirety, FPy’s design allows testing to be done compositionally. Each arithmetic operation in the arithmetic engine can be tested independently of the rounding library; the rounding library can be tested piecewise by verifying each rounding context or core rounding logic in isolation. In particular, I use property-based testing with the Hypothesis library [10].
Ensuring correctness of the core
rounding procedure round_core is especially important,
as it is the basis for all rounding contexts.
To test round_core,
we can verify that it satisfies
the expected properties of rounding.
For example,
we expect round_core to satisfy two properties:
round_core(x, p, None, rm) has at most p significant digits;round_core(x, None, n, rm) has no significant digits below the n+1-th digit.@given(floats(prec_max=256), st.integers(min_value=0, max_value=1024), rounding_modes())
def test_round_p(x, p, rm):
y = round_core(x, p, None, rm)
assert 0 <= y.p <= p
@given(floats(prec_max=256), st.integers(min_value=-1024, max_value=1024), rounding_modes())
def test_round_n(x, n, rm):
y = round_core(x, None, n, rm)
_, lo = y.split(n)
assert lo == 0
The code above checks the two properties.
The methods floats and rounding_modes
are generators for FPy floating-point number values
and rounding modes, respectively.
We could test the claim of the round method
that for hi, lo = x.split(n),
the value of hi is the round to zero result.
Assuming that x is finite,
we can verify this property as follows:
@given(floats(prec_max=256), st.integers(min_value=-1024, max_value=1024))
def test_round_rtz(x, n):
hi, _ = x.split(n)
y = round_core(x, None, n, RoundingMode.RTZ)
assert y == hi
While these three properties aren’t explicitly
checking the numerical correctness of round_core,
they provide confidence that the implementation
behaves as expected.
Similar property tests can be devised
to further test the correctness of round_core
and combined with concrete test cases
or other testing methods to increase confidence
in its correctness.
Clearly,
we would also want to verify the split helper method:
that it does in fact split the number into two non-overlapping
groups of digits at the specified point.
@given(floats(prec_max=256), st.integers(min_value=-1024, max_value=1024))
def test_split(x, n):
hi, lo = x.split(n)
# check we did not lose information
assert hi + lo == x, "split parts must sum to original"
# check we did not gain information
assert hi.p <= x.p, "hi must not have more precision than x"
assert lo.p <= x.p, "lo must not have more precision than x"
# check non-overlapping property
assert hi.exp > n, "LSB of hi must be > n"
assert lo.e <= n, "MSB of lo must be <= n"
Critically,
once round_core is verified,
the effort of verifying each rounding context
mostly involves verifying the unique behavior
of that rounding context;
confidence in the core rounding logic
carries over to each rounding context implementation.
Testing of the arithmetic engine can be done independently of the rounding library. FPy relies on MPFR for round-to-odd arithmetic, which has been extensively tested; the RLibm system also uses MPFR for verifying its transcendental function implementations [4]. Fully verifying custom transcendental function implementations is difficult work, far beyond the scope of property-based testing.
However, we can still verify that each arithmetic operation round-to-odd implementations meet the round-to-odd contract.
@given(floats(prec_max=256), floats(prec_max=256), st.integers(min_value=2, max_value=1024))
def test_rto_mul(x, y, p):
r_rtz = MPFR.mul(x, y, p - 1, MPFR.RTZ)
r_rto = rto_mul(x, y, p)
assert r_rtz.p == p - 1, "MPFR result must have p - 1 precision"
assert r_rto.c % 2 == (1 if r_rtz.inexact else 0), "inexact iff LSB is odd"
For arithmetic,
our custom rto_mul implementation
can be verified against Python’s native Fraction implementation
for finite inputs.
@given(floats(prec_max=256), floats(prec_max=256), st.integers(min_value=2, max_value=1024))
def test_rto_mul(x, y, p):
r_ref = x.as_rational() * y.as_rational()
r_impl = rto_mul(x, y, p)
assert r_ref == r_impl
Verifying both the arithmetic engine and rounding library in a compositional manner increases confidence in the correctness of the overall number library.
Number libraries face an explosion of complexity with new formats, rounding modes, and operations. The two design principles presented here — separating arithmetic from rounding via round-to-odd, and organizing rounding logic into composable contexts — offer a practical solution to this challenge. By decoupling these concerns, number library developers can achieve smaller, more maintainable codebases that are more extensible and easier to verify compared to traditional monolithic designs. As the numerical computing landscape continues to evolve, these principles are one approach to providing a foundation for building maintainable, correct, and extensible number libraries that can adapt to future requirements without overwhelming their maintainers.
Sylvie Boldo, Guillaume Melquiond. When double rounding is odd. 17th IMACS World Congress, Jul 2005, Paris, France. pp.11. ffinria-00070603v2f
Bill Zorn. 2021. Rounding. Ph.D. Dissertation. University of Washington, USA. https://hdl.handle.net/1773/48230
Bill Zorn. 2017. Titanic [GitHub]. https://github.com/billzorn/titanic. Accessed on November 11, 2025
Jay P. Lim and Santosh Nagarakatte. 2022. One Polynomial Approximation to Produce Correctly Rounded Results of an Elementary Function for Multiple Representations and Rounding Modes. Proc. ACM Program. Lang. 6, POPL, Article 3 (January 2022), 28 pages. https://doi.org/10.1145/3498664.
Bita Darvish Rouhani, Ritchie Zhao, Ankit More, Mathew Hall, Alireza Khodamoradi, Summer Deng, Dhruv Choudhary, Marius Cornea, Eric Dellinger, Kristof Denolf, Stosic Dusan, Venmugil Elango, Maximilian Golub, Alexander Heinecke, Phil James-Roxby, Dharmesh Jani, Gaurav Kolhe, Martin Langhammer, Ada Li, Levi Melnick, Maral Mesmakhosroshahi, Andres Rodriguez, Michael Schulte, Rasoul Shafipour, Lei Shao, Michael Siu, Pradeep Dubey, Paulius Micikevicius, Maxim Naumov, Colin Verrilli, Ralph Wittig, Doug Burger, and Eric Chung. 2023. Microscaling Data Formats for Deep Learning. arXiv preprint arXiv:2310.10537 (2023). https://arxiv.org/abs/2310.10537
John L. Gustafson. 2017. Beating Floating Point at its Own Game: Posit Arithmetic. Supercomputing Frontiers and Innovations 4, 2 (2017), 71–86. https://doi.org/10.14529/jsfi170206
Brett Saiki. 2025. FPy [GitHub]. https://github.com/bksaiki/fpy. Accessed on November 12, 2025
MPFR Team. 2025. The GNU MPFR Library. https://www.mpfr.org/. Accessed on November 14, 2025.
Brett Saiki. 2025. Taxonomy of Small Floating-Point Formats. https://uwplse.org/2025/02/17/Small-Floats.html. Accessed on November 12, 2025
Hypothesis Team. 2025. Hypothesis. https://hypothesis.works/. Accessed: 2025-11-12.
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
When writing numerical programs, rounding is a subtle but important aspect that can significantly affect the accuracy and stability of computations. In...
This blog post was split into two parts. This blog post covers rounding in detail and discusses some theory; this blog post focuses on design principles ...
Herbie’s improvement loop has five basic phases:
September 20th marked the one year anniversary of Minim’s creation. It initially started as an pandemic-fueled side project based on an online guide of mak...
Currently, one of the key issues in Minim is the lack of garbage collection. To handle precise memory management, I implemented a weird owner-and-reference ...
Recently, I released version 0.2.0 of Minim, my hobby-language that I’ve been developing since last fall. It’s inspired by my time working with Racket (now m...
As a side effect of recent work, I created an alternate MPFR interface in Racket. I posted in the previous blog, that I was planning on extracting that code ...
Note: This is my first entry for my blog, a compendium of my thoughts, articles, references, etc. Is blog even the correct word? I’m still trying to figure o...