fortran - Strange Function Call behavior

Question

I am quite new to Fortran and I was 'playing' with functions. I found a quite weird behavior in a very simple program that I cannot justify in any way. Here is the simple code:

Real*8 Function gamma(v)
 Implicit None
 Real*8 :: v
 gamma = 1.0_8 / Sqrt(1.0_8 - v ** 2)
End Function gamma

Program test_gamma
 Implicit None
 Real*8 :: v = 0.5_8, gamma
 print *,"Gamma = ", 1.0_8 / Sqrt(1.0_8 - v ** 2)
End Program

This prints the exact result: 1.1547005383792517 but if I use the function call, doing the same calculation, print *,"Gamma = ", gamma(v) I have the unexpected result 1.7724538509055161.

What am I missing here?

Answer 1

The value 1.7724538509055161 corresponds to the (correct) result from the (mathematical) gamma function with argument 0.5. The standard intrinsic function gamma returns result corresponding to this gamma function.

In your program calling gamma you have declared gamma to return a real*8 result, but you haven't given it the external attribute, or made an interface available through other means. So, the intrinsic is to called instead: when compiling the main program the compiler "knows" about no alternative [I'd expect a compiler warning in such circumstances if you ask for it.]

I'd recommend that you call your function something other than gamma rather than add the external attribute. I'd even recommend that for a function which has an interface available in the program.

By "interface available", I mean you could:

• put your gamma function internal to the main program;
• put your gamma function in a module used by the main program;
• [offer an interface block.]

Consult your programming guide for details.

Finally, I'd also recommend one other thing: no real*8. You can read many comments here on SO about that.