a nasty error --> patch to seq.R [nr. 2]
Martin Maechler
Martin Maechler <maechler@stat.math.ethz.ch>
Sat, 5 Sep 1998 19:22:49 +0200
>>>>> "Peter" == Peter Holzer <holzer@stat.math.ethz.ch> writes:
Peter> ... The basic function seq doesn't work properly:
>> seq(1,6,by=3)
Peter> [1] 1 4 7
Peter> Looking at the source code of seq.default I found that strange
Peter> fuzz thing "+ 0.4". Taking that away the seq works fine.
Peter> My question is: Why has this fuzz thing been added?
because, otherwise, in some cases the sequence will be one element too
short..
Peter> If it was for some rounding problems I would suggest to replace
Peter> "0.4" by something much smaller.
yes; and depending on other things (see patch at end).
Actually, I had to sit down write up where which rounding errors can occur
in
n <- as.integer( (to - from) / by )
E.g., major cancellation happens when to & from are almost the same.
Additionally, I wrote a private function
seq.check <- function(from,to,by, print.diag = TRUE)
{
## Purpose: check seq(from, to, by = ...)
## -----------------------------------------------------
## Arguments: from < to; by > 0 [see ?seq ]
## -----------------------------------------------------
## Author: Martin Mächler, Date: 5 Sep 98, 12:08
if(from >= to) stop("seq.check(.) requires from < to")
if(by <= 0) stop("seq.check(.) requires by > 0")
n <- length(ss <- seq(from,to, by=by))
m <- ss[n]
res <- m <= to && to < m+by ##-- iff TRUE: all is fine.
if(print.diag) {
if(res)
cat(".")
else # m > to "last >" or m+by <= to "missed last"
cat("\nfrom=", formatC(from, dig=6,wid=-10),
", to-from=",formatC(to-from,dig=6,wid=-10),
" --> n=", formatC(n, wid=-4),": ",
if(m+by <= to)"missed last" else "last > to",
"; r.err=",
formatC(abs((if(m+by <= from) m+by else m)-to)/to,
wid= 8, dig=4),
sep="")
}
res
}
which revealed that even the version below can `very rarely' end up
one value too short.
However, in that case, you have
from + n * by > to ==> TRUE,
but
from + (n-1)*by + by > to ==> FALSE
[an example of the famous fact the the associative law of addition is
*NOT* guaranteed in computer arithmetic!]
which *is* rare.
One could argue that in that case, it would be preferable that the sequence
had one value more, one very slightly > to.
But I think, we rather want to guarantee
seq(from, to, by) <= to if by > 0
seq(from, to, by) >= to if by < 0
--------
Here, the patch promised.
The new seq.default(.) will be part of tomorrow's (or Monday's)
R-devel.tar.gz snapshot.
[Actually a snapshot that seems pretty stable..]
--- R-0.62.3/src/library/base/R/seq.R Tue Jun 16 14:32:08 1998
+++ R-0.63.0/src/library/base/R/seq.R Sat Sep 5 18:06:03 1998
@@ -16,9 +17,25 @@
if(missing(by))
from:to
else {
- n <- (to - from)/by + 0.4 # fuzz needed
- if(n < 0)
+ n <- (del <- to - from)/by
+ if(is.na(n)) {
+ if(by == 0 && del == 0)
+ return(from)
+ else stop("invalid (to - from)/by in seq(.)")
+ } else if(abs(n) > .Machine$integer.max)
+ stop("'by' argument is much too small")
+ else if(n < 0)
stop("Wrong sign in by= argument")
+
+ eps <- .Machine$double.eps *
+ max(1, max(abs(to),abs(from)) / abs(del))
+ n <- as.integer(n * (1 + eps))
+ if(eps*2*n >= 1)
+ warning(paste("seq.default(f,t,by): n=",n,
+ ": possibly imprecise intervals"))
+ if(by>0) while(from+ n*by > to) n <- n - 1
+ else while(from+ n*by < to) n <- n - 1
+
from + (0:n) * by
}
else if(length.out < 0)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._