• casts and contraction of floating expressions

    From Vincent Lefevre@vincent-news@vinc17.net to comp.std.c on Wed Nov 22 16:35:57 2023
    From Newsgroup: comp.std.c

    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)
    --- Synchronet 3.20a-Linux NewsLink 1.114
  • From Tim Rentsch@tr.17687@z991.linuxsc.com to comp.std.c on Mon Nov 27 11:11:27 2023
    From Newsgroup: comp.std.c

    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.
    --- Synchronet 3.20a-Linux NewsLink 1.114