Thanks for the clarification. I see now that implicit is devoted to types and would be a poor choice.
Ideally we would want something that can be supported by the language standard, so if there is a standards-compliant pragma-like solution then I think that would work well for us - especially if it could be defined at the module level.
Malcolm noticed something I had missed, which saddens me - in that the standard talks about mathematical operations but also notes that computational operations can differ. I agree with Malcolm’s comments. I’m not going to try to discuss the IEEE rules, as this is not a topic I’m well-versed in.
The standard really should stay out of trying to constrain how compilers optimize. The compilers I am familiar with have options to restrict optimizations that can lead to computational differences (FMA, etc.)
I must admit that I have never fully understood that distinction. I understand it when it comes to integer vs. real operations (1/3 is different from 1./3.), but not when it comes to real operations of different kinds, or with extended register intermediate results, or with different hardware (scalar, vector, gpu). I think of all these in a mathematical way, so at that point it is all mathematical to me, but the intention of the standard is to somehow separate some of these cases to allow compilers flexibility, while restricting others. This goes back even to the f77 era, so my general confusion is not new.
I remember in f77, before the epsilon() intrinsic was standard, that programmers would use expressions to evaluate that quantity. One of them was something like
eps = (1./3.)*3. -1.
That works if the compiler doesn’t do anything tricky because of the binary representation of (1./3.). But if the compiler does anything at all to make the statement “more precise”, like extended precision intermediates, or guard bits, etc., then it can fail with a 0.0 result. And by the 1980s, compilers often did exactly those kinds of things.
Algebraic identities aren’t, when operations are implemented in floating-point. One (noted in the Standard) purpose of parentheses is to stop the compiler applying algebraic identities.
Not in our case, but this is primarily because we try to control the order of arithmetic in every expression.
Disabled FMAs (e.g. bit-equivalence after a 90° rotation) are not documented, but I suppose we consider the parentheses to be a light form of documentation. In truth, the FMAs that we accept are probably in equal numbers.
If a solver happens to explicitly require order of operations (e.g. maintain montonicity, avoid division-by-zero) then there is typically a comment. But they don’t follow a standard syntax.
I could message you with some specific examples if you think it would be helpful.
In **Fortran**, parentheses have no numerical meaning other than syntactic grouping. `(A*B)+C` only specifies an existing precedence; it is not a contract pertaining to roundings. Unless you specifically opt out using ieee semantics or divide the expression between statements, the processor may still fuse, reassociate, vectorize, or do anything else since the standard regards `A*B+C` and `(A*B)+C` as one expression tree.
When they alter the tree, such as in `(A+B)+C` vs. `A+(B+C)`, parenthesis *do* matter. Since addition isn’t associative in finite precision, that’s a separate matter. The additional parens there are decorative or for human readers, not a numerical barrier, since `*` binding is already tighter than `+`.
ieee_arithmetic, ieee_real, or a temporary punctuation is not enough to really ban fma or need two roundings. Otherwise, the quickest valid evaluation may be chosen by the compiler, which is exactly what contemporary Fortran compilers will do.
I disagree about IEEE_REAL, Helen. The specification for it reads
The result has the same value as A if that value is representable in the representation method of the result, and is rounded according to the rounding mode otherwise. This shall be consistent with the specification of ISO/IEC 60559:2020 for the convertFromInt operation when A is of type integer, and with the convertFormat operation otherwise.