From p.dalgaard@biostat.ku.dk  Wed Dec 26 14:31:27 2001
From: p.dalgaard@biostat.ku.dk (Peter Dalgaard BSA)
Date: 26 Dec 2001 15:31:27 +0100
Subject: [Rd] bug with var(rep(1e30, 3)) (PR#1228)
In-Reply-To: <200112261207.NAA13242@pubhealth.ku.dk>
References: <200112261207.NAA13242@pubhealth.ku.dk>
Message-ID: <x2heqet5mo.fsf@blueberry.kubism.ku.dk>

paradis@isem.univ-montp2.fr writes:

> There seems to be a bug with var() when the argument is a vector with
> exactly three values of 1e30 (or close to this value). This does not happen
> with twice, four (or more) times this value, or another value.
> 
> 
> > var(rep(1e30, 3))
> [1] 2.971056e+28
> > var(rep(1.2e30, 3))
> [1] 2.971056e+28
> > var(rep(0.9e30, 3))
> [1] 2.971056e+28
> > var(rep(0.8e30, 3))
> [1] 0
> > var(rep(1e29, 3))
> [1] 0
> > var(rep(1e30, 2))
> [1] 0
> > var(rep(1e30, 4))
> [1] 0
> > var(rep(1e31, 3))
> [1] 0
> 
> 
> The bug is repeatable, and I got the same results with R 1.3.1 (binary from
> CRAN) and R 1.4.0 (binary from BD Ripley's site) both on Windows NT, and R
> 1.4.0 on Solaris (compiled from sources).

Not a bug, roundoff:

> (1e30+1e30+1e30)/3-1e30
[1] 1.407375e+14
> 3*((1e30+1e30+1e30)/3-1e30)^2/2
[1] 2.971056e+28

Relative errors of the order 1e-16 are completely expectable in
floating point arithmetic.

For illustration, look at it in octal, after removing a bunch of
trailing zeroes:

> x<-1e30/(8^16)
> structure(x, class="octmode")
[1] "144762623464043165"
> structure(x+x+x, class="octmode")
[1] "456730272634151540"

If you look carefully, you will see that a bit is lost in the process
since 5+5+5 is 17 octal.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[Rd] bug with var(rep(1e30, 3)) (PR#1228)
Peter Dalgaard BSA