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.

```
program mean2
! compute the safe mean of two values that avoids overflow and unnecessary underflow.
implicit none
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)
do
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'
endif
if ( y == huge(x) ) exit
x = t; t=y;
y = nearest(t,1.0)
enddo
write(*,cfmt) 'No errors detected'
contains
function safe_mean2( x, y ) result(a)
implicit none
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 )
return
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
error detected
```

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.

[edit] 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.