Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Posting-Version: version B 2.10.1 6/24/83; site alice.UUCP
Path: utzoo!watmath!clyde!burl!ulysses!allegra!alice!allan
From: allan@alice.UucP (Allan Wilks)
Newsgroups: net.math.stat
Subject: S and NaN on SUNs
Message-ID: <4653@alice.UUCP>
Date: Wed, 4-Dec-85 20:14:28 EST
Article-I.D.: alice.4653
Posted: Wed Dec  4 20:14:28 1985
Date-Received: Thu, 5-Dec-85 09:06:46 EST
Organization: Bell Labs, Murray Hill
Lines: 52

There is a problem using the S statistical system with IEEE floating
point arithmetic on a SUN workstation with release 2.0 of the SUN
software.  People with other kinds of machines need not be concerned.

The problem was pointed out by Steve Peters at MIT:
computations that in the old days would fail due to
a floating point exception will now sometimes silently return NaN.
This is fine, provided that in succeeding computations, S is able
to recognize those NaNs.  It fails to do so, however, because of the
way NAs are implemented in S, viz, a single, hardwired 32-bit pattern is chosen 
for representing NA, and the choice is made in such a way that the pattern
is unlikely to appear either as an integer or as a (single-precision)
float.  All tests for NA go through the single routine $P/natst.r
which simply tests its input (32-bit) argument against the fixed NA
pattern.  This will not work for the current problem as there are
many possible IEEE NANs.  So here is $P/natst.C to replace $P/natst.r:

===============================================
SUPPORT(natst,		test for missing value)
/*
 * Version of natst for IEEE floating point.
 * Tests whether the bit pattern of the input
 * matches one of the bit patterns for NaN,
 * while avoiding a match with small integers.
 */
F77_SUB(natst)(i)
long *i;
{
	int a, b;

	a = (*i & 0x7f000000) == 0x7f000000;
	b = (*i & 0xff800000) == 0xff800000;
	return(a & !b);
}
===============================================

It is a little tricky because many NaNs look like small negative
integers when treated as integers, but this routine should work
pretty reliably.  To use it, install it in $P, move natst.r to
some other file, and as user s do:
	cd $P
	S MAKE natst.o
	ar r $L/psl natst.o
	ranlib $L/psl
To remake an S executive:
	cd $M
	S MAKE NEW.S INSTALL
External functions should also be remade.
Send me mail if you have problems, or if you can think of a better
way to cut through this murk.  Also, please note that this doesn't
seem to be a problem with S under other IEEE implementations.