Monday, March 17, 2008

R special cases x^2 and x^0.5

Pulled from src/main/arithmetic.c:
static SEXP real_unary(ARITHOP_TYPE code, SEXP s1, SEXP lcall)
{
int i, n;
SEXP ans;

switch (code) {
case PLUSOP: return s1;
case MINUSOP:
ans = duplicate(s1);
n = LENGTH(s1);
for (i = 0; i < n; i++)
REAL(ans)[i] = -REAL(s1)[i];
return ans;
default:
errorcall(lcall, _("invalid unary operator"));
}
return s1; /* never used; to keep -Wall happy */
}
Would this be noticeably faster if we redefined the case for MINUSOP as
case MINUSOP:
{
const double* y1 = REAL(s1)[i];
double* yAns;
ans = duplicate(s1);
yAns = REAL(ans);
n = LENGTH(s1);
for (i = 0; i < n; ++i)
yAns[i] = -y1[i];
return ans;
}
?
We have the following defines:
Rinternals.h: #define REAL(x) ((double *) DATAPTR(x))
Rinternals.h: #define DATAPTR(x) (((SEXPREC_ALIGN *) (x)) + 1)
So, it seems like this could theoretically save 2n additions per negation. Woo! Is gcc smart enough to notice the legality of this move and take advantage of it?

Poking around in arithmetic.c also shows that x^2.0 and x^0.5 are special-cased to use more efficient alternatives than pow while there's another special routine R_pow_di to handle raising a double to an integer power. Which is used by pi^3? Setting breakpoints on R_pow and R_pow_di shows R_pow is used. In fact, pi^as.integer(3) also uses R_pow.

I wonder what the distribution of powers used in production R code looks like. 2 is almost certainly the most popular. Would it be worthwhile special casing other integer powers or even including a generic check for an integer power? Do most implementations of pow do this already?

No comments: