Negative exponential should be compiler error - instead gives wrong answer

[caplanr]

Hi,

A colleague of mine noticed a bad behavior of IFX (and IFORT) for Fortran regarding negative exponents without parenthesis.  

From our understanding, the Fortran standard does not allow this:

program negative_exponent
  write(*,*) '10.0**-2 * 2=', 10.0**-2 * 2
end program negative_exponent

However, the above code compiles and runs using gfortran, nvfortran, ifort, and ifx.

On gfortran, it gives a warning about the syntax but then "does the right thing":

```
> gfortran negative_exponent.f90
negative_exponent.f90:2:37:

    2 |   write(*,*) '10.0**-2 * 2=', 10.0**-2 * 2
      |                                     1
Warning: Extension: Unary operator following arithmetic operator (use parentheses) at (1)
> a.out
 10.0**-2 * 2=   1.99999996E-02
```

Using NVIDIA's nvfortran, there is no warning but the code gives the right answer:

```
> nvfortran negative_exponent.f90
> a.out 
 10.0**-2 * 2=   2.0000000E-02
```

However, on ifort or ifx, there is no warning and gives the wrong answer:

```
> ifx negative_exponent.f90
> a.out
 10.0**-2 * 2=  9.9999997E-05
```

It seems to me that this code should ideally throw a compiler error since it is not allowed in Fortran (maybe have a flag to ignore the error?).

However, if it IS going to be allowed, it should give the right answer following standard math PEMDAS in compilation.

 - Ron

[mecej4O]

There is one more point of view that should be given some consideration as part of the decision to leave the compiler as is or to alter its default behavior. If the expression under question is not standard-conforming, but has existed in widely used code for decades, it is likely that the author(s) made reasonable attempts to make the code standard-compliant at the time of creation, given the widely used Fortran machine/compiler combinations available to them and it is equally likely that their colleagues exercised and tested the code segment several times and found it worked as they expected; what would be a reasonable guess for the semantics that they regarded as acceptable?

The authors' guide of the Physics Review Letters journal states:

 "Put in extra bracketing even where convention does not require it, if a likely misreading is thereby avoided. But leave them out where they would merely clutter the picture." 

That strikes me as good advice even for Fortran source code.

[Dave_Allured]

People need to be actively warned when they are using nonstandard ambiguous code that could result in hidden numerical errors with possible severe consequences.  You can not expect the entire user base to use the -strict flag or a code scanner.

I ask the Intel team to reconsider.  Here is a specific refinement of the above suggestion #2 from GaborToth.  Add two new flags to Intel fortran, plus default:

• (default) -- Consecutive operators generate a hard compile error, with informative message.
• -allow nonstandard_consecutive_operators -- Consecutive operators compile to Intel-specific interpretation, and informative warning is printed.
• -allow nonstandard_consecutive_operators_no_warning -- Consecutive operators compile to Intel-specific interpretation, and no warning is printed.

Documentation should be updated with an explicit warning that different compilers have completely different interpretations of this nonstandard construction, and will generate divergent numerical results.

This should provide the right degree of caution, while satisfying everyone's diverse needs.  Since the Intel compiler is already detecting this specific condition, it should be easy to convert existing behavior into this precautionary approach.

I am generally enthusiastic about legacy support extensions in compilers, but this one is dangerous and scary.  I would be glad to be warned that any of my legacy code is using a dangerous nonstandard construct, and to add a handful of extra parentheses.

[mecej4O]

It has been noted that reliance on mnemonics such as PEMDAS could create more problems than it solves (https://math.berkeley.edu/~gbergman/misc/numbers/ord_ops.html)

In addition, 'PEMDAS' is replaced by 'BOMDAS' or 'BIDMAS' in some English speaking countries such as Canada and New Zealand. I had never come across any such mnemonic as a high school pupil, even though the language of instruction was English. In German, the mnemonic is 'Punkt vor Strich'. Even in typeset mathematics, the risk of incorrect evaluation persists. Consider, for example, a/bc .The multiplication of b and c must be done first, and then a should be divided by the result.

[CarlBidwell]

Interesting catch, looks like the compiler should flag that. Maybe worth reporting as a bug for stricter standard checks.

[AlHill]

@CarlBidwell   Have you bothered reading this thread, or just spewing responses for some reason like all your other responses?

Doc (not an Intel employee or contractor)
[AI is the same as snake-oil]

[CarlBidwell]

I get where you’re coming from, I did read through the thread, and I’m just trying to add something useful rather than repeat what’s already been said. Definitely open to corrections if I’ve missed a key point.

[mecej4O]

It has been noted that reliance on mnemonics such as PEMDAS could create more problems than it solves (https://math.berkeley.edu/~gbergman/misc/numbers/ord_ops.html)

In addition, 'PEMDAS' is replaced by 'BOMDAS' or 'BIDMAS' in some English speaking countries such as Canada and New Zealand. I had never come across any such mnemonic as a high school pupil, even though the language of instruction was English. In German, the mnemonic is Punkt vor Strich, see Punktrechnung vor Strichrechnung

[mecej4O]

One way to avoid the ambiguities in the test code is to use named constants, using the PARAMETER attribute:

real,parameter :: base=10.0, expnt=-2.0, factor=2.0
print *, *factor*(base**expnt)

[Steve_Lionel]

Your friend is mistaken. Consecutive operators were disallowed in Fortran 77. Here's a quote from that standard:

Note that these formation rules do not permit
expressions containing two consecutive arithmetic
operators, such as A**-B or A+-B. However, expressions
such as A**(-B) and A+(-B) are permitted.
(FORTRAN 77 Full Language, 6.1.2.4 Arithmetic Expression)

It wasn't Intel that created this extension, and I doubt that DEC created it either. It probably came from IBM in the 1970s.

As I explained earlier, Fortran does not treat "-2" as a single entity, whereas the other languages you mention do. Intel Fortran, like its predecessors back into the 1970s, is simply applying Fortran standard precedence rules. If you don't like it, use parentheses the way the standard requires you to do.

[GaborToth]

Indeed, I looked up the actual Fortran 77 standard https://nvlpubs.nist.gov/nistpubs/Legacy/FIPS/fipspub69-1.pdf, and it does not allow it. But many of the specific implementations do, including this one: FORTRAN 77 Language Reference, FORTRAN 77 Version 5.0 (by Sun corporation). Sounds pretty official, right? And it does say the same thing as the Intel interpretation. Since this was _never_ part of the standard, insisting on this odd behavior is even less justified.

Note:

perl -e 'print -4**-0.5*3," ",-(4**(-0.5))*3," ",(-4)**(-0.5)*3'
-1.5 -1.5 NaN

so -4 is _not_ interpreted as an "entity". The unary - operator is evaluated _after_ the **, but _before_ the * and /, which is exactly what mathematics demands. That's the proper rule, and the idea that left-to-right and some obscure ordering of operators is more important than mathematics conventions is just that, an idea. And not a good one. 

[PeteRiley]

I'd like to push back on this line of argument (not directed at you in particular, Steve, but the conversation has been meandering off course in my opinion). This really isn’t about who is “right” or “wrong” in interpreting a decades-old extension; it's a much more pragmatic question, at least to me - what best serves users today? The bottom line is that the current behaviour silently produces divergent numerical results across compilers, which is a serious risk in scientific and engineering applications.

The issue has been well-documented here, and my superficial searches have found a number of other cases, so I think it is important enough to be addressed and not just dismissed. If Intel won't address it, then they should defend that decision adequately. And I think that, at least in part, that should include why choosing not to address it serves their users better. Thus far, the 'reasons' have included: it's just a "corner case", referencing an arcane historical point, or making it a user problem (to require the use of parens). Now, these might get you to "I'm right", but they sound more like excuses to not have to deal with the issue. Most importantly, they don't address the question of whether it serves the user better.

So I'd support the proposal (voiced several times in this conversation) that the practical path forward is to make ambiguous consecutive-operator forms an error by default, while providing a clear opt-in flag for legacy code. That way, safety and portability become the default, but continuity is still preserved for those who need it. Do you think that the majority of users would support or oppose this change? 

[jimdempseyatthecove]

 FORTRAN 77 Version 5.0 (by Sun corporation). Sounds pretty official, right?

No, what came out of the standards committee is official.

Jim Dempsey

[jimdempseyatthecove]

>>make ambiguous consecutive-operator forms an error by default, while providing a clear opt-in flag for legacy code. 

I agree with this, but with an additional feature.

This would be following the error message (sans opt-in flag), would be how the opt-in would interpret the statement.

Jim Dempsey

[Steve_Lionel]

I'd be in favor of making it an error by default with an option to allow it. What I would not support is silently changing the documented behavior.

[Arjen_Markus]

I would like to point out another peculiarity of the minus sign to the already large mix:

In languages that consider a minus sign as part of the number, what would they do with expressions like: -3.0 ** x? Surely, this would mean - if they are consistent - (-3.0) ** x? But that would make the interpretation of -y ** x different from a case with a signed constant. Would they have to use an exception and parse the first expression as -(3.0) ** x? (Languages like C solve it by not having a power operation, I suppose ;)).

In view of this sort of ambiguities, the choice Fortran made to NOT consider the minus sign as part of the number may not be that bad after all.

[mecej4O]

Arjen, reading your post makes me wish (once more) for a handy card with a BNF or EBNF (Backus-Naur or Extended Backus-Naur) formulation of the Fortran syntax rules. I remember reading the rules for Algol and Pascal in those forms, and they make things pleasantly clear!

I have to wonder what Edsger W. Dijkstra, the renowned Dutch computer scientist, would have said after reading the posts in this thread.

For a startling example of why we need clear definitions and rules, try to pin down the exact year in which Queen Elizabeth I of England died -- 1602 or 1603?

[GaborToth]

Now that I advanced to the novice status, this is how pI is defined in the SWMF/share/Library/src/ModNumConst.f90  (note that the name of the module is not "Base", which is not very meaningful, but ModNumConst that describes what it contains (we also have a ModConst module that contains physics constants. So here we go:

real, parameter:: cPi = 3.1415926535897932384626433832795

real(Real8_), parameter:: cPi8=3.1415926535897932384626433832795

and Real8_ is defined in SWMF/share/Library/src/ModKind.f90 as

integer, parameter :: Real8_ = selected_real_kind(12,100)

The first cPi definition works either for compilation with default precision, or compilation with -r8 (or similar) flag, the second cPi8 insists on double precision. There is no function call, so this is way simpler to use and implement than the pi() function. The number of digits exceeds the double precision, so there is no issue with accuracy either. 

Defining 1e-6 as Var3 = 1/100000.0 does not achieve anything useful. If double precision is the goal, use 1d-6. 

Best regards,

Novice

[mecej4O]

You provide a string of 30 digits for Pi, but that is going to have a precision of only 7 decimal digits. To see the difference, add "d0" to the end of the string and compare the results.

[GVautier]

But with that value does cos(cpi8) returns -1. ?

[andrew_4619]

program pi_test
    implicit none
    real(8), parameter:: cPi8 = 3.1415926535897932384626433832795_8
    real(8), parameter:: myPi = acos(-1.0_8)
    real(8) :: vv1, vv2
    vv1 = cos(cPi8)
    vv2 = cos(myPi)
    print *, 'cos(cPi8)', vv1
    print *, 'cos(myPi)', vv2
    print *, vv1-vv2
    print *, cPi8-myPi
end program pi_test

cos(cPi8) -1.00000000000000
cos(myPi) -1.00000000000000
0.000000000000000E+000
0.000000000000000E+000

Of course myPi is less typing and less human error prone so wins every time

[GaborToth]

OK, there is always something new to learn. Since Fortran 2003, acos(-1.0) can be put into the parameter declaration. Fortran 95 still did not allow that. When we started the SWMF in 2003, there were no working Fortran 2003 compilers. So we did not have this option, but now it is certainly a good solution. Note that I was only comparing setting a parameter with the function pi(), which I find an overcomplicated approach.

On the other hand, real(8) is not correct, because 8 does not mean 8-byte real in general. The only general approach I know of is 

integer, parameter :: Real8_ = selected_real_kind(12,100)

and then "Real8_" can be used instead of "8". An alternative is "double precision", but that used to mean something else (16-byte) on a Cray vs. all other computers (Cray still did this in 2003). I guess these days "double precision" means 8-byte real for all compilers. But I find typing "double precision" tiresome, and also it does not allow changing the precision with a compiler flag. So we use "real" and get the compiler flag to change the precision. 

Here comes the interesting question: when using a constant like 3.1415926535897932384626433832795 to set a double precision constant, is it  4-byte or 8-byte precision?

gfortran seems to read all the digits, but nagfor does not, at least with the default settings. When the compiler is used with the double precision compiler flags, all reals are upgraded to double precision, including constants. So our cPi8 definition may be wrong when we compile in single precision (which we don't do much, and cPi8 is apparently not used anywhere in our code, whew). But the correct definition would indeed be 

real(Real8_), parameter:: cPi8 = 3.1415926535897932384626433832795d0

where the "d0" takes care of making the number double precision, or using Fortran 2003 

real(Real8_), parameter:: cPi8 = acos(-1d0)

Our cPi definition is OK, but it could be replaced with

real, parameter:: cPi = acos(-1.0)

which is correct for both single and double precision compilation. In fact, I replaced our definitions with this. Progress is made one commit at a time.