Correctly rounded floating-point maths functions

1. Motivation

Floating-point maths functions are rounded after calculation. The rounding mode used is part of the floating-point state which is maintained per-thread. This introduces two problems:

  1. Changing rounding modes, for example for calculations that require correct rounding in a set of optimised expression evaluations, is unergonomic.
  2. There is a burden on library writers to restore the floating-point state to the condition it was in when a library function was called.

This is not a new problem, and the C standard in Annex F reserves the prefix cr_ for functions fully matching the IEC 60559 mathematical operations. This paper proposes adding five overload sets to the standard library for addition, subtraction, multiplication, division and square root calculation.

These functions guarantee that the operation will be carried out as if with infinite precision, and rounded using the round-nearest-to-even rounding mode, as specified in 60559:2020. The parameter type must satisfy the type trait std::numeric_limits<T>::is_iec559.

Other functions, such as trig functions and complex maths functions, are not available for correctly rounded implementations. This would require specifying an algorithm, which is outside the scope of this standard. It is hoped that cross-implementation third-party libraries will be developed which will offer solutions to this problem.

This solves problem 1 directly by providing a more ergonomic way of expressing intention: the client can explicitly state that they require correctly rounded calculation without having to ensure that the floating-point state is set appropriately.

This does not solve problem 2 directly, since this does not necessarily affect the floating-point state at all. However, it reduces the sources of error when correct rounding is important.

2. Draft wording

This draft wording is based on N4950.

In the <cmath> synopsis, after the declaration of sph_neumannl, insert the following:

// 28.7.7 Correctly rounded functions constexpr floating-point-type cr_add(floating-point-type x, floating-point-type y); constexpr floating-point-type cr_sub(floating-point-type x, floating-point-type y); constexpr floating-point-type cr_mul(floating-point-type x, floating-point-type y); constexpr floating-point-type cr_div(floating-point-type x, floating-point-type y); constexpr floating-point-type cr_sqrt(floating-point-type x);

Amend 27.7.1 paragraph 1 to read:

The contents and meaning of the header <cmath> are the same as the C standard library header <math.h>, with the addition of a three-dimensional hypotenuse function (28.7.3), a linear interpolation function (28.7.4), the mathematical special functions described in 28.7.6, and five correctly-rounded math functions (28.7.7).

Introduce a new section as follows:

28.7.7 Correctly rounded functions [cmath.crfunc]
constexpr floating-point-type cr_add(floating-point-type x, floating-point-type y);
1
#
Requires: std::numeric_limits<floating-point-type>::is_iec559 = true.
2
#
Returns: x + y, calculated as if using infinite precision, with the result rounded using round-nearest-to-even, as specified in 60559:2020.
constexpr floating-point-type cr_sub(floating-point-type x, floating-point-type y);
3
#
Requires: std::numeric_limits<floating-point-type>::is_iec559 = true.
4
#
Returns: x - y, calculated as if using infinite precision, with the result rounded using round-nearest-to-even, as specified in 60559:2020.
constexpr floating-point-type cr_mul(floating-point-type x, floating-point-type y);
5
#
Requires: std::numeric_limits<floating-point-type>::is_iec559 = true.
6
#
Returns: x * y, calculated as if using infinite precision, with the result rounded using round-nearest-to-even, as specified in 60559:2020.
constexpr floating-point-type cr_div(floating-point-type x, floating-point-type y);
7
#
Requires: std::numeric_limits<floating-point-type>::is_iec559 = true.
8
#
Returns: x / y, calculated as if using infinite precision, with the result rounded using round-nearest-to-even, as specified in 60559:2020.
9
#
Remarks: if y = 0 then the result is undefined.
constexpr floating-point-type cr_sqrt(floating-point-type x);
10
#
Requires: std::numeric_limits<floating-point-type>::is_iec559 = true.
11
#
Returns: √x, calculated as if using infinite precision, with the result rounded using round-nearest-to-even, as specified in 60559:2020.
12
#
Remarks: if x < 0 then the result is undefined.

Proposal P3375 seeks to introduce reliable reproducibility to floating-point operations regardless of platform. This proposal partially addresses this problem by reducing the places where implementations can diverge. However, it does not necessarily solve the divergence introduced by inlining or optimisation choices, which are declared per translation unit.

This can be mitigated, at the cost of performance, by creating a struct which contains a value of the appropriate type as a single member, and implementing the arithmetic operators in a separate library, linked at runtime, in terms of the correctly rounded functions. In this way a client can mix correctly rounded and optimised operations in a single translation unit.

This author currently believes that more performant reproducibility can only be achieved by introducing a new type which is defined to be immune to the floating-point optimisations specified in 60559:2020.

4. References and bibliography

[60559] ISO/IEC JTC 1/SC 25 (May 2020). ISO/IEC 60559:2020 — Information technology — Microprocessor Systems — Floating-Point arithmetic. ISO. pp. 1–74.

[C++] C++23 Working Draft https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2023/n4950.pdf

[LIBMCR] Correctly rounded libm https://github.com/simonbyrne/libmcr

[OBILTSCHNIG] Cross-Platform Issues With Floating-Point Arithmetics in C++, Günter Obiltschnig, https://www.appinf.com/download/FPIssues.pdf

[P3375R3] Reproducible floating-point results, Guy Davidson, https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2025/p3375r3.html