correctly port FLT_EVAL_METHOD

[vks]

FLT_EVAL_METHOD is defined by the C99 standard and determines the precision of intermediate results.

Currently we assume FLT_EVAL_METHOD = 0. Do we want to support this feature correctly?

[hanna-kruppe]

What does "supporting the feature correctly" even mean? Do the MUSL algorithms inspect FLT_EVAL_METHOD to compensate for possible differences in rounding introduced by methods other than 0, or do they exploit the less strict methods to give faster but less-deterministically-rounded results themselves?

[vks]

I don't think it is the latter, because the default (FLT_EVAL_METHOD = 0) is the behavior you would expect. Increasing FLT_EVAL_METHOD should result in more correctly rounded results. From Wikipedia:

Expression evaluation is defined to be performed in one of three well-defined methods, indicating whether floating point variables are first promoted to a more precise format in expressions: FLT_EVAL_METHOD == 2 indicates that all internal intermediate computations are performed by default at high precision (long double) where available (e.g., 80 bit double extended), FLT_EVAL_METHOD == 1 performs all internal intermediate expressions in double precision (unless an operand is long double), while FLT_EVAL_METHOD == 0 specifies each operation is evaluated only at the precision of the widest operand of each operator. The intermediate result type for operands of a given precision are summarized in the table on the right.

FLT_EVAL_METHOD	float	double	long double
0	float	double	long double
1	double	double	long double
2	long double	long double	long double

[hanna-kruppe]

Increasing FLT_EVAL_METHOD should result in more correctly rounded results.

It usually will, on average, but not always, especially since inconsistent rounding of intermediate results can have far-reaching consequences. As a contrived example, consider:

float z = x * x;
float L2 = sqrt(z - y * y);

Under FLT_EVAL_METHOD == 0, L2 will be a non-negative number if x and y are finite and squaring them does not overflow, and in particular it will be zero if x and y are also equal. On the other hand, if intermediate results have higher precision, the subtraction can yield a (very small) negative number and thus make L2 a NaN. Specifically, if x == y and x * x is rounded down, y * y is not similarly rounded before being subtracted from z, so (in the higher intermediate precision) z < (y * y) despite z == x * x and x == y.

---

But this is actually an argument for the first option (libm trying to give correctly rounded results even in the face of inconsistent rounding). After doing some reading, I'm not even sure whether FLT_EVAL_METHOD is even allowed to affect the results of math functions 🤷‍♂️

I searched and found commit messages such as the following, which do fall in the first category:

• https://git.musl-libc.org/cgit/musl/commit/src/math?id=79ca86094d70f43252b683c3a3ccb572d462cf28
• https://git.musl-libc.org/cgit/musl/commit/src/math?id=4b539a826b64d846d02a6cde9d6bcff22472af88
• https://git.musl-libc.org/cgit/musl/commit/src/math?id=8eead3ef18ea71a64ef3cbff8c09bac3b82f1242

As also noted in the last one, the only relevant architecture where FLT_EVAL_METHOD != 0 is useful is i368 (i.e., x87 FPU and no SSE). rustc does support targeting that, but it's not typically a no_std target so it seems unlikely anyone will want to use this crate on that target.