minval/maxval with IEEE infinities
Thomas Koenig
real :: a(1)
a(1) = Infinity ! Assign IEEE Infinity, by whatever means
print *,minval(a)
Is it OK for an implementation to print the value of HUGE(a)?
OK to print "Infinity"?
OK to start WW III ? :-)
James Giles
I'd think the difficult case would be NaN and not infinity.
It seems to me that minval(a) where size(a) is one should be
the value of that one element. For multiple elements, you
return the one the compares less than or equal to all others.
Those two rules cover infinities (of either sign). But they
don't yield a clear result if the array contains NaNs. I believe
the answer should be NaN if the array contains any NaNs.
But, I think there is a standard somewhere that requires that
NaNs be disregarded in cases like this (you return the min
of the non-NaN values). In that case, I suppose you only
return NaN if all the elements of the array are NaN.
--
J. Giles
"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare
Dick Hendrickson
For that particular program, the answer is all three are OK. The
IEEE support in Fortran is funny. To get any of it, as far as the
language is concerned, you need to access one or more of the IEEE
support modules. Here, you really need
USE IEEE_ARITHMETIC
and the module has to have support for IEEE infinities. Then,
the answer might still be "who knows". Fortran doesn't really define
how IEEE arithmetic, etc., is supposed to work, that comes from the
IEEE standard and, as far as I know, the only intrinsic they formally
define is SQRT. So, I don't know what the IEEE standard says is
supposed to happen. Common sense says infinity should be returned,
but "common sense" is just a synonym for "processor dependent." ;)
Oddly enough, Fortran does say that if infinity is returned, the
processor must raise either OVERFLOW or DIVIDE_BY_ZERO.
Fortran explicitly says that the max/min-like functions must "in
effect" use the intrinsic relational operators for characters.
But, it is silent about methods for all of the other types. I
think this is force a collating sequence for character max/min
and at the same time, allow processors to use efficient
compare instructions for other types.
So, it looks like there is little in Fortran to answer your question.
It was a conscious decision by J3 not to specify how processors
support IEEE and not to specify what happens in cases where the
IEEE standard is silent.
Dick Hendrickson
Thomas Koenig
>Thomas Koenig wrote:
>> Assume a system with IEEE arithmetic, and
>>
>> real :: a(1)
>> a(1) = Infinity ! Assign IEEE Infinity, by whatever means
>> print *,minval(a)
>>
>> Is it OK for an implementation to print the value of HUGE(a)?
>> OK to print "Infinity"?
>> OK to start WW III ? :-)
>For that particular program, the answer is all three are OK.
Is this also the case if the processer specifies a(1) > huge(a(1)) ?
Depends on what you call the "minimum", maybe...
James Giles
The only description of the result of minval that applies to this
program is:
Case (i): The result of MINVAL (ARRAY) has a value equal to
the minimum value of all the elements of ARRAY if the size of
ARRAY is not zero. [...]
As far as I can tell, this requires the result to be equal to (at least)
one of the elements of the argument array. Hence neither WW III
nor A(1) (if HUGE(A(1) doesn't equal A(1)) qualify under this
requirement. Regardless of the "standard speak" explanations to
the contrary, if your implementation doesn't do that report it as
an error. I'm sure the vendor will agree with you that it's an error
and will fix it.
As I said before, the difficult case is if the array contains NaNs since
it doesn't compare equal to anything.
Dick Hendrickson
> Thomas Koenig wrote:
>> Dick Hendrickson <dick.hen...@att.net> wrote:
>>> Thomas Koenig wrote:
>>>> Assume a system with IEEE arithmetic, and
>>>>
>>>> real :: a(1)
>>>> a(1) = Infinity ! Assign IEEE Infinity, by whatever means
>>>> print *,minval(a)
>>>>
>>>> Is it OK for an implementation to print the value of HUGE(a)?
>>>> OK to print "Infinity"?
>>>> OK to start WW III ? :-)
>>> For that particular program, the answer is all three are OK.
>> Is this also the case if the processer specifies a(1) > huge(a(1)) ?
>> Depends on what you call the "minimum", maybe...
>
> The only description of the result of minval that applies to this
> program is:
>
> Case (i): The result of MINVAL (ARRAY) has a value equal to
> the minimum value of all the elements of ARRAY if the size of
> ARRAY is not zero. [...]
>
The intro to 13.7 (the function descriptions) says
"A program is prohibited from invoking an intrinsic procedure under
circumstances where a value to be returned in a subroutine argument
or function result is outside the range of values representable by
objects of the specified type and type parameters, unless the intrinsic
module IEEE ARITHMETIC (section 14) is accessible and there is support
for an infinite or a NaN result, as appropriate."
That propagates through all of the intrinsics. So, unless a processor
supports infinities, you can't use them.
It's my recollection that the IEEE standard does define the compare
order for the exceptional values. If I'm remembering correctly,
then a processor which has "support for an infinite result" must do
whatever the IEEE people say they should do. The standard specifically
doesn't say anything about how infinities and NaNs propagate
through computations. It relies on the IEEE standard for that.
Dick Hendrickson
FX
> yield a clear result if the array contains NaNs. I believe the answer
> should be NaN if the array contains any NaNs.
I'd think it is the most useful implementation.
--
FX
glen herrmannsfeldt
I believe in some statistics problems, and including the R
language, NaN is used for NA, that is, values for which no data
is available. In many cases the right thing is to just
ignore them.
-- glen
Richard Maine
But that's inconsistent with getting NaN as a result from invalid
operations such as 0./0. I have seen that kind of interpretation of NaN,
but I tend to think of it as attempting to shoehorn another meaning onto
NaN - one which doesn't quite fit.
At least if you always return a NaN, it will be obvious that something
is wrong. Quitely ignoring NaNs has the potential to be an unnoticed
error if that was not the intended meaning.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
glen herrmannsfeldt
> glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:
(snip)
>>I believe in some statistics problems, and including the R
>>language, NaN is used for NA, that is, values for which no data
>>is available. In many cases the right thing is to just
>>ignore them.
> But that's inconsistent with getting NaN as a result from invalid
> operations such as 0./0. I have seen that kind of interpretation of NaN,
> but I tend to think of it as attempting to shoehorn another meaning onto
> NaN - one which doesn't quite fit.
Well, R is mostly an interpreted language, and NA is much more common
in statistics problems than 0./0.. I believe, though, that it uses
a specific NaN (specific values of the lower bits) to distinguish
them when needed. Maybe just a meaning that IEEE didn't think about.
> At least if you always return a NaN, it will be obvious that something
> is wrong. Quitely ignoring NaNs has the potential to be an unnoticed
> error if that was not the intended meaning.
It might, for example, display a warning message on the screen before
ignoring it. Many R built in functions have options on what to do with
NA and NaN.
But even 0./0. needs to be ignored sometimes.
It is, for example, the mean of 0 numbers. I sometimes wonder about
this when I hear a baseball score in the middle of the first inning.
The home team hasn't even been up yet, so maybe the score should be
0 to NaN. (0 innings to score 0 runs.) Also, the batting average
of a player with 0 at bats.
-- glen
James Giles
I see Maine has already made the opposite case. My point
(in wording that has been deleted) was that NaNs are
problematical. For the majority of cases I'd prefer an
operation in which *any* operand was NaN to return NaN.
I concede there are also cases where that's a bad answer.
Whatever choice is made should be explicit though. Leaving
it "implementation defined" is the least desirable choice.
Of course, programmers can avoid the problem with
minval(array, mask=.not.isnan(array)))
Which is what you write if you want NaNs ignored. Or
you could write:
merge(NaN, minval(array), any(isnan(array)))
Which is what you write if you want NaN as the result
whenever any element is NaN.
(In standard Fortran 2003, ISNAN is spelled much more
verbosely and requires that you USE an intrinsic module.
I'd recommend renaming it to just ISNAN.)
And, of course, you can write a MYMINVAL function
whose definition is whichever of the above you prefer.
That makes calls more legible. Of couse, if the language
explicitly required the behavior of whichever of the above
was most commonly desired (no, I have no statistics), the
user would only have to redefine if (s)he wanted the opposite
behavior.
If this expression is part of a time-critical portion of you code
you might have to write an explicit FORALL or DO-loop to
avoid multiple passes throught the array. On the other hand,
it's plausible a good compiler could do the above (either form)
with one pass through the data already.
Tobias Burnus
> I'd think it is the most useful implementation.
This does not seem to be what IEEE 754 suggests.
At least in the draft at http://754r.ucbtest.org/drafts/754r.pdf
states:
"Operations on infinite operands are usually exact and therefore
signal no
exceptions."
"Every general-computational and quiet-computational operation
involving one or
more input NaNs, none of them signaling, shall signal no exception,
except
fusedMultiplyAdd (see 9.2) . For an operation with quiet NaN inputs
other than
max and min operations, if a floating-point result is to be delivered,
the
result shall be a quiet NaN, which should be one of the input NaNs."
"minNum(x,y) is the canonical floating-point number x if x < y, y if y
< x, the
canonicalized floatingpoint number if one operand is a floating-point
number
and the other a NaN. Otherwise it is either x or y.
maxNum(x,y) is the canonical floating-point number y if x < y, x if y
< x, the
canonicalized floatingpoint number if one operand is a floating-point
number
and the other a NaN. Otherwise it is either x or y."
If I understand this correctly:
maxval( [ 4.0, NAN ])
should return 4.0 while
maxval( [ NAN ])
should return NAN.
Tobias
James Giles
> On Feb 5, 11:01 pm, "FX" <coud...@alussinan.org> wrote:
>> I'd think it is the most useful implementation.
>
> This does not seem to be what IEEE 754 suggests.
> At least in the draft at http://754r.ucbtest.org/drafts/754r.pdf
> states:
[... quotations from the proposal ...]
> If I understand this correctly:
> maxval( [ 4.0, NAN ])
> should return 4.0 while
> maxval( [ NAN ])
> should return NAN.
I thought I'd heard that a standard somewhere required that.
It's not what I'd prefer. Still, having a standard meaning is
better than leaving it up to the implementation. It does mean
that for most of *my* uses I'd have to use an alternative:
merge (NaN, maxval(array), any(isnan(array)))
Of course, if this were in a speed-critical part of the program
I'd research ways to guarantee that the array contained no
NaNs in the first place.
Thomas Koenig
>If I understand this correctly:
> maxval( [ 4.0, NAN ])
>should return 4.0
I agree (for what it's worth).
>while
> maxval( [ NAN ])
>should return NAN.
Maybe it could also return -HUGE, I'm not sure.
How does the following sample implementation look of a 1-d minval
look?
real function real_masked_minval (a, n, mask)
integer, intent(in) :: n
real, intent(in) :: a(n)
logical, intent(in) :: mask(n)
logical :: seen_first
integer :: i,j
real :: result
seen_first = .false.
do i=1, n
if (mask(i)) then
result = a(i)
seen_first = .true.
exit
end if
end do
if (.not. seen_first) then
real_masked_minval = huge(a)
return
end if
do j=i+1, n
if (mask(j)) result = min(a(j), result)
end do
real_masked_minval = result
end function real_masked_minval