The C standard says in 6.5:

A floating expression may be contracted, that is, evaluated as
though it were a single operation, thereby omitting rounding errors
implied by the source code and the expression evaluation method.

with the following note:

The intermediate operations in the contracted expression are
evaluated as if to infinite range and precision, while the final
operation is rounded to the format determined by the expression
evaluation method. [...]

But it says nothing about casts, which are explicit conversion
operations that normally remove any intermediate extended precision,
as if they were assimilated to a real rounding operation (rather
than something that would just give an additional rounding error),
just like the rint() function.

So my question is:

In the following code

volatile double a = 0x1.00000001p+0;
volatile double b = -1.0;

printf ("%a\n", (double) (a * a) + b);
printf ("%a\n", (float) (a * a) + b);

may the (double) (a * a) + b and (float) (a * a) + b FP expressions
be contracted with the use of a FMA, i.e. as if fma(a,a,b) were used?

Based on my above remark, I would say "no".

But GCC (at least until a 14.0.0 20231017 snapshot), Clang (at least
until version 17) and Intel's oneAPI compiler 2021.3.0.20210619 all
contract (double) (a * a) + b with the use of a FMA; thus they output
0x1.000000008p-31 instead of 0x1p-31.

However, these compilers do not contract (float) (a * a) + b:
they output 0x1p-31. It is rather surprising that they treat it
differently.

Notes:
* With GCC, at least the -O2 optimization level is needed for
contraction of FP expressions, and options like -std=c17 must not
be used, because GCC disables the contraction of FP expressions
if such an option is present, due to the fact that GCC doesn't
support the STDC FP_CONTRACT pragma.
* A compiler option like -march=native may be needed in case the
default architecture is assumed not to have a hardware FMA.
* Intel's oneAPI compiler 2021.3.0.20210619 is actually buggy as
it contracts FP expressions even with -ffp-contract=off and the
STDC FP_CONTRACT pragma set to OFF.

--
Vincent Lefèvre <vincent@vinc17.net> - Web: <https://www.vinc17.net/>
100% accessible validated (X)HTML - Blog: <https://www.vinc17.net/blog/>
Work: CR INRIA - computer arithmetic / AriC project (LIP, ENS-Lyon)

Vincent Lefevre <vincent-news@vinc17.net> writes:

> The C standard says in 6.5:
>
> A floating expression may be contracted, that is, evaluated as
> though it were a single operation, thereby omitting rounding errors
> implied by the source code and the expression evaluation method.
>
> with the following note:
>
> The intermediate operations in the contracted expression are
> evaluated as if to infinite range and precision, while the final
> operation is rounded to the format determined by the expression
> evaluation method. [...]
>
> But it says nothing about casts, which are explicit conversion
> operations that normally remove any intermediate extended precision,
> as if they were assimilated to a real rounding operation (rather
> than something that would just give an additional rounding error),
> just like the rint() function.
>
> So my question is:
>
> In the following code
>
> volatile double a = 0x1.00000001p+0;
> volatile double b = -1.0;
>
> printf ("%a\n", (double) (a * a) + b);
> printf ("%a\n", (float) (a * a) + b);
>
> may the (double) (a * a) + b and (float) (a * a) + b FP expressions
> be contracted with the use of a FMA, i.e. as if fma(a,a,b) were used?
>
> Based on my above remark, I would say "no".
>
> But GCC (at least until a 14.0.0 20231017 snapshot), Clang (at least
> until version 17) and Intel's oneAPI compiler 2021.3.0.20210619 all
> contract (double) (a * a) + b with the use of a FMA; thus they output
> 0x1.000000008p-31 instead of 0x1p-31.
>
> However, these compilers do not contract (float) (a * a) + b:
> they output 0x1p-31. It is rather surprising that they treat it
> differently.

I can see how an argument might be made that contraction is
allowed for the (double) cast line, and not allowed for the
line with the (float) cast. A good test to do would be the
same exression but with a (long double) cast in place of the
(double) or (float) cast. I think an argument could be made
either way that contraction is allowed in that situation.

I think I'm mostly on your side, the cast not preventing the
contraction is unexpected and at least at first blush at
odds with what else the standard says. I completely agree
that the point needs clarifying.