From Martin Maechler <maechler@stat.math.ethz.ch>  Wed Jan  7 08:00:24 1998
From: Martin Maechler <maechler@stat.math.ethz.ch> (Martin Maechler)
Date: Wed, 7 Jan 1998 09:00:24 +0100
Subject: IEEE 754 Style Arithmetic
In-Reply-To: <199801062351.MAA09299@stat1.stat.auckland.ac.nz> (message from
 Ross Ihaka on Wed, 7 Jan 1998 12:51:15 +1300 (NZDT))
Message-ID: <199801070800.JAA03967@sophie.ethz.ch>

>>>>> "Ross" == Ross Ihaka <ihaka@stat.auckland.ac.nz> writes:

    Ross> I have been looking at the R mathematical library with a view to
    Ross> making changes so that it will handle IEEE 754 entities like NaN
    Ross> and +/- Inf.  This appears to be not too hard and I am fairly
    Ross> well down the path to converting the existing code (and
    Ross> simultaneously converting some of the more suspect algorithms to
    Ross> something more solid).
Great !  
Thank you -- this has been my longest lasting wish for R...

(has Heiner Schwarte contacted you, yet?
 during holiday season, he has almost finished doing the same
 while successfully porting R 0.61 to Windows NT (!!!) (MS C compiler)
)


    Ross> Q1:

    Ross> I have looked at the Splus implementation and I have a bit of a
    Ross> problem with some of the decisions that have been made.  For
    Ross> example

    Ross> 	pnorm(0, 1, Inf) -> 0.5

    Ross> This seems VERY dangerous to me.  It's the kind of thing which
    Ross> could make errors almost impossible to find.  I would much rather
    Ross> that this generated NaN (it's a domain error in my view).

I disagree quite strongly.

You are right, that some errors are harder to track down when  Inf & -Inf
work properly.
But in my oppinion the advantages really do abound.

for(s in 10^(0:9))cat("sd=",formatC(s,wid=5),"pnorm(0,1,sd)=",pnorm(0,1,s),"\n")
sd=     1 pnorm(0,1,sd)= 0.1586553 
sd=    10 pnorm(0,1,sd)= 0.4601722 
sd=   100 pnorm(0,1,sd)= 0.4960106 
sd=  1000 pnorm(0,1,sd)= 0.4996011 
sd= 1e+04 pnorm(0,1,sd)= 0.4999601 
sd= 1e+05 pnorm(0,1,sd)= 0.499996 
sd= 1e+06 pnorm(0,1,sd)= 0.4999996 
sd= 1e+07 pnorm(0,1,sd)= 0.5 
sd= 1e+08 pnorm(0,1,sd)= 0.5 
sd= 1e+09 pnorm(0,1,sd)= 0.5 


I think one of the merits of useful IEEE implementations is
continuity (in the mathematical sense)
and ``automatic fulfillment'' of limit theorems.

This has been very useful for me in several situations of S-plus code
development.  
And I have been longing for these to be available in R, as well.

    Ross> There are some other nits like this too.  Would it be a serious
    Ross> problem if I were to "fix" these?

yes (see above).


    Ross> Q2:

    Ross> At present, the functions returning information about discrete
    Ross> distributions work on integer (32-bit C int) arguments.  This
    Ross> limits the domain for these functions.  Has anyone found that
    Ross> that this is a problem?

not yet.
How much speed would be lost, if they worked for all double (and Inf, -Inf)
arguments?

- Martin
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


IEEE 754 Style Arithmetic
Martin Maechler