This is one of those cases where I wish the standard were a little more specific. As written, if the
scale() operation underflows, then the processor could initiate WWIII. I would like the consequences to be a little more confined than that. That clause could mean that the result is undefined in the same way that the equivalent multiplication operation
(X*B**I) is undefined when underflow occurs. That is, it could be set to 0.0 (as I would want in the case we are discussing), or it could signal an error. I would not want this operation to result in some random bit pattern with no warning, or start WWIII, or anything like that.
Here is a program that tests the original problem with all possible floating point numbers.
! compute the safe mean of two values that avoids overflow and unnecessary underflow.
real :: x, y, t, a
character(*), parameter :: cfmt='(*(a,1x,es15.8,1x))'
x = -huge(x)
t = nearest(x,1.0)
y = nearest(t,1.0)
a = safe_mean2(x,y)
if ( a .ne. t ) then
write(*,cfmt) 'error for x=', x, 't=', t, 'a=', a, ' y=', y
error stop 'error detected'
if ( y == huge(x) ) exit
x = t; t=y;
y = nearest(t,1.0)
write(*,cfmt) 'No errors detected'
function safe_mean2( x, y ) result(a)
real :: a
real, intent(in) :: x, y
integer :: k
intrinsic :: max, exponent, scale
k = max( exponent(x), exponent(y) )
a = scale( (scale(x,-k) + scale(y,-k)), k-1 )
end function safe_mean2
end program mean2
I have verified that this works for a few combinations of compilers and with various compiler options. However, with ifort it fails.
error for x= -1.17549449E-38 t= -1.17549435E-38 a= 0.00000000E+00 y= -1.17549421E-38
This is when
y=-tiny(x), and the problem does seem to trace back to corner cases of
scale() where the exponent changes from -125 to -126.
gfortran does not fail in this case, but it does seem to step through all of the denormal values between
-tiny() and 0.0, which I was not expecting. I was rather expecting
nearest(-tiny(x),1.0) to step over the denormals and go right to 0.0.
At least the undefined behavior did not start WWIII when I ran the program.
 A little experimentation reveals the following.
$ ifort -fp-model precise mean2.F90 && a.out
No errors detected
The problem without the option seems to be that the
nearest() function moves through the denormal numbers between -tiny(x) and +tiny(x), but the
scale() function flushes the denormals to zero. With the option,
nearest() does the same thing, but
scale() changes its behavior and produces denormal numbers (the correct ones, apparently).
This behavior of nearest() is a little surprising. I would have expected
nearest(-tiny(x),1.0) to be
0.0, and then I would expect
nearest(0.0,1.0) to be
tiny(x). That is, I would have expected
nearest() to skip over the denormal numbers.