/*
* R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
* Copyright (C) 1998-2006 The R Development Core Team.
* Copyright (C) 2004 The R Foundation
*
* 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
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street Fifth Floor, Boston, MA 02110-1301 USA
*/
/* <UTF8> char here is handled as a whole string.
Does rely on strcoll being correct.
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <Defn.h> /* => Utils.h with the protos from here */
#include <Rmath.h>
/*--- Part I: Comparison Utilities ---*/
static int icmp(int x, int y, Rboolean nalast)
{
if (x == NA_INTEGER && y == NA_INTEGER) return 0;
if (x == NA_INTEGER)return nalast?1:-1;
if (y == NA_INTEGER)return nalast?-1:1;
if (x < y) return -1;
if (x > y) return 1;
return 0;
}
static int rcmp(double x, double y, Rboolean nalast)
{
int nax = ISNAN(x), nay = ISNAN(y);
if (nax && nay) return 0;
if (nax) return nalast?1:-1;
if (nay) return nalast?-1:1;
if (x < y) return -1;
if (x > y) return 1;
return 0;
}
static int ccmp(Rcomplex x, Rcomplex y, Rboolean nalast)
{
int nax = ISNAN(x.r), nay = ISNAN(y.r);
/* compare real parts */
if (nax && nay) return 0;
if (nax) return nalast?1:-1;
if (nay) return nalast?-1:1;
if (x.r < y.r) return -1;
if (x.r > y.r) return 1;
/* compare complex parts */
nax = ISNAN(x.i); nay = ISNAN(y.i);
if (nax && nay) return 0;
if (nax) return nalast?1:-1;
if (nay) return nalast?-1:1;
if (x.i < y.i) return -1;
if (x.i > y.i) return 1;
return 0; /* equal */
}
static int scmp(SEXP x, SEXP y, Rboolean nalast)
{
if (x == NA_STRING && y == NA_STRING) return 0;
if (x == NA_STRING) return nalast?1:-1;
if (y == NA_STRING) return nalast?-1:1;
return STRCOLL(CHAR(x), CHAR(y));
}
Rboolean isUnsorted(SEXP x)
{
int n, i;
if (!isVectorAtomic(x))
error(_("only atomic vectors can be tested to be sorted"));
n = LENGTH(x);
if(n >= 2)
switch (TYPEOF(x)) {
/* NOTE: x must have no NAs {is.na(.) in R};
hence be faster than `rcmp()', `icmp()' for these two cases */
case LGLSXP:
case INTSXP:
for(i = 0; i+1 < n ; i++)
if(INTEGER(x)[i] > INTEGER(x)[i+1])
return TRUE;
break;
case REALSXP:
for(i = 0; i+1 < n ; i++)
if(REAL(x)[i] > REAL(x)[i+1])
return TRUE;
break;
case CPLXSXP:
for(i = 0; i+1 < n ; i++)
if(ccmp(COMPLEX(x)[i], COMPLEX(x)[i+1], TRUE) > 0)
return TRUE;
break;
case STRSXP:
for(i = 0; i+1 < n ; i++)
if(scmp(STRING_ELT(x, i ),
STRING_ELT(x,i+1), TRUE) > 0)
return TRUE;
break;
default:
UNIMPLEMENTED_TYPE("isUnsorted", x);
}
return FALSE;/* sorted */
}
SEXP attribute_hidden do_isunsorted(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP ans;
checkArity(op, args);
ans = allocVector(LGLSXP, 1);
LOGICAL(ans)[0] = isUnsorted(CAR(args));
return ans;
}
/*--- Part II: Complete (non-partial) Sorting ---*/
/* SHELLsort -- corrected from R. Sedgewick `Algorithms in C'
* (version of BDR's lqs():*/
#define sort_body \
Rboolean nalast=TRUE; \
int i, j, h; \
\
for (h = 1; h <= n / 9; h = 3 * h + 1); \
for (; h > 0; h /= 3) \
for (i = h; i < n; i++) { \
v = x[i]; \
j = i; \
while (j >= h && TYPE_CMP(x[j - h], v, nalast) > 0) \
{ x[j] = x[j - h]; j -= h; } \
x[j] = v; \
}
void R_isort(int *x, int n)
{
int v;
#define TYPE_CMP icmp
sort_body
#undef TYPE_CMP
}
void R_rsort(double *x, int n)
{
double v;
#define TYPE_CMP rcmp
sort_body
#undef TYPE_CMP
}
void R_csort(Rcomplex *x, int n)
{
Rcomplex v;
#define TYPE_CMP ccmp
sort_body
#undef TYPE_CMP
}
/* used in platform.c */
void attribute_hidden ssort(SEXP *x, int n)
{
SEXP v;
#define TYPE_CMP scmp
sort_body
#undef TYPE_CMP
}
void rsort_with_index(double *x, int *indx, int n)
{
double v;
int i, j, h, iv;
for (h = 1; h <= n / 9; h = 3 * h + 1);
for (; h > 0; h /= 3)
for (i = h; i < n; i++) {
v = x[i]; iv = indx[i];
j = i;
while (j >= h && rcmp(x[j - h], v, TRUE) > 0)
{ x[j] = x[j - h]; indx[j] = indx[j-h]; j -= h; }
x[j] = v; indx[j] = iv;
}
}
void revsort(double *a, int *ib, int n)
{
/* Sort a[] into descending order by "heapsort";
* sort ib[] alongside;
* if initially, ib[] = 1...n, it will contain the permutation finally
*/
int l, j, ir, i;
double ra;
int ii;
if (n <= 1) return;
a--; ib--;
l = (n >> 1) + 1;
ir = n;
for (;;) {
if (l > 1) {
l = l - 1;
ra = a[l];
ii = ib[l];
}
else {
ra = a[ir];
ii = ib[ir];
a[ir] = a[1];
ib[ir] = ib[1];
if (--ir == 1) {
a[1] = ra;
ib[1] = ii;
return;
}
}
i = l;
j = l << 1;
while (j <= ir) {
if (j < ir && a[j] > a[j + 1]) ++j;
if (ra > a[j]) {
a[i] = a[j];
ib[i] = ib[j];
j += (i = j);
}
else
j = ir + 1;
}
a[i] = ra;
ib[i] = ii;
}
}
SEXP attribute_hidden do_sort(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP ans;
Rboolean decreasing;
checkArity(op, args);
decreasing = asLogical(CADR(args));
if(decreasing == NA_LOGICAL)
error(_("'decreasing' must be TRUE or FALSE"));
if(CAR(args) == R_NilValue) return R_NilValue;
if(!isVectorAtomic(CAR(args)))
errorcall(call, _("only atomic vectors can be sorted"));
if(TYPEOF(CAR(args)) == RAWSXP)
errorcall(call, _("raw vectors cannot be sorted"));
/* we need consistent behaviour here, including dropping attibutes,
so as from 2.3.0 we always duplicate. */
ans = duplicate(CAR(args));
SET_ATTRIB(ans, R_NilValue); /* this is never called with names */
sortVector(ans, decreasing);
return(ans);
}
/* faster versions of shellsort, following Sedgewick (1986) */
static const int incs[16] = {1073790977, 268460033, 67121153, 16783361, 4197377,
1050113, 262913, 65921, 16577, 4193, 1073, 281, 77,
23, 8, 1};
#define sort2_body \
for (h = incs[t]; t < 16; h = incs[++t]) \
for (i = h; i < n; i++) { \
v = x[i]; \
j = i; \
while (j >= h && x[j - h] less v) { x[j] = x[j - h]; j -= h; } \
x[j] = v; \
}
static void R_isort2(int *x, int n, Rboolean decreasing)
{
int v;
int i, j, h, t;
for (t = 0; incs[t] > n; t++);
if(decreasing)
#define less <
sort2_body
#undef less
else
#define less >
sort2_body
#undef less
}
static void R_rsort2(double *x, int n, Rboolean decreasing)
{
double v;
int i, j, h, t;
for (t = 0; incs[t] > n; t++);
if(decreasing)
#define less <
sort2_body
#undef less
else
#define less >
sort2_body
#undef less
}
static void R_csort2(Rcomplex *x, int n, Rboolean decreasing)
{
Rcomplex v;
int i, j, h, t;
for (t = 0; incs[t] > n; t++);
for (h = incs[t]; t < 16; h = incs[++t])
for (i = h; i < n; i++) {
v = x[i];
j = i;
if(decreasing)
while (j >= h && (x[j - h].r < v.r ||
(x[j - h].r == v.r && x[j - h].i < v.i)))
{ x[j] = x[j - h]; j -= h; }
else
while (j >= h && (x[j - h].r > v.r ||
(x[j - h].r == v.r && x[j - h].i > v.i)))
{ x[j] = x[j - h]; j -= h; }
x[j] = v;
}
}
static void ssort2(SEXP *x, int n, Rboolean decreasing)
{
SEXP v;
int i, j, h, t;
for (t = 0; incs[t] > n; t++);
for (h = incs[t]; t < 16; h = incs[++t])
for (i = h; i < n; i++) {
v = x[i];
j = i;
if(decreasing)
while (j >= h && scmp(x[j - h], v, TRUE) < 0)
{ x[j] = x[j - h]; j -= h; }
else
while (j >= h && scmp(x[j - h], v, TRUE) > 0)
{ x[j] = x[j - h]; j -= h; }
x[j] = v;
}
}
void sortVector(SEXP s, Rboolean decreasing)
{
int n = LENGTH(s);
if (n >= 2 && (decreasing || isUnsorted(s)))
switch (TYPEOF(s)) {
case LGLSXP:
case INTSXP:
R_isort2(INTEGER(s), n, decreasing);
break;
case REALSXP:
R_rsort2(REAL(s), n, decreasing);
break;
case CPLXSXP:
R_csort2(COMPLEX(s), n, decreasing);
break;
case STRSXP:
ssort2(STRING_PTR(s), n, decreasing);
break;
default:
UNIMPLEMENTED_TYPE("sortVector", s);
}
}
/*--- Part III: Partial Sorting ---*/
/*
Partial sort so that x[k] is in the correct place, smaller to left,
larger to right
NOTA BENE: k < n required, and *not* checked here but in do_psort();
----- infinite loop possible otherwise!
*/
#define psort_body \
Rboolean nalast=TRUE; \
int L, R, i, j; \
\
for (L = lo, R = hi; L < R; ) { \
v = x[k]; \
for(i = L, j = R; i <= j;) { \
while (TYPE_CMP(x[i], v, nalast) < 0) i++; \
while (TYPE_CMP(v, x[j], nalast) < 0) j--; \
if (i <= j) { w = x[i]; x[i++] = x[j]; x[j--] = w; }\
} \
if (j < k) L = i; \
if (k < i) R = j; \
}
static void iPsort2(int *x, int lo, int hi, int k)
{
int v, w;
#define TYPE_CMP icmp
psort_body
#undef TYPE_CMP
}
static void rPsort2(double *x, int lo, int hi, int k)
{
double v, w;
#define TYPE_CMP rcmp
psort_body
#undef TYPE_CMP
}
static void cPsort2(Rcomplex *x, int lo, int hi, int k)
{
Rcomplex v, w;
#define TYPE_CMP ccmp
psort_body
#undef TYPE_CMP
}
static void sPsort2(SEXP *x, int lo, int hi, int k)
{
SEXP v, w;
#define TYPE_CMP scmp
psort_body
#undef TYPE_CMP
}
/* elements of ind are 1-based, lo and hi are 0-based */
static void Psort(SEXP x, int lo, int hi, int k)
{
/* Rprintf("looking for index %d in (%d, %d)\n", k, lo, hi);*/
switch (TYPEOF(x)) {
case LGLSXP:
case INTSXP:
iPsort2(INTEGER(x), lo, hi, k);
break;
case REALSXP:
rPsort2(REAL(x), lo, hi, k);
break;
case CPLXSXP:
cPsort2(COMPLEX(x), lo, hi, k);
break;
case STRSXP:
sPsort2(STRING_PTR(x), lo, hi, k);
break;
default:
UNIMPLEMENTED_TYPE("Psort", x);
}
}
static void Psort0(SEXP x, int lo, int hi, int *ind, int k)
{
if(k < 1 || hi-lo < 1) return;
if(k <= 1)
Psort(x, lo, hi, ind[0]-1);
else {
/* Look for index nearest the centre of the range */
int i, this = 0, mid = (lo+hi)/2, z;
for(i = 0; i < k; i++)
if(ind[i]-1 <= mid) this = i;
z = ind[this]-1;
Psort(x, lo, hi, z);
Psort0(x, lo, z-1, ind, this);
Psort0(x, z+1, hi, ind+this+1, k-this-1);
}
}
/* Needed for mistaken decision to put these in the API */
void iPsort(int *x, int n, int k)
{
iPsort2(x, 0, n-1, k);
}
void rPsort(double *x, int n, int k)
{
rPsort2(x, 0, n-1, k);
}
void cPsort(Rcomplex *x, int n, int k)
{
cPsort2(x, 0, n-1, k);
}
/* FUNCTION psort(x, indices) */
SEXP attribute_hidden do_psort(SEXP call, SEXP op, SEXP args, SEXP rho)
{
int i, k, n;
int *l;
checkArity(op, args);
if (!isVectorAtomic(CAR(args)))
errorcall(call, _("only atomic vectors can be sorted"));
if(TYPEOF(CAR(args)) == RAWSXP)
errorcall(call, _("raw vectors cannot be sorted"));
n = LENGTH(CAR(args));
SETCADR(args, coerceVector(CADR(args), INTSXP));
l = INTEGER(CADR(args));
k = LENGTH(CADR(args));
for (i = 0; i < k; i++) {
if (l[i] == NA_INTEGER)
errorcall(call, _("NA index"));
if (l[i] < 1 || l[i] > n)
errorcall(call, _("index %d outside bounds"), l[i]);
}
SETCAR(args, duplicate(CAR(args)));
SET_ATTRIB(CAR(args), R_NilValue); /* remove all attributes */
Psort0(CAR(args), 0, n - 1, l, k);
return CAR(args);
}
/*--- Part IV : Rank & Order ---*/
static int equal(int i, int j, SEXP x, Rboolean nalast)
{
int c=-1;
switch (TYPEOF(x)) {
case LGLSXP:
case INTSXP:
c = icmp(INTEGER(x)[i], INTEGER(x)[j], nalast);
break;
case REALSXP:
c = rcmp(REAL(x)[i], REAL(x)[j], nalast);
break;
case CPLXSXP:
c = ccmp(COMPLEX(x)[i], COMPLEX(x)[j], nalast);
break;
case STRSXP:
c = scmp(STRING_ELT(x, i), STRING_ELT(x, j), nalast);
break;
default:
UNIMPLEMENTED_TYPE("equal", x);
break;
}
if (c == 0)
return 1;
return 0;
}
static int greater(int i, int j, SEXP x, Rboolean nalast, Rboolean decreasing)
{
int c = -1;
switch (TYPEOF(x)) {
case LGLSXP:
case INTSXP:
c = icmp(INTEGER(x)[i], INTEGER(x)[j], nalast);
break;
case REALSXP:
c = rcmp(REAL(x)[i], REAL(x)[j], nalast);
break;
case CPLXSXP:
c = ccmp(COMPLEX(x)[i], COMPLEX(x)[j], nalast);
break;
case STRSXP:
c = scmp(STRING_ELT(x, i), STRING_ELT(x, j), nalast);
break;
default:
UNIMPLEMENTED_TYPE("greater", x);
break;
}
if (decreasing) c = -c;
if (c > 0 || (c == 0 && j < i)) return 1; else return 0;
}
/* listgreater(): used as greater_sub in orderVector() in do_order(...) */
static int listgreater(int i, int j, SEXP key, Rboolean nalast,
Rboolean decreasing)
{
SEXP x;
int c = -1;
while (key != R_NilValue) {
x = CAR(key);
switch (TYPEOF(x)) {
case LGLSXP:
case INTSXP:
c = icmp(INTEGER(x)[i], INTEGER(x)[j], nalast);
break;
case REALSXP:
c = rcmp(REAL(x)[i], REAL(x)[j], nalast);
break;
case CPLXSXP:
c = ccmp(COMPLEX(x)[i], COMPLEX(x)[j], nalast);
break;
case STRSXP:
c = scmp(STRING_ELT(x, i), STRING_ELT(x, j), nalast);
break;
default:
UNIMPLEMENTED_TYPE("listgreater", x);
}
if (decreasing) c = -c;
if (c > 0)
return 1;
if (c < 0)
return 0;
key = CDR(key);
}
if (c == 0 && i < j) return 0; else return 1;
}
/* Needs indx set to 1...n initially */
static void orderVector(int *indx, int n, SEXP key, Rboolean nalast,
Rboolean decreasing, int greater_sub())
{
int i, j, h, t;
int itmp;
for (t = 0; incs[t] > n; t++);
for (h = incs[t]; t < 16; h = incs[++t])
for (i = h; i < n; i++) {
itmp = indx[i];
j = i;
while (j >= h &&
greater_sub(indx[j - h], itmp, key, nalast^decreasing,
decreasing)) {
indx[j] = indx[j - h];
j -= h;
}
indx[j] = itmp;
}
}
#define sort2_with_index \
for (h = incs[t]; t < 16; h = incs[++t]) \
for (i = lo + h; i <= hi; i++) { \
itmp = indx[i]; \
j = i; \
while (j >= lo + h && less(indx[j - h], itmp)) { \
indx[j] = indx[j - h]; j -= h; } \
indx[j] = itmp; \
}
/* Needs indx set to 1...n initially.
Also used by do_options.
*/
void attribute_hidden
orderVector1(int *indx, int n, SEXP key, Rboolean nalast, Rboolean decreasing)
{
int c, i, j, h, t, lo = 0, hi = n-1;
int itmp, *isna, numna = 0;
int *ix = NULL /* -Wall */;
double *x = NULL /* -Wall */;
Rcomplex *cx = NULL /* -Wall */;
SEXP *sx = NULL /* -Wall */;
switch (TYPEOF(key)) {
case LGLSXP:
case INTSXP:
ix = INTEGER(key);
break;
case REALSXP:
x = REAL(key);
break;
case STRSXP:
sx = STRING_PTR(key);
break;
case CPLXSXP:
cx = COMPLEX(key);
break;
}
/* First sort NAs to one end */
isna = (int *) malloc(n * sizeof(int));
switch (TYPEOF(key)) {
case LGLSXP:
case INTSXP:
for (i = 0; i < n; i++) isna[i] = (ix[i] == NA_INTEGER);
break;
case REALSXP:
for (i = 0; i < n; i++) isna[i] = ISNAN(x[i]);
break;
case STRSXP:
for (i = 0; i < n; i++) isna[i] = (sx[i] == NA_STRING);
break;
case CPLXSXP:
for (i = 0; i < n; i++) isna[i] = ISNAN(cx[i].r) || ISNAN(cx[i].i);
break;
default:
UNIMPLEMENTED_TYPE("orderVector1", key);
}
for (i = 0; i < n; i++) numna += isna[i];
if(numna)
switch (TYPEOF(key)) {
case LGLSXP:
case INTSXP:
case REALSXP:
case STRSXP:
case CPLXSXP:
if (!nalast) for (i = 0; i < n; i++) isna[i] = !isna[i];
for (t = 0; incs[t] > n; t++);
#define less(a, b) (isna[a] > isna[b] || (isna[a] == isna[b] && a > b))
sort2_with_index
#undef less
if(nalast) hi -= numna; else lo += numna;
}
/* Shell sort isn't stable, so add test on index */
for (t = 0; incs[t] > hi-lo+1; t++);
switch (TYPEOF(key)) {
case LGLSXP:
case INTSXP:
if (decreasing) {
#define less(a, b) (ix[a] < ix[b] || (ix[a] == ix[b] && a > b))
sort2_with_index
#undef less
} else {
#define less(a, b) (ix[a] > ix[b] || (ix[a] == ix[b] && a > b))
sort2_with_index
#undef less
}
break;
case REALSXP:
if (decreasing) {
#define less(a, b) (x[a] < x[b] || (x[a] == x[b] && a > b))
sort2_with_index
#undef less
} else {
#define less(a, b) (x[a] > x[b] || (x[a] == x[b] && a > b))
sort2_with_index
#undef less
}
break;
case CPLXSXP:
if (decreasing) {
#define less(a, b) (ccmp(cx[a], cx[b], 0) < 0 || (cx[a].r == cx[b].r && cx[a].i == cx[b].i && a > b))
sort2_with_index
#undef less
} else {
#define less(a, b) (ccmp(cx[a], cx[b], 0) > 0 || (cx[a].r == cx[b].r && cx[a].i == cx[b].i && a > b))
sort2_with_index
#undef less
}
break;
case STRSXP:
if (decreasing)
#define less(a, b) (c=STRCOLL(CHAR(sx[a]),CHAR(sx[b])), c < 0 || (c == 0 && a > b))
sort2_with_index
#undef less
else
#define less(a, b) (c=STRCOLL(CHAR(sx[a]),CHAR(sx[b])), c > 0 || (c == 0 && a > b))
sort2_with_index
#undef less
break;
default:
#define less(a, b) greater(a, b, key, nalast^decreasing, decreasing)
sort2_with_index
#undef less
}
free(isna);
}
/* FUNCTION order(...) */
SEXP attribute_hidden do_order(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP ap, ans;
int i, n = -1, narg = 0;
Rboolean nalast, decreasing;
nalast = asLogical(CAR(args));
if(nalast == NA_LOGICAL)
error(_("invalid '%s' value"), "na.last");
args = CDR(args);
decreasing = asLogical(CAR(args));
if(decreasing == NA_LOGICAL)
error(_("'decreasing' must be TRUE or FALSE"));
args = CDR(args);
if (args == R_NilValue)
return R_NilValue;
if (isVector(CAR(args)))
n = LENGTH(CAR(args));
for (ap = args; ap != R_NilValue; ap = CDR(ap), narg++) {
if (!isVector(CAR(ap)))
errorcall(call, _("argument %d is not a vector"), narg + 1);
if (LENGTH(CAR(ap)) != n)
errorcall(call, _("argument lengths differ"));
}
ans = allocVector(INTSXP, n);
if (n != 0) {
for (i = 0; i < n; i++) INTEGER(ans)[i] = i;
if(narg == 1)
orderVector1(INTEGER(ans), n, CAR(args), nalast, decreasing);
else
orderVector(INTEGER(ans), n, args, nalast, decreasing, listgreater);
for (i = 0; i < n; i++) INTEGER(ans)[i]++;
}
return ans;
}
/* FUNCTION: rank(x) */
SEXP attribute_hidden do_rank(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP rank, indx, x;
int *in;
double *rk;
int i, j, k, n;
char *ties_str;
enum {AVERAGE, MAX, MIN} ties_kind = AVERAGE;
checkArity(op, args);
if (args == R_NilValue)
return R_NilValue;
x = CAR(args);
if (!isVectorAtomic(x))
errorcall(call, _("argument is not an atomic vector"));
if(TYPEOF(x) == RAWSXP)
errorcall(call, _("raw vectors cannot be sorted"));
n = LENGTH(x);
PROTECT(indx = allocVector(INTSXP, n));
PROTECT(rank = allocVector(REALSXP, n));
UNPROTECT(2);
ties_str = CHAR(STRING_ELT(coerceVector(CADR(args), STRSXP), 0));
if(!strcmp(ties_str, "average")) ties_kind = AVERAGE;
else if(!strcmp(ties_str, "max")) ties_kind = MAX;
else if(!strcmp(ties_str, "min")) ties_kind = MIN;
else error(_("invalid ties.method for rank() [should never happen]"));
if (n > 0) {
in = INTEGER(indx);
rk = REAL(rank);
for (i = 0; i < n; i++)
in[i] = i;
orderVector1(in, n, x, TRUE, FALSE);
i = 0;
while (i < n) {
j = i;
while ((j < n - 1) && equal(in[j], in[j + 1], x, TRUE))
j++;
if (i != j) { /* ties */
switch(ties_kind) {
case AVERAGE:
for (k = i; k <= j; k++)
rk[in[k]] = (i + j + 2) / 2.; break;
case MAX:
for (k = i; k <= j; k++) rk[in[k]] = j+1; break;
case MIN:
for (k = i; k <= j; k++) rk[in[k]] = i+1; break;
}
}
else
rk[in[i]] = i + 1;
i = j + 1;
}
}
return rank;
}
#include <R_ext/RS.h>
SEXP attribute_hidden do_radixsort(SEXP call, SEXP op, SEXP args, SEXP rho)
{
SEXP x, ans;
Rboolean nalast, decreasing;
unsigned int *cnts;
int i, n, tmp, xmax = NA_INTEGER, xmin = NA_INTEGER, off, napos;
checkArity(op, args);
x = CAR(args);
nalast = asLogical(CADR(args));
if(nalast == NA_LOGICAL)
error(_("invalid '%s' value"), "na.last");
decreasing = asLogical(CADDR(args));
if(decreasing == NA_LOGICAL)
error(_("'decreasing' must be TRUE or FALSE"));
off = nalast^decreasing ? 0 : 1;
n = LENGTH(x);
PROTECT(ans = allocVector(INTSXP, n));
for(i = 0; i < n; i++) {
tmp = INTEGER(x)[i];
if(tmp == NA_INTEGER) continue;
if(tmp < 0) errorcall(call, _("negative value in 'x'"));
if(xmax == NA_INTEGER || tmp > xmax) xmax = tmp;
if(xmin == NA_INTEGER || tmp < xmin) xmin = tmp;
}
if(xmin == NA_INTEGER) { /* all NAs, so nothing to do */
for(i = 0; i < n; i++) INTEGER(ans)[i] = i+1;
UNPROTECT(1);
return ans;
}
xmax -= xmin;
if(xmax > 100000) errorcall(call, _("too large a range of values in 'x'"));
napos = off ? 0 : xmax + 1;
off -= xmin;
/* alloca is fine here: we know this is small */
cnts = (unsigned int *) alloca((xmax+1)*sizeof(unsigned int));
R_CheckStack();
for(i = 0; i <= xmax+1; i++) cnts[i] = 0;
for(i = 0; i < n; i++) {
if(INTEGER(x)[i] == NA_INTEGER) cnts[napos]++;
else cnts[off+INTEGER(x)[i]]++;
}
for(i = 1; i <= xmax+1; i++) cnts[i] += cnts[i-1];
if(decreasing)
for(i = 0; i < n; i++){
tmp = INTEGER(x)[i];
INTEGER(ans)[n-(cnts[(tmp==NA_INTEGER) ? napos : off+tmp]--)] = i+1;
}
else
for(i = n-1; i >= 0; i--) {
tmp = INTEGER(x)[i];
INTEGER(ans)[--cnts[(tmp==NA_INTEGER) ? napos : off+tmp]] = i+1;
}
UNPROTECT(1);
return ans;
}
syntax highlighted by Code2HTML, v. 0.9.1