Correctly rounded floating-point maths functions

Contents

1. Introduction

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 optimized 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 [60559:2020] mathematical operations; see 7.33.8 [Mathematics <math.h>]:

Function names that begin with cr_ are potentially reserved identifiers and may be added to the <math.h> header. The cr_ prefix is intended to indicate a correctly rounded version of the function.

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 roundTiesToEven rounding mode, as specified in [60559:2020]. The parameter type must satisfy the type trait std::numeric_­limits<T>::­is_iec559.

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.

1.1. Related papers

Proposal [P3375R3] 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 optimization 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 optimized 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 optimizations specified in [60559:2020].

2. Design considerations

2.1. Error handling

The error handling for the proposed functions should match that of the existing mathematical functions. That is:

• NaN is propagated.
• A range errors occur on overflow.
• Infinity may be propagated without error.
• A pole error occurs when dividing by zero.
• Underflow or denormalization is ignored.
• A domain error occurs when the result is not mathematically defined.

2.2. cr_sqrt for -0 arguments

If the argument to cr_sqrt is negative zero, the result should also be negative zero. This matches the behavior of sqrt in typical implementations (although this behavior is not mandated by the C standard), and it matches the behavior specified for squareRoot(-0) in [60559:2020].

3. Wording

The proposed changes are based on [N5014].

[version.syn]

In [version.syn], insert a feature-test macro as follows:

[cmath.syn]

In [cmath.syn], prior to the declaration of the mathematical special functions, insert the following:

Amend [cmath.syn] paragraph 1 to read:

Amend [cmath.syn] paragraph 2 to read:

[c.math.crfunc]

Immediately preceding [cf.cmath], introduce a new subclause:

4. References and bibliography