a nasty error --> patch to seq.R [nr. 2]

[Martin Maechler]

>>>>> "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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._