Fortran Power Operator (**) Precedence Error

[K__Jin]

I am filing a potential bug about the Intel Fortran Compiler. When a negative value is used as a power, the power operation precedence might be wrong.

For example, I tested the following code on a Windows 10 x64 machine with Intel(R) Visual Fortran Compiler 17.0.6.270 [IA-32].

program testpower
implicit none
double precision :: s, a

s = 2**-0.5*2
write(*,*) '2**-0.5*2 = ', s

s = 2**(-0.5)*2
write(*,*) '2**(-0.5)*2 = ', s

end program

The program produces the following output (no errors or warnings during compilation):

 2**-0.5*2 =   0.500000000000000
 2**(-0.5)*2 =    1.41421341896057

For the first case, it gives a result as if it was 2**(-0.5*2).

This error happens if there is a "-" sign after the power operand. Or, maybe a parenthesis is required in this situation?

[K__Jin]

BTW, I also tested that s=2**+0.5*2 produces s=2.0. So, looks like we will have a problem if we have a sign after the **. I tested with gfotran and it gave an error message about missing parenthesis.

[JohnNichols]

program Test

    Implicit None

    Real(kind=8)      :: mu    = 1.68875694639187
    Real(kind=8)      :: sigma = 0.103093409616129
    Real(kind=8)      :: xi    = -1.07666780088544
    Real(kind=8)      :: x     = 1.76
    Real(kind=8)      :: tx, term_1, term_2, term_a, term_b

    call subroutine_test(mu, sigma, xi, x, tx)

    tx =           (1.0+xi*((x-mu)/sigma))**(-(1.0/xi))
    term_1 =       (1.0+xi*((x-mu)/sigma))**(-(1.0/xi))**(xi+1.0)
    term_a =                                         tx**(1.0+xi)
    term_2 = exp(-(1.0+xi*((x-mu)/sigma))**(-(1.0/xi)))
    term_b = exp(-tx)
    write (6,*) 'tx:     ', tx
    write (6,*) 'term_1: ', term_1
    write (6,*) 'term_a: ', term_a
    write (6,*) 'term_2: ', term_2
    write (6,*) 'term_b: ', term_b
    End

    subroutine subroutine_test(mu, sigma, xi, x, tx)
    Implicit None


    Real(kind=8), intent(in) :: mu
    Real(kind=8), intent(in) :: sigma
    Real(kind=8), intent(in) :: xi
    Real(kind=8), intent(in) :: x
    Real(kind=8), intent(inout) :: tx
    real(kind=8) :: a,b,c,d
 !    tx =           (1.0+xi*((x-mu)/sigma))**(-(1.0/xi))
    a = x-mu
    write(*,*) 'a: ', a
    b = a/sigma
    write(*,*) 'b: ', b
    c = -(1.0/xi)
    write(*,*) 'c: ', c
    d = 1.0+(xi*b)
    write(*,*) 'd: ', d

    tx =  d**c         
    write(*,*) 'tx: ', tx
    return
    end subroutine subroutine_test

Did you intend for terma and term1 to be different, because the only difference in the equations is a () in term1

It is to easy to make a mistake with long equations, unless you have an answer for a test sample then it can be blind flight.  

You know I had an old Physics Professor who had only two power points in his office, he plugged in a cable and ran a bare copper wire 240 volts set around the ceiling, it worked, but it was dangerous.  

My first Pure Math prof once did an inversion on the board, huge old fashioned boards that covered a huge wall.  I said to my friend, he made a mistake, friend said tell him, I looked at friend and said I did not have a death wish, in those days professors were gods. 

He finished did not get the answer, looked at it, I raised my hand and said, you made this mistake and explained the line of mistakes.  

I got a pass in the class.   But you calculate a more than once. 

In engineering, Mr Betts who was the best structural engineer ever, lost the entire class in every class, I loved it, at the end he had an eraser in his hand and if he spotted a mistake he started rubbing it out before I could point it out.   He gave me a distinction. I still cannot do the last question in the final, I can on a computer, but not by hand. I often thought it was a humility lesson.  

Of course there was the red headed kid in Physics 1, argued with the Professor, got 99 in the class where the average grade was 24.  

Simplify the code and you are less likely to make mistakes.  

Just an observation, it is your code.  

Using kind would be better. 

[JohnNichols]

I would be interested in @mecej4O  opinion.  

[TiborKibedi1]

Here is the definition of the PDF: 

Let us consider the case, when xi /= 0.   The expression of the PDF lends itself to calculate t(x) first, however turned out it is not the correct way in calculating the term in the middle:

t(x)**(xi+1.0)  

and 

(1.0+xi*(x-mu)/sigma)**(-(1.0/xi))**(xi+1.0)

give different results. The later one is the correct one. It is due to the way how consecutive operators with equal precedence (**, power) are evaluated. One possible way to break the expression of the PDF into shorter terms and getting the correct value is:

a = 1.0/sigma
b = 1.0+xi*(x-mu)/sigma
c = -(1.0/xi)
d = (xi+1.0)
e = exp(-(b**c))
PDF = a * b**c**d * e

Just looking the definition taken from the https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution page it might be easy to overlook, that inserting the expression  of t(x) into PDF, there are two consecutive power operators and breaking up the expression into shorter term needs care! Even inserting brackets could produce incorrect results.

Yes, read the manual first ....

[JohnNichols]

I think you found the answer. Good luck.  

How often is eta zero?

mathematics gives you if you are lucky answers, statistics just gives you more questions.  

This picture from your wiki entry is really nice, if you have say 3 million points you might get a curve like this, if you have ten points then you dream you have three million points.  

I have a machine we collect 10000 points per day, it is at least ten days till the t stat passes 2 and we can decide if it is Gaussian or non-Gaussian,  we have 8000 variables, maybe 20 are non-gaussian and each one will be a different distributions,  after 8 years of continuous data collection we know the answer on one bridge.  

Unless you have a lot of points, the stat is interesting it points you in a direction, but it is limited mainly by your imagination hoping one fits.  

These distributions are just mental gym, the data is never true to one stat. 

 

[Steve_Lionel]

D:\Projects>ifort /stand t.f90
Intel(R) Visual Fortran Intel(R) 64 Compiler for applications running on Intel(R) 64, 
Version 19.0.0.117 Build 20180804
Copyright (C) 1985-2018 Intel Corporation.  All rights reserved.

t.f90(5): warning #8810: Consecutive operators are an extension to the Fortran 2008 standard.
s = 2**-0.5*2
-------^

It's not a bug. The Intel documentation says:

Ordinarily, the exponentiation operator would be evaluated first in the following example. However, because Intel Fortran allows the combination of the exponentiation and minus operators, the exponentiation operator is not evaluated until the minus operator is evaluated:

A**-B*C is evaluated as A**(-(B*C))

[K__Jin]

Thanks Steve. I see the warning message after adding "/stand", and also find the explanation from Intel Fortran Reference. For completeness, I am just adding another example about consecutive operators (from Intel Fortran Reference):

When consecutive operators are used with constants, the unary plus or minus before the constant is treated the same as any other operator. This can produce unexpected results. In the following example, the multiplication operator is evaluated first, since it takes precedence over the minus operator:
X/-15.0*Y is evaluated as X/-(15.0*Y)

[mecej4]

It is unfortunate that this extension is enabled by default (i.e., if /stand is not specified). Consider this bizarre case:

program strange
real x,y,z
x=3.2
y=4.0
print *,x/+1.0*y,x/1.0*y
end

In the first expression, the mere presence of the normally superfluous but harmless '+' sign has the unfortunate effect of moving the y from the numerator to the denominator. The output of the program is

  0.8000000       12.80000

Even if you specify /stand, all that you are rewarded with is just a warning, which implies that we have to go hunting for a compiler option that has the effect of "when I say /stand, I really mean /stand -- do not produce an OBJ file for non-standard code". In a large program with many files, it is easy to lose sight of a warning if it comes in the midst of hundreds of other warnings.

[Steve_Lionel]

It's the general case of consecutive operators being allowed, otherwise you couldn't do something like a*-b. This extension goes back 40 years or more.

The option you want is /warn:stderrors

[K__Jin]

Thanks mecej4 and Steve for your input on this topic. I agree that the extension on consecutive operators is convenient in some cases. I also understand mecej4's point on this, because it may give unexpected results. Also, I do see lots of warnings in my program after using "/stand". I will just add parenthesis to avoid consecutive operations.

[mecej4]

Few compilers support this extension. I was surprised to see Gfortran doing so, but its interpretation of the semantics is quite different. For the program in #5, it says

S:\LANG>gfortran strange.f90
strange.f90:5:11:

 print *,x/+1.0*y,x/1.0*y
           1
Warning: Extension: Unary operator following arithmetic operator (use parentheses) at (1)

S:\LANG>a
   12.8000002       12.8000002

[K__Jin]

Seems like gfortran thinks unary operators have higher precedence than binary operators (at least for this case, it calculates the unary '-' operator, then division and multiplication from left to right) and I think this implementation (for this example) is more intuitive to me. But, how about -3**2? In this case, I am expecting a result of -9. I tested with gfortran and Intel, both gave -9 as expected. But, in this case if gfortran strictly obeyed the rule of unary operator has high precedence of binary operator, then we should get 9. But, looks like here gfortran thinks the "-" sign is not an unary operator but a binary minus operator and treat the expression as something like 0-3**2 which gives -9 (this is just my guess). So, my guess is that gfortran only interpret the "+/-" sign as an unary operator in the case of consecutive operation.

Sorry, in comment #2 I made a mistake by saying gfortran gave an error message. It was a warning message like mecej4 showed in the above comment. Thanks.

[mecej4]

Jin, Kai wrote:
... But, how about -3**2? In this case, I am expecting a result of -9. I tested with gfortran and Intel, both gave -9 as expected. But, in this case if gfortran strictly obeyed the rule of unary operator has high precedence of binary operator, then we should get 9.

This case does not require any extensions to standard Fortran. The standard precedence rules (Table 7.1 of Fortran 2008) tell us how to interpret the expression. In fact, the example in the 2008 standard that is given after that table is:

For example, in the expression
-A ** 2
the exponentiation operator (**) has precedence over the negation operator (-); therefore, the operands of
the exponentiation operator are combined to form an expression that is used as the operand of the negation
operator. The interpretation of the above expression is the same as the interpretation of the expression
- (A ** 2)

[JohnNichols]

If it is not intuitive and you have to grab the manual then it is to complicated. 

[TiborKibedi1]

I am using VS22 17.14.30 and Intel Fortran 2026.0. I run a similar problem as discussed a few years back in this thread. Run into trouble evaluating expressions involving exponentiation operator (**).  The following 3 lines are extractions from the attached code example:

tx = (1.0+xi*((x-mu)/sigma))**(-(1.0/xi))
term_1 = (1.0+xi*((x-mu)/sigma))**(-(1.0/xi))**(xi+1.0)
term_a = tx**(1.0+xi)

tx= 0.282048140135952. term_1=0.253991657520195 is the correct value, but term_a=1.10190080410454 is incorrect.

Using the /stand compiler switch does not report any missing bracket.

Any advice would be welcomed.

[JohnNichols]

 Implicit None

  Real(kind=8)      :: mu    = 1.68875694639187
  Real(kind=8)      :: sigma = 0.103093409616129
  Real(kind=8)      :: xi    = -1.07666780088544
  Real(kind=8)      :: x     = 1.76
  Real(kind=8)      :: tx, term_1, term_2, term_a, term_b
  
  tx =           (1.0+xi*((x-mu)/sigma))**(-(1.0/xi))
  term_1 =       (1.0+xi*((x-mu)/sigma))**(-(1.0/xi))**(xi+1.0)
  term_a =                                         tx**(1.0+xi)
  term_2 = exp(-(1.0+xi*((x-mu)/sigma))**(-(1.0/xi)))
  term_b = exp(-tx)
  write (6,*) 'tx:     ', tx
  write (6,*) 'term_1: ', term_1
  write (6,*) 'term_a: ', term_a
  write (6,*) 'term_2: ', term_2
  write (6,*) 'term_b: ', term_b
##   tx:       0.282048140135952
## term_1:   0.253991657520195
## term_a:    1.10190080410454
## term_2:   0.754237374569360
## term_b:   0.754237374569360
  Stop
  End

This type of code was ok in 1968 as they had limited memory, but today it is a recipe for a disaster. 

Break each micro step in the equations down so  a = x - mu

b = a/sigma and so on so that you can check the steps.  The point is to get the correct answer not the most obscure code that is hard to check.  Ok LISP people do this, but they are generally mathematicians who understand precedence and the use () a lot. 

plus xi  is negative so it makes it harder to look at it and make sense of it.  

[TiborKibedi1]

Thanks for the response. In my code I wanted to calculate the Probability Density Function, PDF of the so called Generalized Extreme Value distribution . I followed the description from the Wikipedia web cite: https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution.

The PDF is defined as:

PDF(x) = 1/sigma * t(x)**(xi+1) * exp(-t(x)) 

However evaluating t(x) first and then inserting to above equation gives the incorrect answer.

Combining all terms for xi/=0 case the PDF should be evaluated as 

PDF = (1.0/sigma) * (1.0+xi*((x-mu)/sigma))**(-(1.0/xi)) **(xi+1.0) * exp(-(1.0+xi*((x-mu)/sigma))**(-(1.0/xi)))

and this gives the correct value.

The source of the discrepancy is in the middle - equivalent to (t(x)**(xi+1) ). There are two consecutive power operators. The Intel Fortan compiler manual states:

Operators with equal precedence are evaluated in left-to-right order. However, exponentiation is evaluated from right to left. For example, A**B**C is evaluated as A**(B**C). B**C is evaluated first, then A is raised
to the resulting power.

The above long expression can be divided into sub-terms (micro steps), but with care!

In the calculation of the PDF and evaluating t(x) is equivalent to a**b, then using t(x) in the final expression of PDF(x) it would be equivalent to (a**b)**c. 

Finally working on complex calculations in nuclear physics, it is not uncommon to see very long equations, often spreading over a number of lines.