473,883 Members | 2,619 Online
Bytes | Software Development & Data Engineering Community
+ Post

Home Posts Topics Members FAQ

max(NaN,0) should be NaN

After tracking down a bug in my Fortran program, I found that it
assumed
max(NaN,0.) = 0.

This makes no sense, as the outcome of the operation is undefined and
should be NaN.
max(NaN,0.) = NaN

After researching, it appears the first outcome is accepted behavior,
and might be included in the revised IEEE 754 standard, which affects
not only Fortran. The discussion posted at
http://www.cs.berkeley.edu/~ejr/Proj...21.html#minmax
suggests that "There is no mathematical reason to prefer one reason to
another."

But I think otherwise, for the following reason. Suppose the NaN is
produced by x/y, where x=0 came from an underflow and y=0 came from an
underflow. Then x/y would be a well-defined number that could be
postive or negative. The convetion max(NaN,0.) = 0. is wrong at least
half the time.

Aug 28 '06
61 8691
In comp.lang.fortr an P.J. Plauger <pj*@dinkumware .comwrote:
(snip)
Well, the topic of this thread is one example. I assumed that NaN
always meant, "Abandon hope, all meaning is lost." The idea never
occurred to me that it might simply mean, "Pay no attention to me,
just go compute a useful value from other arguments." (And I still
have trouble with that viewpoint.)
The R language, mostly used for statistics work, has both NA and NaN.
NA when you don't know anything about the value, often used
in input data where no data exists. NaN usually as the result
of computations.

I believe the distiction could also make sense in other languages,
though I don't expect it to appear anytime soon.

-- glen

Aug 30 '06 #31
P.J. Plauger wrote:
"Nick Maclaren" <nm**@cus.cam.a c.ukwrote in message
news:ed******** **@gemini.csx.c am.ac.uk...
....
> or exceptions (as it specifies them)

I guess the parenthetic quibble gives you some wiggle room, but
I still have trouble with it. Our fenv.h tests seem to indicate
that the popular architectures *do* meet the requirements of
IEEE 754 in this area.
I suppose you could regard trap handlers as part of "IEEE exceptions".
I've never seen an implementation of IEEE traps.

--
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
Aug 31 '06 #32
P.J. Plauger wrote:
[...]

The code was also intentionally cavailer about preserving, and
generating, -0 results. (And I still think that -0 is of dubious
value.)
I hesitate to stick my head up in such august company, being so
nonaugust, but I was impressed by how easy it was to replicate
Kahan's Borda's mouthpiece example in C99-based ANS Forth
implementations (pfe and gforth) of complex mathematics:

http://www-personal.umich.edu/~willi...lex/borda.html

And furthermore how easy it was to get the right answers
generally for the complex elementary functions above and below
their principal value branch cuts:

http://www-personal.umich.edu/~willi...x-ieee-test.fs

The link just above is not comprehensible outside of a niche of
a niche. :-)

-- David
Aug 31 '06 #33

In article <cI************ *****@newssvr13 .news.prodigy.c om>,
"David N. Williams" <wi******@umich .eduwrites:
|P.J. Plauger wrote:
| >
| The code was also intentionally cavailer about preserving, and
| generating, -0 results. (And I still think that -0 is of dubious
| value.)
|>
|I hesitate to stick my head up in such august company, being so
|nonaugust, ...

Don't worry - as in politics, it is often the least august people who
see most clearly.

|but I was impressed by how easy it was to replicate
|Kahan's Borda's mouthpiece example in C99-based ANS Forth
|implementation s (pfe and gforth) of complex mathematics:
|>
|http://www-personal.umich.edu/~willi...lex/borda.html
|>
|And furthermore how easy it was to get the right answers
|generally for the complex elementary functions above and below
|their principal value branch cuts:
|>
|http://www-personal.umich.edu/~willi...x-ieee-test.fs

Well, yes, but one of the questions is whether that is useful. Now, I
quite agree that a numerical analyst as good as Kahan can use branch
cut information directly in real applications, and I date from the days
when we were taught to do so (though I doubt that I still can), but I
know of nobody under 45 who can.

My experience is that Bill Plauger is understating the case - the IEEE 754
handling of zeroes (and hence infinities) causes a huge number of errors
that would not have occurred in older arithmetics. People REALLY don't
expect the "gotchas" it introduces - and C99 makes that a lot worse.

|The link just above is not comprehensible outside of a niche of
|a niche. :-)

Very true. In a related niche, I have a simple test that shows up the
chaos with C99's complex division and overflow. It really is VERY
nasty :-(
Regards,
Nick Maclaren.
Aug 31 '06 #34

In article <Q9************ ******@bgtnsc04-news.ops.worldn et.att.net>,
"James Giles" <ja********@wor ldnet.att.netwr ites:
|>
|I suppose you could regard trap handlers as part of "IEEE exceptions".
|I've never seen an implementation of IEEE traps.

I have. I have even tested them. They didn't work, because few (if
any) modern systems support interrupt recovery by user code properly.
However, Plauger has mangled what I said so appallingly that I hope
that you don't think that I said what he implied I said.

What I said was:

Many or most don't support denormalised numbers or exceptions
(as it specifies them) in some of their modes (often their defaults),
and some older and embedded systems didn't/don't support NaNs or
infinities.

That statement is true and I stand by it. I did not say that any
modern general-purpose desktop and server systems don't support
IEEE 754 at all[*], though many do not support C99 - for good reasons,
including customer demand.

[*] zOS supports only a subset - see:

http://www-306.ibm.com/software/awdtools/czos/features/
Regards,
Nick Maclaren.
Aug 31 '06 #35

In article <hO************ *************** ***@giganews.co m>,
"P.J. Plauger" <pj*@dinkumware .comwrites:
|>
| Secondly, IEEE 754 has NOT gained widespread acceptance, and almost all
| "IEEE 754" systems use its format and not its semantics, or go in for
| some simplification or variant of it.
|>
|You could say the same thing about Standard C, if you're sufficiently
|vague about degree of conformance. The fact is that *most* architectures
|in widespread use support a pretty good approximation to IEEE 754. If
|you want to savor all the small ways they don't quite conform, then ask
|Fred Tydeman to give you his list of grievances. But in real life they
|hardly matter to the practicing programmer.

You cannot now say the same about C90, though you could up to about 1995.
You can, indeed, say the same about C99 - and a large part of the reason
for that is customer demand for C90 (for very good reasons).

Unfortunately, they do matter very much to anyone who is attempting to
write robust, portable code or is in the position of assisting users
who develop on one system and then find that their code doesn't work
on another. I have considerable experience with both.

| Many or most don't support
| denormalised numbers
|>
|I called you on this in Berlin and you temporized. The statement
|IME is simply untrue. And ME involves supporting Standard C99 and
|C++ libraries on about a dozen popular architectures, with about
|as many floating-point processors.

That is a falsehood. As you have done here, you quoted me out of
context, completely changing the meaning of what I said to something
that I did not say and was false, and I was not given a chance to
respond by the chairman. I wondered whether to object on a point
of order, but that seemed excessive.

| or exceptions (as it specifies them)
|>
|I guess the parenthetic quibble gives you some wiggle room, but
|I still have trouble with it. Our fenv.h tests seem to indicate
|that the popular architectures *do* meet the requirements of
|IEEE 754 in this area.
|>
| in some of
| their modes (often their defaults),
|>
|Another potential quibble, which I still have trouble believing,
|from direct experience.

For heaven's sake! Quoting individual phrases (not even clauses!) out
of context is a low political trick. I have requoted the full sentence
in another posting, and shall not do so here.

My statement was, however, based on detailed investigations of four
Unices on 4 totally separate architectures, and brief ones on a fair
number of others. In particular, it is true for at least AIX, IRIX
Solaris and Linux/gcc.

| Thirdly, only Java has attempted to include it as part of its
| specification,
| and Kahan has written a diatribe against Java.
|>
|IIRC, the diatribe did *not* complain that Java slavishly followed
|IEEE 754 to its detriment. There was this little thing about favoring
|Sun architecture over Intel, then patching thing up in a half-assed
|manner...

Yes, and there was a much larger thing about how it was essential to
provide both the flags and the values in order to get reliable exception
handling.

| C99 has some support, but
| words fail me when describing how unusable it is.
|>
|The detailed words always seem to fail you, but you're consistently
|quick with the poisonous adjectives.

As you know, that is another falsehood. I wrote a great many detailed
descriptions of the problems for the SC22WG14 reflector, often including
solutions, and some of them were raised by the BSI as official National
Body comments. All were either ignored or responded to entirely by
ad hominem attacks, as you are doing here.

I have a fair number of very detailed documents describing the issues,
and often solutions to the problems, which have been widely circulated.
I posted one to the SC22WG21 list, which you managed to get ignored
without consideration at a meeting at which I was not present, apparently
by claiming that it was false. However, I should point out that I gave
an independent reference that my statements were correct (Python).

[ To anyone else: please Email me if you want copies, and tell me
what aspects you are interested in. I don't guarantee to be able to
find everything! ]

|In fact, C99 has two extensive
|annexes (F and G) describing in exquisite detail how IEEE 754 and
|Standard C should play together.

Well, no, and that was one of the reasons that the BSI voted "no"
to C99 and many customers have explicitly specified C90. As I gave
enough reasons in a previous response in this thread, I shan't repeat
them.

|My current company Dinkumware, Ltd.
|has endeavored to conform completely to those annexes, and aside
|from a few aggressively perverse tests in Tydeman's floating-point
|suite we do just fine. While I don't always agree with some of the
|more creative choices made in those annexes (particularly G), I had
|little trouble following their guidance. Hence, I have to disagree
|that C99 support for IEEE 754 is unusable, either to us or our
|customers.

There is a difference between "unimplementabl e" and "unusable". I
never said that Annex F or G were unimplementable on systems with
hardware that supports IEEE 754 in at least one mode. They are
unusable for real applications that need robustness, portability
and efficiency (and, to some extent, any of those).

Enough, already. I am not going to respond to your attacks further.
If you ask me to justify what I have actually said, I may respond, but
that is unlikely if you misquote to the above extent.
Regards,
Nick Maclaren.
Aug 31 '06 #36
nm**@cus.cam.ac .uk (Nick Maclaren) writes:
[...]
Very true. In a related niche, I have a simple test that shows up the
chaos with C99's complex division and overflow. It really is VERY
nasty :-(
Can you post it, or a link to it if it's too large? I think a lot of
comp.lang.c folks would be interested (not sure about
comp.lang.fortr an).

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <* <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.
Aug 31 '06 #37
"Nick Maclaren" <nm**@cus.cam.a c.ukwrote in message
news:ed******** **@gemini.csx.c am.ac.uk...
...
|I called you on this in Berlin and you temporized. The statement
|IME is simply untrue. And ME involves supporting Standard C99 and
|C++ libraries on about a dozen popular architectures, with about
|as many floating-point processors.

That is a falsehood.
That's one.
...
As you have done here, you quoted me out of
context, completely changing the meaning of what I said to something
that I did not say and was false, and I was not given a chance to
respond by the chairman. I wondered whether to object on a point
of order, but that seemed excessive.
You're responding now, by questioning my integrity and motives.
...
|Another potential quibble, which I still have trouble believing,
|from direct experience.

For heaven's sake! Quoting individual phrases (not even clauses!) out
of context is a low political trick.
That's two.
...
| C99 has some support,
but
| words fail me when describing how unusable it is.
|>
|The detailed words always seem to fail you, but you're consistently
|quick with the poisonous adjectives.

As you know, that is another falsehood.
That's three strikes, and you're completely out of line. It's one
thing to accuse me of stating falsehoods, and of various other
tricks; it's quite another to say that I deliberately made a false
statement.
...
Enough, already. I am not going to respond to your attacks further.
If you ask me to justify what I have actually said, I may respond, but
that is unlikely if you misquote to the above extent.
Not to worry. The conversation is over.

P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com
Aug 31 '06 #38
Nick Maclaren wrote:
In article <cI************ *****@newssvr13 .news.prodigy.c om>,
"David N. Williams" <wi******@umich .eduwrites:
[...]
|but I was impressed by how easy it was to replicate
|Kahan's Borda's mouthpiece example in C99-based ANS Forth
|implementation s (pfe and gforth) of complex mathematics:
|>
|>
http://www-personal.umich.edu/~willi...lex/borda.html
|>
|And furthermore how easy it was to get the right answers
|generally for the complex elementary functions above and below
|their principal value branch cuts:
|>
|>
http://www-personal.umich.edu/~willi...x-ieee-test.fs
>
Well, yes, but one of the questions is whether that is useful. [...]
Indeed!
My experience is that Bill Plauger is understating the case - the
IEEE 754
handling of zeroes (and hence infinities) causes a huge number of errors
that would not have occurred in older arithmetics. People REALLY don't
expect the "gotchas" it introduces - and C99 makes that a lot worse.
I really don't have much experience here.
|The link just above is not comprehensible outside of a niche of
|a niche. :-)

Very true. In a related niche, I have a simple test that shows up the
chaos with C99's complex division and overflow. It really is VERY
nasty :-(
I'd be curious to see it. The stuff I was talking about doesn't
use the C99 complex math library at all, just the underlying
gcc support for IEEE 754.

-- David
Aug 31 '06 #39

In article <U1************ ****@newssvr25. news.prodigy.ne t>,
"David N. Williams" <wi******@umich .eduwrites:
| >
| Very true. In a related niche, I have a simple test that shows up the
| chaos with C99's complex division and overflow. It really is VERY
| nasty :-(
|>
|I'd be curious to see it. The stuff I was talking about doesn't
|use the C99 complex math library at all, just the underlying
|gcc support for IEEE 754.

Mine doesn't, either. Yes, it was designed with malice aforethought,
but I was checking whether the situation was really as bad as a look
at the code indicated it was. It shows what happens if you use the
'native' complex divide, the example in Annex G, or calculate the real
and imaginary parts separately. The results on several systems are
along the lines of:

0.00e+00 (1.000,1.000) (1.000,1.000) (1.000,1.000)
1.00e+307 (1.089,0.891) (1.089,0.891) (1.089,0.891)
2.00e+307 (1.154,0.769) (1.154,0.769) (1.154,0.769)
3.00e+307 (1.193,0.642) (1.193,0.642) (1.193,0.642)
4.00e+307 (1.207,0.517) (1.207,0.517) (1.207,0.517)
5.00e+307 (1.200,0.400) (1.200,0.400) (1.200,0.400)
6.00e+307 (1.176,0.294) (1.176,0.294) (1.176,0.294)
7.00e+307 (1.141,0.201) (inf,nan) (1.141,0.201)
8.00e+307 (inf,nan) (inf,nan) (inf,0.122)
9.00e+307 (nan,nan) (inf,nan) (nan,0.000)
1.00e+308 (nan,nan) (inf,nan) (nan,0.000)
1.10e+308 (nan,nan) (inf,nan) (nan,-0.000)
1.20e+308 (nan,nan) (inf,nan) (nan,-0.000)
1.30e+308 (0.000,-0.000) (inf,nan) (0.000,-0.000)
1.40e+308 (0.000,-0.000) (inf,nan) (0.000,-0.000)
1.50e+308 (0.000,-0.000) (inf,nan) (0.000,-0.000)
1.60e+308 (0.000,-0.000) (inf,nan) (0.000,-0.000)
1.70e+308 (0.000,-0.000) (nan,nan) (0.000,-0.000)
inf (nan,nan) (nan,nan) (0.000,-0.000)
inf (nan,nan) (nan,nan) (0.000,-0.000)

Note that the result should be decreasing, but becomes infinite before
it drops to zero. The example code does rather better, but the rule
that (inf,nan) is an infinity causes significant trouble - as I knew it
would, and was investigating. I don't know how much better is possible,
because this sort of issue is very nasty.
Regards,
Nick Maclaren.


#pragma STDC CX_LIMITED_RANG E OFF
#pragma STDC FP_CONTRACT OFF
#include <math.h>
#include <stdio.h>
#include <complex.h>

#ifndef TRUST_C99
#define creal(x) ((double)(x))
#define cimag(x) ((double)(-I*(x)))
#define INFINITY (1.0/0.0)
#define isfinite(x) (! isinf(x) && ! isnan(x))
extern double fmax(double,dou ble);
#endif

double complex cdivd(double complex z, double complex w)
{
double a, b, c, d, logbw, denom, x, y;
int ilogbw = 0;
a = creal(z); b = cimag(z);
c = creal(w); d = cimag(w);
logbw = logb(fmax(fabs( c), fabs(d)));
if (isfinite(logbw )) {
ilogbw = (int)logbw;
c = scalbn(c, -ilogbw); d = scalbn(d, -ilogbw);
}
denom = c * c + d * d;
x = scalbn((a * c + b * d) / denom, -ilogbw);
y = scalbn((b * c - a * d) / denom, -ilogbw);
if (isnan(x) && isnan(y)) {
if ((denom == 0.0) &&
(!isnan(a) || !isnan(b))) {
x = copysign(INFINI TY, c) * a;
y = copysign(INFINI TY, c) * b;
}
else if ((isinf(a) || isinf(b)) &&
isfinite(c) && isfinite(d)) {
a = copysign(isinf( a) ? 1.0 : 0.0, a);
b = copysign(isinf( b) ? 1.0 : 0.0, b);
x = INFINITY * ( a * c + b * d );
y = INFINITY * ( b * c - a * d );
}
else if (isinf(logbw) &&
isfinite(a) && isfinite(b)) {
c = copysign(isinf( c) ? 1.0 : 0.0, c);
d = copysign(isinf( d) ? 1.0 : 0.0, d);
x = 0.0 * ( a * c + b * d );
y = 0.0 * ( b * c - a * d );
}
}
return x + I * y;
}

int main() {
int i;
double d, r, z;
double complex c1, c2;

for (i = 0; i < 20; ++i) {
d = i*0.1e308;
c1 = (1.0e308+I*1.0e 308)/(1.0e308+I*d);
c2 = cdivd(1.0e308+I *1.0e308,1.0e30 8+I*d);
if (1.0e308 d) {
r = d/1.0e308;
z = 1.0e308+d*r;
printf("%.2e (%.3f,%.3f) (%.3f,%.3f) (%.3f,%.3f)\n", d,
creal(c1),cimag (c1),creal(c2), cimag(c2),
(1.0e308+1.0e30 8*r)/z,(1.0e308-1.0e308*r)/z);
} else {
r = 1.0e308/d;
z = d+1.0e308*r;
printf("%.2e (%.3f,%.3f) (%.3f,%.3f) (%.3f,%.3f)\n", d,
creal(c1),cimag (c1),creal(c2), cimag(c2),
(1.0e308*r+1.0e 308)/z,(1.0e308*r-1.0e308)/z);
}
}
return 0;
}
Aug 31 '06 #40

This thread has been closed and replies have been disabled. Please start a new discussion.

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.