Fortran - Testing if a real is 'NaN'

Author

Testing if a real is 'NaN'

Martin Buchmann

2004-05-19, 11:32 am

Hi,

i'm having some problems finding out if a function returns a actual real
or something like NaN or INF/-INF.
Is there an easy way testing this?

Best regards
Martin

--
Never argue with an idiot. They pull you down to their level,
then beat you with experience.

B52B

2004-05-19, 11:32 am

Hi,

It is compiler dependent problem. Refer to the documentation
(topics: data types representation, bit patterns etc.)

Regards,
B52B

Martin Buchmann wrote:

> Hi,
>
> i'm having some problems finding out if a function returns a actual real
> or something like NaN or INF/-INF.
> Is there an easy way testing this?
>
> Best regards
> Martin
>

David Ham

2004-05-19, 11:32 am

If you have a compiler which implements the floating point technical
report (a sort of official extension to Fortran 95) then the standard
way to do it is as follows:

use ieee_arithmetic

.....

if (ieee_is_finite(x)) then

......

end if

ieee_is_finite returns true for any "normal" real and false for any sort
of inf or nan.

If you don't have a compiler which supports the TR then I would
recommend writing your own (incomplete) ieee_arithmetic module which
implements the ieee calls you need as wrappers around the non-standard
library calls of your compiler. That way, if your code needs porting
then at the worst all your non-standard ieee calls are in one place to
be replaced and if you are lucky, the next compiler will support the TR
and you can use the real ieee_arithmetic module.

The good news is that the TR has been incorporated in the f2003 draft.
The bad news is that in a fit of idiocy, the committee made it
optional. (Apologies to the committee - you do a great job but that call
was just daft).

David

On Wed, 19 May 2004 16:14:08 +0200
Martin_Buchmann@gmx.net (Martin Buchmann) wrote:

> Hi,
>
> i'm having some problems finding out if a function returns a actual
> real or something like NaN or INF/-INF.
> Is there an easy way testing this?
>
> Best regards
> Martin
>
> --
> Never argue with an idiot. They pull you down to their level,
> then beat you with experience.

Richard Maine

2004-05-19, 12:32 pm

Martin_Buchmann@gmx.net (Martin Buchmann) writes:

> i'm having some problems finding out if a function returns a actual real
> or something like NaN or INF/-INF.
> Is there an easy way testing this?

If your compiler is an f95 compiler that supports the IEEE TR, then see
the IEEE intrinsics. You could use either IEEE_CLASS or, likely more
convenient, IEEE_IS_FINITE.

Otherwise, as B52B says, it is compiler-dependent. Many compilers
have intrinsics for this. Or you could play games looking at the
bit patterns yourself, but that's a bit of a pain (and may end up
with dependencies on things like byte sex anyway).

I strongly recommend against some of the "tricks" like testing the
result of x.eq.x. As discussed here in the past, those tricks are
*NOT* robust. If you use them, you'll be fighting compiler
dependencies (and possible dependencies on particular compiler
options) anyway - better to identify and document the compiler
dependency than to just find something that appears to work in one
test and then gets burried in the middle of the code to become a bug
at some later time.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Richard Maine

2004-05-19, 12:32 pm

David Ham <d.a.ham@citg.tudelft.nl> writes:

> The good news is that the TR has been incorporated in the f2003 draft.
> The bad news is that in a fit of idiocy, the committee made it
> optional. (Apologies to the committee - you do a great job but that call
> was just daft).

We did? Oh yes, I guess so. No, I don't recall all the details of
the debate on why. Undoubtedly had something to do with some
architectures having various "difficulties" with it. Tales of
slowdowns like factors of 100 come to mind if you put an alpha in
the mode to do it all. No, I won't try to defend that number - I
just heard it mentioned. Anyway, trying to drag up memories from
about a decade ago, I think the optionality was likely a compromise
needed to avoid vehement objection from some vendors. Were it to
be done from scratch today, I'm not sure whether that would be
needed or not.

The other good news, though, is that I'm betting that all the f2003
compilers will implement it anyway, even before they implement all
of the mandatory parts of the f2003 standard.

--
Richard Maine | Good judgment comes from experience;
email: my first.last at org.domain | experience comes from bad judgment.
org: nasa, domain: gov | -- Mark Twain

Tim Prince

2004-05-19, 5:36 pm

"Martin Buchmann" <Martin_Buchmann@gmx.net> wrote in message
news:1ge1jko.cectxu14o7p68N%Martin_Buchmann@gmx.net...
> Hi,
>
> i'm having some problems finding out if a function returns a actual real
> or something like NaN or INF/-INF.
> Is there an easy way testing this?
>
If you want it portable and optimization-proof, use the IEEE TR mentioned in
other posts. If you can't do that, but are able to turn off optimizations
which violate IEEE,
if( x /= x) would detect NaN,
and various examples could detect INF, such as
temp = x * 0
if( temp /= temp .and. 1/x == 0)
If you have an f95 compiler, the SIGN() function should distinguish between
+-INF.
In gnu fortran, those should work as long as you don't set -ffast-math.
With other compilers, it's not often so easy.
In order to use ordinary optimizations in your code, you may need to perform
these IEEE-dependent tests in functions which are compiled separately, with
special de-optimization.

Robert Corbett

2004-05-19, 9:31 pm

David Ham <d.a.ham@citg.tudelft.nl> wrote in message news:<20040519165434.763038d2.d.a.ham@citg.tudelft.nl>...
>
> The good news is that the TR has been incorporated in the f2003 draft.
> The bad news is that in a fit of idiocy, the committee made it
> optional. (Apologies to the committee - you do a great job but that call
> was just daft).

Since there are still implementations of Fortran for machines that
implement floating-point arithmetic that is not compatible with
IEC 60559, the modules must be optional.

Sincerely,
Bob Corbett

Gerry Thomas

2004-05-19, 9:31 pm

Richard Maine

2004-05-19, 10:31 pm

robert.corbett@sun.com (Robert Corbett) writes:

> Since there are still implementations of Fortran for machines that
> implement floating-point arithmetic that is not compatible with
> IEC 60559, the modules must be optional.

Well, that doesn't quite follow as the modules don't require IEC 60559;
they are designed to be compatible with IEC 60559, but not to strictly
require it.

However, you may be right in the larger sense that there are probably
implementations of Fortran on machines that not only aren't compatible
with IEC 60559, but also aren't compatible with the less stringent
requirements of a minimal implementation of those modules. A minimal
one doesn't actually require a whole lot...but it does require something.
I don't recall the details without checking.

--
Richard Maine
email: my last name at domain
domain: sumertriangle dot net

Gerry Thomas

2004-05-20, 12:32 am

"Robert Corbett" <robert.corbett@sun.com> wrote in message
news:cb977dbc.0405191629.345adf2d@posting.google.com...
> David Ham <d.a.ham@citg.tudelft.nl> wrote in message
news:<20040519165434.763038d2.d.a.ham@citg.tudelft.nl>...
call[color=darkred]
>
> Since there are still implementations of Fortran for machines that
> implement floating-point arithmetic that is not compatible with
> IEC 60559, the modules must be optional.

It's time that these relics were exploded for the archeological posterity
of Fortran's past, present, and future, Marley in May no less, never mind.

--
E&OE

Ciao,
Gerry T.
______
"It is dangerous to be right in matters on which the established
authorities are wrong." -- Voltaire.

Eric K.

2004-05-20, 3:31 am

"Gerry Thomas" <gfthomas@sympatico.ca> wrote in message news:<MaTqc.68647$325.1332103@news20.bellglobal.com>...
> "Martin Buchmann" <Martin_Buchmann@gmx.net> wrote in message
> news:1ge1jko.cectxu14o7p68N%Martin_Buchmann@gmx.net...
>
> In Fortran, no. You could debug the code if the problem persists. Anyway,
> even if you happen on NaNs and Infs, then what?, never mind, Fortran has no
> error handling capabilities either.

Why the snotty attitude? Infs and NaNs don't always represent errors, and
even when they do, there's no portable and reliable method to trap or handle
them in C and C++, either.

In any case, Richard Maine is right: there's no foolproof way to check for
infinities and NaNs unless your compiler provides full IEEE support (most
don't). However, it's possible to exploit the IEEE representations for
floating-point values and the Fortran model for floating-point to obtain
a test that, in my own experience, isn't too likely to be foiled by
poorly conceived optimizers.

As I explained in the recent thread "Basic Real Exponent Question", the
IEEE single-precision floating-point representation has an exponent range
from -126 to +127. In the Fortran floating-point model, this maps to
minimum and maximum exponent values of -125 and +128, respectively. In
general, an exponent value of e in the IEEE model corresponds to an
exponent of e+1 in the Fortran model for base-2 arithmetic.

In IEEE arithmetic, +infinity corresponds to an exponent field of 255 (i.e.,
exponent equal to +128 after the bias of 127 is subtracted) and a zero
significand. A NaN has an exponent field of 255 (e=+128) and a nonzero
significand. An IEEE exponent value of +128 corresponds to a Fortran
model exponent of +129.

In the Fortran floating-point model, a zero IEEE significand corresponds
to a fractional part that is exactly equal to 1/2. So, to test for
an infinity or NaN, you first check whether the exponent field is larger
than the maximum value provided in the Fortran floating-point model.
Then you inspect the significand to see whether it's equal to 1/2 in the
Fortran model. The following little program illustrates the idea:

read(5,*) x,y
x=x/y
k=exponent(x)
if(k.gt.maxexponent(x)) then
s=fraction(x)
if(s.eq.0.5) then
write(6,*) 'infinity'
else
write(6,*) 'NaN'
endif
else
write(6,*) x
endif
stop
end

If you give this program inputs like x=1 and y=2, then the quotient is
finite and the program prints 0.5. If x=1 and y=0, then, on an
IEEE-compliant processor, x/y yields +infinity. If x=0 and y=0, then
IEEE defines the result to be NaN.

Hope this helps,
--Eric

Gerry Thomas

2004-05-20, 4:31 am

"Eric K." <ejko123@yahoo.com> wrote in message
news:4ccbcc1c.0405192231.2849ef6e@posting.google.com...
> "Gerry Thomas" <gfthomas@sympatico.ca> wrote in message
news:<MaTqc.68647$325.1332103@news20.bellglobal.com>...
real[color=darkred]
Anyway,[color=darkred]
has no[color=darkred]
>
> Why the snotty attitude?

Snotty is of your making.

> Infs and NaNs don't always represent errors,

No they don't, they're welcome info.

>and even when they do, there's no portable and reliable method to trap or
handle
> them in C and C++, either.
>

Well your superlatively wrong (but I'm sure you're used to that), there's
no way to do it in pure Fortran except via the 'system', aka C/C++.

> In any case, Richard Maine is right: there's no foolproof way to check
for
> infinities and NaNs unless your compiler provides full IEEE support (most
> don't). However, it's possible to exploit the IEEE representations for
> floating-point values and the Fortran model for floating-point to obtain
> a test that, in my own experience, isn't too likely to be foiled by
> poorly conceived optimizers.
>

You and Herr Maine are equally wrong. So go share a wombat or each other!

> As I explained in the recent thread "Basic Real Exponent Question", the
> IEEE single-precision floating-point representation has an exponent range
> from -126 to +127. In the Fortran floating-point model, this maps to
> minimum and maximum exponent values of -125 and +128, respectively. In
> general, an exponent value of e in the IEEE model corresponds to an
> exponent of e+1 in the Fortran model for base-2 arithmetic.
>
> In IEEE arithmetic, +infinity corresponds to an exponent field of 255
(i.e.,
> exponent equal to +128 after the bias of 127 is subtracted) and a zero
> significand. A NaN has an exponent field of 255 (e=+128) and a nonzero
> significand. An IEEE exponent value of +128 corresponds to a Fortran
> model exponent of +129.
>
> In the Fortran floating-point model, a zero IEEE significand corresponds
> to a fractional part that is exactly equal to 1/2. So, to test for
> an infinity or NaN, you first check whether the exponent field is larger
> than the maximum value provided in the Fortran floating-point model.
> Then you inspect the significand to see whether it's equal to 1/2 in the
> Fortran model. The following little program illustrates the idea:
>
> read(5,*) x,y
> x=x/y
> k=exponent(x)
> if(k.gt.maxexponent(x)) then
> s=fraction(x)
> if(s.eq.0.5) then
> write(6,*) 'infinity'
> else
> write(6,*) 'NaN'
> endif
> else
> write(6,*) x
> endif
> stop
> end
>
> If you give this program inputs like x=1 and y=2, then the quotient is
> finite and the program prints 0.5. If x=1 and y=0, then, on an
> IEEE-compliant processor, x/y yields +infinity. If x=0 and y=0, then
> IEEE defines the result to be NaN.
>
> Hope this helps,

No it doesn't.

Firstly, I don't stalk your every word but from now on I will, you being
such fun.
Second, whatever your point, you merely expose your naive perspective that
most here ought to have outgrown, although I have my doubts, never mind. So
grow up, learn something as you present as somewhat of an empty head, and
keep your juvenile crap to yourself.

--
E&OE

Ciao,
Gerry T.
______
"Ah, Klinger, my constant reminder that Darwin was right!." -- Maj Charles
Winchester III, the 4077th M*A*S*H

Martin Buchmann

2004-05-20, 8:32 am

David and all the others,

thanks for answering. I think at the moment my compiler is just to old
to implement the TR but i will check if more recent version of Absoft's
ProFortran have it. I'm planning to update anyway.

Futhermore, i guess the best way would be to check the calling arguments
of my function so that the whole thing is avoided at all. An idea which
came to my mind after posting. Some things are just too simple ;-)

Best regards
Martin

P.S.: Oh, Eric i'm sorry that my posting was the cause that you were
insulted by Mr. Thomas. I don't understand at all why some people abuse
usenet postings to find compensation for whatever went wrong in their
lifes.

--
"If Microsoft is ever going to produce something that does not suck,
it is very likely a vacuum cleaner."

Herman D. Knoble

2004-05-20, 9:34 am

Compilable, runable example code for detecting NaN's may be found at:
http://ftp.cac.psu.edu/pub/ger/fortran/hdk/nan.f90

Skip Knoble

On Wed, 19 May 2004 16:14:08 +0200, Martin_Buchmann@gmx.net (Martin Buchmann) wrote:

-|Hi,
-|
-|i'm having some problems finding out if a function returns a actual real
-|or something like NaN or INF/-INF.
-|Is there an easy way testing this?
-|
-|Best regards
-| Martin

Herman D. (Skip) Knoble, Research Associate
(a computing professional for 38 years)
Email: SkipKnobleLESS at SPAMpsu dot edu
Web: http://www.personal.psu.edu/hdk
Penn State Information Technology Services
Academic Services and Emerging Technologies
Graduate Education and Research Services
Penn State University
214C Computer Building
University Park, PA 16802-21013
Phone:+1 814 865-0818 Fax:+1 814 863-7049

Gerry Thomas

2004-05-20, 9:34 am

"Martin Buchmann" <Martin_Buchmann@gmx.net> wrote in message
news:1ge357j.128hjz21cpgosgN%Martin_Buchmann@gmx.net...

> David and all the others,
>
> thanks for answering. I think at the moment my compiler is just to old
> to implement the TR but i will check if more recent version of Absoft's
> ProFortran have it. I'm planning to update anyway.
>

You decrepit cheapskate, pay the wombat or get out of the sty.

> Futhermore, i guess the best way would be to check the calling arguments
> of my function so that the whole thing is avoided at all. An idea which
> came to my mind after posting. Some things are just too simple ;-)

Not really, but for remedial you, who's counting?, :-)

--
Ciao,
Gerry T.
_____________
"We want you!" -- Oberfeldwebel 'click your heels and stand erect' R.
Maine, Fortran recruitment.

Eric K.

2004-05-20, 8:31 pm

Martin_Buchmann@gmx.net (Martin Buchmann) wrote in message news:<1ge357j.128hjz21cpgosgN%Martin_Buchmann@gmx.net>...
> Futhermore, i guess the best way would be to check the calling arguments
> of my function so that the whole thing is avoided at all. An idea which
> came to my mind after posting. Some things are just too simple ;-)

If this is fast and easy to do, then it's certainly the best way to proceed.

For what it's worth, I'll post the basic fragment of my general function
for checking for IEEE infinities. The function is elemental, so you can
apply it to arrays (useful for checking for garbage input). Herman Knoble's
example code also will work, but the advantage of my approach is that
you don't have to worry about magic bit patterns, which can be problematic
when porting codes to machines with different endian conventions. (Endian
issues apply to the ordering of bits within bytes as well as to bytes within
words.) Also, my code is the same for both single and double precision;
just change the declaration of x. You can package this up into a module and
add a generic interface (email me for details if you don't know how to do this).

The following function returns one of 0, 1, or 2, according as its argument
is a normalized real, an infinity, or not-a-number, respectively.

elemental integer function typeof_real(x)
real,intent(in):: x
integer,parameter:: FINITE_REAL=0, INFINITE_REAL=1, NAN_REAL=2

if(exponent(x).gt.maxexponent(x)) then
typeof_real=merge(INFINITE_REAL, NAN_REAL, abs(fraction(x)).eq.0.5)
else
typeof_real=FINITE_REAL
endif
return
end function

You can extend this idea in a similar way to check for denormalized numbers,
too. The test is
if(exponent(x).lt.minexponent(x).and.abs(fraction(x)).ne.0.5) then
! x is denormalized
endif

Again, I make no claim that these functions are foolproof, but in my
experience, optimizers aren't likely to elide the tests and render them
useless. Note also that they apply *only* to IEEE single and double format.
The formats for quadruple-precision numbers are not (yet) standardized.
I'm not sure that the IEEE standard specifies how infinities and NaNs are
represented in the double-extended ("real*10") format, so these functions
may or may not work in that case (I haven't tested).

> P.S.: Oh, Eric i'm sorry that my posting was the cause that you were
> insulted by Mr. Thomas. I don't understand at all why some people abuse
> usenet postings to find compensation for whatever went wrong in their
> lifes.

Don't worry about it. In a technical forum like comp.lang.fortran, it's
pretty easy to spot the cranks and the fools. I just ignore them and advise
you to do the same.

--Eric

Martin Buchmann

2004-05-21, 6:31 am

Eric,

> The following function returns one of 0, 1, or 2, according as its argument
> is a normalized real, an infinity, or not-a-number, respectively.

thanks for the code. I will have a look at it.

Martin

--
"The Macintosh may only have 10% of the market,
but it is clearly the top 10%" - Douglas Adams 1952-2001

keith bierman

2004-05-21, 3:32 pm

Richard Maine wrote:
>
>
> I strongly recommend against some of the "tricks" like testing the
> result of x.eq.x. As discussed here in the past,

While you are, of course, correct .. I think there is a happy middle
ground. If you don't have a compiler which implements the TR, just code
up the routine you need (perhaps giving it the same name, or a name of
your own) and in the machine specific implementation you may well be
able to use things like x.ne.x ... you can compile the platform
specific parts with optimization disabled increasing the likihood that
you won't have to do assembly language coding for most platforms.

No doubt you were thinking of something like that ... but I thought it
might be worth making explicit.

David Ham

2004-05-24, 6:36 am

On Wed, 19 May 2004 18:03:21 -0700
Richard Maine <nospam@see.signature> wrote:

> robert.corbett@sun.com (Robert Corbett) writes:
>
>
> Well, that doesn't quite follow as the modules don't require IEC
> 60559; they are designed to be compatible with IEC 60559, but not to
> strictly require it.
>
> However, you may be right in the larger sense that there are probably
> implementations of Fortran on machines that not only aren't compatible
> with IEC 60559, but also aren't compatible with the less stringent
> requirements of a minimal implementation of those modules. A minimal
> one doesn't actually require a whole lot...but it does require
> something. I don't recall the details without checking.

As far as I can work out it doesn't require anything. The ieee functions
can be split into roughly three groups:

Inquiry - you need these but you can return "no support" in one way or
another for all of them.

Mode setting - these only need to exist or do anything useful (can't
quite work out which) if you support the relevant features.

Functions based on the values of variables - these are only required for
real kinds for which ieee_support_datatype(x) is true. You can therefore
dispose of them for all other kinds - which might be all kinds.

David

>
> --
> Richard Maine
> email: my last name at domain
> domain: sumertriangle dot net

Sponsored Links