Skip to content

Commit

Permalink
seq.default(from, to, *) and seq.int(*) when |from-to| is Inf
Browse files Browse the repository at this point in the history 
git-svn-id: https://svn.r-project.org/R/trunk@80483 00db46b3-68df-0310-9c12-caf00c1e9a41
  • Loading branch information
maechler committed Jun 11, 2021
1 parent a172045 commit 4086ae8
Show file tree
Hide file tree
Showing 4 changed files with 121 additions and 38 deletions.
4 changes: 4 additions & 0 deletions doc/NEWS.Rd
View file
Expand Up @@ -157,6 +157,10 @@

\item \code{.First} and \code{.Last} can again be set from the site
profile.

\item \code{seq.int(from, to, *)} and \code{seq.default(..)} now work
better in large range cases where \code{from-to} is infinite where
the two boundaries are finite.
}
}
}
Expand Down
29 changes: 22 additions & 7 deletions src/library/base/R/seq.R
View file
@@ -1,7 +1,7 @@
# File src/library/base/R/seq.R
# Part of the R package, https://www.R-project.org
#
# Copyright (C) 1995-2019 The R Core Team
# Copyright (C) 1995-2021 The R Core Team
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
Expand Down Expand Up @@ -66,7 +66,10 @@ seq.default <-
int <- FALSE
else if(!int)
storage.mode(by) <- "double"
n <- del/by # of length 1, as {from, to, by} are
n <- if(finite.del <- is.finite(del))
del/by # of length 1, as {from, to, by} are
else
to/by - from/by
if(!is.finite(n)) {
if(!is.na(by) && by == 0 && del == 0)
return(from)
Expand All @@ -77,15 +80,18 @@ seq.default <-
if(n > .Machine$integer.max)
stop("'by' argument is much too small")

dd <- abs(del)/max(abs(to), abs(from))
if (dd < 100*.Machine$double.eps) return(from)
if (finite.del && abs(del)/max(abs(to), abs(from)) < 100*.Machine$double.eps)
return(from) ## 100 is a fudge factor
if (int) {
n <- as.integer(n) # truncates
if (n >= 2L) cumsum(rep.int(c(from, by), c(1L, n))) else
from + (0L:n) * by
} else {
n <- as.integer(n + 1e-10)
x <- from + (0L:n) * by
x <- if(finite.del)
from + (0L:n) * by
else
(from/4 + (0L:n) * (by/4))*4
## correct for possible overshot because of fuzz
if(by > 0) pmin(x, to) else pmax(x, to)
}
Expand Down Expand Up @@ -119,8 +125,17 @@ seq.default <-
}
else {
if (intdel) storage.mode(from) <- "double"
by <- (to - from) / n1
as.vector(c(from, from + seq_len(length.out - 2L) * by, to))
del <- to - from
if(is.finite(del)) {
## by = del/n1
as.vector(c(from, from + seq_len(length.out - 2L) * (del/n1), to))
} else { # |del| = Inf, when from,to are ok (and large!)
from <- from/4 # and by = (to/n1 - from/n1) / 4
to <- to / 4
as.vector(c(from,
from + seq_len(length.out - 2L) * ((to-from)/n1),
to)) * 4
}
}
}
else as.vector(c(from, to))[seq_len(length.out)]
Expand Down
85 changes: 55 additions & 30 deletions src/main/seq.c
View file
@@ -1,6 +1,6 @@
/*
* R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1998-2020 The R Core Team.
* Copyright (C) 1998-2021 The R Core Team.
* Copyright (C) 1995-1998 Robert Gentleman and Ross Ihaka
*
* This program is free software; you can redistribute it and/or modify
Expand Down Expand Up @@ -352,7 +352,7 @@ SEXP attribute_hidden do_rep_int(SEXP call, SEXP op, SEXP args, SEXP rho)

if (DispatchOrEval(call, op, "rep", args, rho, &a, 0, 0))
return(a);

if (!isVector(ncopy))
error(_("invalid type (%s) for '%s' (must be a vector)"),
type2char(TYPEOF(ncopy)), "times");
Expand Down Expand Up @@ -428,7 +428,7 @@ SEXP attribute_hidden do_rep_len(SEXP call, SEXP op, SEXP args, SEXP rho)
}
UNPROTECT(1);
}

if (!isVector(s) && s != R_NilValue)
error(_("attempt to replicate non-vector"));

Expand Down Expand Up @@ -801,10 +801,8 @@ SEXP attribute_hidden do_rep(SEXP call, SEXP op, SEXP args, SEXP rho)
SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP ans = R_NilValue /* -Wall */, from, to, by, len, along;
int nargs = length(args), lf;
Rboolean One = (nargs == 1);
R_xlen_t i, lout = NA_INTEGER;
static SEXP do_seq_formals = NULL;
Rboolean One = (length(args) == 1); // *before* messing with args ..

/* DispatchOrEval internal generic: seq.int */
if (DispatchOrEval(call, op, "seq", args, rho, &ans, 0, 1))
Expand All @@ -814,6 +812,7 @@ SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)
We pretend this is
seq(from, to, by, length.out, along.with, ...)
*/
static SEXP do_seq_formals = NULL;
if (do_seq_formals == NULL)
do_seq_formals = allocFormalsList6(install("from"), install("to"),
install("by"), install("length.out"),
Expand All @@ -828,9 +827,8 @@ SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)
Rboolean
miss_from = (from == R_MissingArg),
miss_to = (to == R_MissingArg);

if(One && !miss_from) {
lf = length(from);
int lf = length(from);
if(lf == 1 && (TYPEOF(from) == INTSXP || TYPEOF(from) == REALSXP)) {
double rfrom = asReal(from);
if (!R_FINITE(rfrom))
Expand Down Expand Up @@ -860,7 +858,7 @@ SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)
}

if(lout == NA_INTEGER) {
double rfrom, rto, rby = asReal(by);
double rfrom, rto;
if(miss_from) rfrom = 1.0;
else {
if(length(from) != 1) errorcall(call, _("'%s' must be of length 1"), "from");
Expand All @@ -877,24 +875,29 @@ SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)
}
if(by == R_MissingArg)
ans = seq_colon(rfrom, rto, call);
else {
else { // 'by' specified
if(length(by) != 1) errorcall(call, _("'%s' must be of length 1"), "by");
double del = rto - rfrom;
if(del == 0.0 && rto == 0.0) {
ans = to; // is *not* missing in this case
goto done;
}
/* printf("from = %f, to = %f, by = %f\n", rfrom, rto, rby); */
double n = del/rby;
double n, rby = asReal(by);
Rboolean finite_del = R_FINITE(del);
if(finite_del) {
n = del/rby;
} else { // overflow in (to - from) when both are finite
n = rto/rby - rfrom/rby;
}
if(!R_FINITE(n)) {
if(del == 0.0 && rby == 0.0) {
ans = miss_from ? ScalarReal(rfrom) : from;
goto done;
} else
errorcall(call, _("invalid '(to - from)/by'"));
}
double dd = fabs(del)/fmax2(fabs(rto), fabs(rfrom));
if(dd < 100 * DBL_EPSILON) {
// inherited from seq.default() but "fudgy"
if(finite_del && fabs(del)/fmax2(fabs(rto), fabs(rfrom)) < 100 * DBL_EPSILON) {
ans = miss_from ? ScalarReal(rfrom) : from;
goto done;
}
Expand All @@ -909,8 +912,8 @@ SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)

R_xlen_t nn;
if((miss_from || TYPEOF(from) == INTSXP) &&
(miss_to || TYPEOF(to) == INTSXP) &&
TYPEOF(by) == INTSXP) {
(miss_to || TYPEOF(to) == INTSXP) && TYPEOF(by) == INTSXP)
{
int *ia, ifrom = miss_from ? (int)rfrom : asInteger(from),
iby = asInteger(by);
/* With the current limits on integers and FEPS
Expand All @@ -928,8 +931,15 @@ SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)
nn = (int)(n + FEPS);
ans = allocVector(REALSXP, nn+1);
double *ra = REAL(ans);
for(i = 0; i <= nn; i++)
ra[i] = rfrom + (double)i * rby;
if(finite_del)
for(i = 0; i <= nn; i++)
ra[i] = rfrom + (double)i * rby;
else { // |from - to| is infinite, but n = (from-to)/by is not
rfrom /= 4.;
rby /= 4.;
for(i = 0; i <= nn; i++) // ldexp(Y, 2) := Y * 2^2 = 4 Y
ra[i] = ldexp(rfrom + (double)i * rby, 2);
}
/* Added in 2.9.0 */
if (nn > 0)
if((rby > 0 && ra[nn] > rto) || (rby < 0 && ra[nn] < rto))
Expand All @@ -940,38 +950,53 @@ SEXP attribute_hidden do_seq(SEXP call, SEXP op, SEXP args, SEXP rho)
ans = allocVector(INTSXP, 0);
} else if (One) {
ans = seq_colon(1.0, (double)lout, call);
} else if (by == R_MissingArg) { // and len := length.out specified
} else if (by == R_MissingArg) { // and len := length.out specified, >= 1
double rfrom = asReal(from), rto = asReal(to), rby = 0; // -Wall
if(miss_to) rto = rfrom + (double)lout - 1;
if(miss_from) rfrom = rto - (double)lout + 1;
if(!R_FINITE(rfrom)) errorcall(call, _("'%s' must be a finite number"), "from");
if(!R_FINITE(rto)) errorcall(call, _("'%s' must be a finite number"), "to");
if(lout > 2) rby = (rto - rfrom)/(double)(lout - 1);
Rboolean finite_del = 0;
if(lout > 2) { // only then, use 'by'
double nint = (double)(lout - 1);
if((finite_del = R_FINITE(rby = (rto - rfrom))))
rby /= nint;
else // overflow in (to - from), nint >= 2 => finite 'by'
rby = (rto/nint - rfrom/nint);
}
if(rfrom <= INT_MAX && rfrom >= INT_MIN &&
rto <= INT_MAX && rto >= INT_MIN &&
rfrom == (int)rfrom &&
(lout <= 1 || rto == (int)rto) &&
(lout <= 2 || rby == (int)rby)) {
ans = allocVector(INTSXP, lout);
if(lout > 0) INTEGER(ans)[0] = (int)rfrom;
INTEGER(ans)[0] = (int)rfrom;
if(lout > 1) INTEGER(ans)[lout - 1] = (int)rto;
if(lout > 2)
for(i = 1; i < lout-1; i++) {
// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();
INTEGER(ans)[i] = (int)(rfrom + (double)i*rby);
}
} else {
ans = allocVector(REALSXP, lout);
if(lout > 0) REAL(ans)[0] = rfrom;
if(lout > 1) REAL(ans)[lout - 1] = rto;
if(lout > 2) {
rby = (rto - rfrom)/(double)(lout - 1);
for(i = 1; i < lout-1; i++) {
// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();
REAL(ans)[i] = rfrom + (double)i*rby;
ans = allocVector(REALSXP, lout);
REAL(ans)[0] = rfrom;
if(lout > 1) REAL(ans)[lout - 1] = rto;
if(lout > 2) {
if(finite_del)
for(i = 1; i < lout-1; i++) {
// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();
REAL(ans)[i] = rfrom + (double)i*rby;
}
else { // del:=(from - to) is infinite, but n = del/by is not
rfrom /= 4.; // or ldexp(*, -2) for speed
rby /= 4.;
for(i = 1; i < lout-1; i++) {
// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();
REAL(ans)[i] = ldexp(rfrom + (double)i*rby, 2); // ldexp(y,2) = y * 4
}
}
}
}
}
} else if (miss_to) {
double rfrom = asReal(from), rby = asReal(by), rto;
if(miss_from) rfrom = 1.0;
Expand Down
41 changes: 40 additions & 1 deletion tests/reg-tests-1d.R
View file
Expand Up @@ -5071,7 +5071,46 @@ if (l10n_info()$"Latin-1" && localeToCharset()=="ISO8859-1") {
# l10n_info() would report Latin-1 when that is the code page
y <- "\xfc"
stopifnot(y == encodeString(y))
}
}


## seq(from, to, *) i.e. seq.default() *and* seq.int(..) in case of large
## from & to, notably *infinite (to - from) :
seq (-1.5e308, 1e308, by=1e307) # gave error in R <= 4.1.0
seq (-1.5e308, 1.6e308, length.out=33)# all Inf apart from first & last
## and these two where identical to seq(), aka seq.default()
seq.int(-1.5e308, 1e308, by=1e307)
seq.int(-1.5e308, 1.6e308, length.out=33)
## more systematically: ------------------------------------------------------
## test a series
B <- .Machine$double.xmax; B. <- 1.79769e308
Lby <- lapply(1:25, function(N) seq.int(-.99*B, B , by = N*2e306) / B)
LbyR <- lapply(1:25, function(N) seq.default(-.99*B, B , by = N*2e306) / B)
Llen <- lapply(2:26, function(N) seq.int(- B, B., length.out = N) / B)
LleR <- lapply(2:26, function(N) seq.default(- B, B., length.out = N) / B)
## first diff should be constant
relEdiff <- function(L)
vapply(lapply(L, diff), function(x) {m <- mean(x); max(abs(x - m) / m) }, 1.23)
by <- 1e307
stopifnot(exprs = {
## C = R : seq.int() <==> seq.default :
all.equal(Lby , LbyR, tol=1e-15)
all.equal(Llen, LleR, tol=1e-15)
## by :
abs(diff(s <- seq.int(-1.5e308, 1e308, by=by))/by - 1) < 1e-14
is.matrix(rng <- vapply(Lby, range, numeric(2)))
all.equal(rep(-.99, 25), rng[1,])
0.79 <= rng[2,] ; rng[2,] <= 0.991
## length.out
is.matrix(rng <- vapply(Llen, range, numeric(2)))
all.equal(rep( -1, 25), rng[1,])
all.equal(rep(B./B, 25), rng[2,])
## first diff s:
abs(relEdiff(Llen) / lengths(Llen)^1.5) < 2e-16 # see max ~3e-17
abs(relEdiff(Llen) / lengths(Llen)^1.5) < 2e-16 # (ditto)
})



## keep at end
rbind(last = proc.time() - .pt,
Expand Down

0 comments on commit 4086ae8

Please sign in to comment.