| /* |
| * R : A Computer Language for Statistical Data Analysis |
| * Copyright (C) 1998--2019 The R Core Team. |
| * Copyright (C) 2003--2019 The R Foundation |
| * Copyright (C) 1995--1997 Robert Gentleman and Ross Ihaka |
| * |
| * 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, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| */ |
| |
| #ifdef HAVE_CONFIG_H |
| #include <config.h> |
| #endif |
| |
| // LDBL_EPSILON |
| #include <float.h> |
| |
| /* interval at which to check interrupts, a guess */ |
| #define NINTERRUPT 10000000 |
| |
| |
| #ifdef __OpenBSD__ |
| /* for definition of "struct exception" in math.h */ |
| # define __LIBM_PRIVATE |
| #endif |
| #include <Defn.h> /*-> Arith.h -> math.h */ |
| #ifdef __OpenBSD__ |
| # undef __LIBM_PRIVATE |
| #endif |
| |
| #include <Internal.h> |
| |
| #define R_MSG_NA _("NaNs produced") |
| #define R_MSG_NONNUM_MATH _("non-numeric argument to mathematical function") |
| |
| #include <Rmath.h> |
| |
| #include <R_ext/Itermacros.h> |
| |
| #include "arithmetic.h" |
| |
| #include <errno.h> |
| |
| #ifdef HAVE_MATHERR |
| |
| /* Override the SVID matherr function: |
| the main difference here is not to print warnings. |
| */ |
| #ifndef __cplusplus |
| int matherr(struct exception *exc) |
| { |
| switch (exc->type) { |
| case DOMAIN: |
| case SING: |
| errno = EDOM; |
| break; |
| case OVERFLOW: |
| errno = ERANGE; |
| break; |
| case UNDERFLOW: |
| exc->retval = 0.0; |
| break; |
| /* |
| There are cases TLOSS and PLOSS which are ignored here. |
| According to the Solaris man page, there are for |
| trigonometric algorithms and not needed for good ones. |
| */ |
| } |
| return 1; |
| } |
| #endif |
| #endif |
| |
| typedef union |
| { |
| double value; |
| unsigned int word[2]; |
| } ieee_double; |
| |
| /* gcc had problems with static const on AIX and Solaris |
| Solaris was for gcc 3.1 and 3.2 under -O2 32-bit on 64-bit kernel */ |
| #ifdef _AIX |
| #define CONST |
| #elif defined(sparc) && defined (__GNUC__) && __GNUC__ == 3 |
| #define CONST |
| #else |
| #define CONST const |
| #endif |
| |
| #ifdef WORDS_BIGENDIAN |
| static CONST int hw = 0; |
| static CONST int lw = 1; |
| #else /* !WORDS_BIGENDIAN */ |
| static CONST int hw = 1; |
| static CONST int lw = 0; |
| #endif /* WORDS_BIGENDIAN */ |
| |
| |
| static double R_ValueOfNA(void) |
| { |
| /* The gcc shipping with Fedora 9 gets this wrong without |
| * the volatile declaration. Thanks to Marc Schwartz. */ |
| volatile ieee_double x; |
| x.word[hw] = 0x7ff00000; |
| x.word[lw] = 1954; |
| return x.value; |
| } |
| |
| int R_IsNA(double x) |
| { |
| if (isnan(x)) { |
| ieee_double y; |
| y.value = x; |
| return (y.word[lw] == 1954); |
| } |
| return 0; |
| } |
| |
| int R_IsNaN(double x) |
| { |
| if (isnan(x)) { |
| ieee_double y; |
| y.value = x; |
| return (y.word[lw] != 1954); |
| } |
| return 0; |
| } |
| |
| /* ISNAN uses isnan, which is undefined by C++ headers |
| This workaround is called only when ISNAN() is used |
| in a user code in a file with __cplusplus defined */ |
| |
| int R_isnancpp(double x) |
| { |
| return (isnan(x)!=0); |
| } |
| |
| |
| /* Mainly for use in packages */ |
| int R_finite(double x) |
| { |
| #ifdef HAVE_WORKING_ISFINITE |
| return isfinite(x); |
| #else |
| return (!isnan(x) & (x != R_PosInf) & (x != R_NegInf)); |
| #endif |
| } |
| |
| |
| /* Arithmetic Initialization */ |
| |
| void attribute_hidden InitArithmetic() |
| { |
| R_NaInt = INT_MIN; |
| R_NaReal = R_ValueOfNA(); |
| // we assume C99, so |
| #ifndef OLD |
| R_NaN = NAN; |
| R_PosInf = INFINITY; |
| R_NegInf = -INFINITY; |
| #else |
| R_NaN = 0.0/R_Zero_Hack; |
| R_PosInf = 1.0/R_Zero_Hack; |
| R_NegInf = -1.0/R_Zero_Hack; |
| #endif |
| } |
| |
| |
| #if HAVE_LONG_DOUBLE && (SIZEOF_LONG_DOUBLE > SIZEOF_DOUBLE) |
| # ifdef __PPC64__ |
| // PowerPC 64 (when gcc has -mlong-double-128) fails constant folding with LDOUBLE |
| # define q_1_eps (1 / LDBL_EPSILON) |
| # else |
| static LDOUBLE q_1_eps = 1 / LDBL_EPSILON; |
| # endif |
| #else |
| static double q_1_eps = 1 / DBL_EPSILON; |
| #endif |
| |
| /* Keep myfmod() and myfloor() in step */ |
| static double myfmod(double x1, double x2) |
| { |
| if (x2 == 0.0) return R_NaN; |
| if(fabs(x2) > q_1_eps && R_FINITE(x1) && fabs(x1) <= fabs(x2)) { |
| return |
| ((x1 < 0 && x2 > 0) || |
| (x1 > 0 && x2 < 0)) |
| ? x1+x2 // differing signs |
| : x1 ; // "same" signs (incl. 0) |
| } |
| double q = x1 / x2; |
| if(R_FINITE(q) && (fabs(q) > q_1_eps)) |
| warning(_("probable complete loss of accuracy in modulus")); |
| LDOUBLE tmp = (LDOUBLE)x1 - floor(q) * (LDOUBLE)x2; |
| return (double) tmp - floorl(tmp/x2) * x2; |
| } |
| |
| static double myfloor(double x1, double x2) |
| { |
| double q = x1 / x2; |
| if (x2 == 0.0 || fabs(q) > q_1_eps || !R_FINITE(q)) |
| return q; |
| if(fabs(q) < 1) |
| return (q < 0) ? -1 |
| : ((x1 < 0 && x2 > 0) || |
| (x1 > 0 && x2 < 0) // differing signs |
| ? -1 : 0); |
| LDOUBLE tmp = (LDOUBLE)x1 - floor(q) * (LDOUBLE)x2; |
| return (double) floor(q) + floorl(tmp/x2); |
| } |
| |
| double R_pow(double x, double y) /* = x ^ y */ |
| { |
| /* squaring is the most common of the specially handled cases so |
| check for it first. */ |
| if(y == 2.0) |
| return x * x; |
| if(x == 1. || y == 0.) |
| return(1.); |
| if(x == 0.) { |
| if(y > 0.) return(0.); |
| else if(y < 0) return(R_PosInf); |
| else return(y); /* NA or NaN, we assert */ |
| } |
| if (R_FINITE(x) && R_FINITE(y)) { |
| /* There was a special case for y == 0.5 here, but |
| gcc 4.3.0 -g -O2 mis-compiled it. Showed up with |
| 100^0.5 as 3.162278, example(pbirthday) failed. */ |
| #ifdef USE_POWL_IN_R_POW |
| // this is used only on 64-bit Windows (so has powl). |
| return powl(x, y); |
| #else |
| return pow(x, y); |
| #endif |
| } |
| if (ISNAN(x) || ISNAN(y)) |
| return(x + y); |
| if(!R_FINITE(x)) { |
| if(x > 0) /* Inf ^ y */ |
| return (y < 0.)? 0. : R_PosInf; |
| else { /* (-Inf) ^ y */ |
| if(R_FINITE(y) && y == floor(y)) /* (-Inf) ^ n */ |
| return (y < 0.) ? 0. : (myfmod(y, 2.) != 0 ? x : -x); |
| } |
| } |
| if(!R_FINITE(y)) { |
| if(x >= 0) { |
| if(y > 0) /* y == +Inf */ |
| return (x >= 1) ? R_PosInf : 0.; |
| else /* y == -Inf */ |
| return (x < 1) ? R_PosInf : 0.; |
| } |
| } |
| return R_NaN; // all other cases: (-Inf)^{+-Inf, non-int}; (neg)^{+-Inf} |
| } |
| |
| double R_pow_di(double x, int n) |
| { |
| double xn = 1.0; |
| |
| if (ISNAN(x)) return x; |
| if (n == NA_INTEGER) return NA_REAL; |
| |
| if (n != 0) { |
| if (!R_FINITE(x)) return R_POW(x, (double)n); |
| |
| Rboolean is_neg = (n < 0); |
| if(is_neg) n = -n; |
| for(;;) { |
| if(n & 01) xn *= x; |
| if(n >>= 1) x *= x; else break; |
| } |
| if(is_neg) xn = 1. / xn; |
| } |
| return xn; |
| } |
| |
| |
| /* General Base Logarithms */ |
| |
| SEXP R_unary(SEXP, SEXP, SEXP); |
| SEXP R_binary(SEXP, SEXP, SEXP, SEXP); |
| static SEXP logical_unary(ARITHOP_TYPE, SEXP, SEXP); |
| static SEXP integer_unary(ARITHOP_TYPE, SEXP, SEXP); |
| static SEXP real_unary(ARITHOP_TYPE, SEXP, SEXP); |
| static SEXP real_binary(ARITHOP_TYPE, SEXP, SEXP); |
| static SEXP integer_binary(ARITHOP_TYPE, SEXP, SEXP, SEXP); |
| |
| #if 0 |
| static int naflag; |
| static SEXP lcall; |
| #endif |
| |
| /* Integer arithmetic support */ |
| |
| /* The tests using integer comparisons are a bit faster than the tests |
| using doubles, but they depend on a two's complement representation |
| (but that is almost universal). The tests that compare results to |
| double's depend on being able to accurately represent all int's as |
| double's. Since int's are almost universally 32 bit that should be |
| OK. */ |
| |
| #ifndef INT_32_BITS |
| /* configure checks whether int is 32 bits. If not this code will |
| need to be rewritten. Since 32 bit ints are pretty much universal, |
| we can worry about writing alternate code when the need arises. |
| To be safe, we signal a compiler error if int is not 32 bits. */ |
| # error code requires that int have 32 bits |
| #endif |
| |
| #define INTEGER_OVERFLOW_WARNING _("NAs produced by integer overflow") |
| |
| #define CHECK_INTEGER_OVERFLOW(call, ans, naflag) do { \ |
| if (naflag) { \ |
| PROTECT(ans); \ |
| warningcall(call, INTEGER_OVERFLOW_WARNING); \ |
| UNPROTECT(1); \ |
| } \ |
| } while(0) |
| |
| #define R_INT_MAX INT_MAX |
| #define R_INT_MIN -INT_MAX |
| // .. relying on fact that NA_INTEGER is outside of these |
| |
| static R_INLINE int R_integer_plus(int x, int y, Rboolean *pnaflag) |
| { |
| if (x == NA_INTEGER || y == NA_INTEGER) |
| return NA_INTEGER; |
| |
| if (((y > 0) && (x > (R_INT_MAX - y))) || |
| ((y < 0) && (x < (R_INT_MIN - y)))) { |
| if (pnaflag != NULL) |
| *pnaflag = TRUE; |
| return NA_INTEGER; |
| } |
| return x + y; |
| } |
| |
| static R_INLINE int R_integer_minus(int x, int y, Rboolean *pnaflag) |
| { |
| if (x == NA_INTEGER || y == NA_INTEGER) |
| return NA_INTEGER; |
| |
| if (((y < 0) && (x > (R_INT_MAX + y))) || |
| ((y > 0) && (x < (R_INT_MIN + y)))) { |
| if (pnaflag != NULL) |
| *pnaflag = TRUE; |
| return NA_INTEGER; |
| } |
| return x - y; |
| } |
| |
| #define GOODIPROD(x, y, z) ((double) (x) * (double) (y) == (z)) |
| static R_INLINE int R_integer_times(int x, int y, Rboolean *pnaflag) |
| { |
| if (x == NA_INTEGER || y == NA_INTEGER) |
| return NA_INTEGER; |
| else { |
| int z = x * y; // UBSAN will warn if this overflows (happens in bda) |
| if (GOODIPROD(x, y, z) && z != NA_INTEGER) |
| return z; |
| else { |
| if (pnaflag != NULL) |
| *pnaflag = TRUE; |
| return NA_INTEGER; |
| } |
| } |
| } |
| |
| static R_INLINE double R_integer_divide(int x, int y) |
| { |
| if (x == NA_INTEGER || y == NA_INTEGER) |
| return NA_REAL; |
| else |
| return (double) x / (double) y; |
| } |
| |
| static R_INLINE SEXP ScalarValue1(SEXP x) |
| { |
| if (NO_REFERENCES(x)) |
| return x; |
| else |
| return allocVector(TYPEOF(x), 1); |
| } |
| |
| static R_INLINE SEXP ScalarValue2(SEXP x, SEXP y) |
| { |
| if (NO_REFERENCES(x)) |
| return x; |
| else if (NO_REFERENCES(y)) |
| return y; |
| else |
| return allocVector(TYPEOF(x), 1); |
| } |
| |
| /* Unary and Binary Operators */ |
| |
| SEXP attribute_hidden do_arith(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP ans, arg1, arg2; |
| int argc; |
| |
| if (args == R_NilValue) |
| argc = 0; |
| else if (CDR(args) == R_NilValue) |
| argc = 1; |
| else if (CDDR(args) == R_NilValue) |
| argc = 2; |
| else |
| argc = length(args); |
| arg1 = CAR(args); |
| arg2 = CADR(args); |
| |
| if (ATTRIB(arg1) != R_NilValue || ATTRIB(arg2) != R_NilValue) { |
| if (DispatchGroup("Ops", call, op, args, env, &ans)) |
| return ans; |
| } |
| else if (argc == 2) { |
| /* Handle some scaler operations immediately */ |
| if (IS_SCALAR(arg1, REALSXP)) { |
| double x1 = SCALAR_DVAL(arg1); |
| if (IS_SCALAR(arg2, REALSXP)) { |
| double x2 = SCALAR_DVAL(arg2); |
| ans = ScalarValue2(arg1, arg2); |
| switch (PRIMVAL(op)) { |
| case PLUSOP: SET_SCALAR_DVAL(ans, x1 + x2); return ans; |
| case MINUSOP: SET_SCALAR_DVAL(ans, x1 - x2); return ans; |
| case TIMESOP: SET_SCALAR_DVAL(ans, x1 * x2); return ans; |
| case DIVOP: SET_SCALAR_DVAL(ans, x1 / x2); return ans; |
| } |
| } |
| else if (IS_SCALAR(arg2, INTSXP)) { |
| int i2 = SCALAR_IVAL(arg2); |
| double x2 = i2 != NA_INTEGER ? (double) i2 : NA_REAL; |
| ans = ScalarValue1(arg1); |
| switch (PRIMVAL(op)) { |
| case PLUSOP: SET_SCALAR_DVAL(ans, x1 + x2); return ans; |
| case MINUSOP: SET_SCALAR_DVAL(ans, x1 - x2); return ans; |
| case TIMESOP: SET_SCALAR_DVAL(ans, x1 * x2); return ans; |
| case DIVOP: SET_SCALAR_DVAL(ans, x1 / x2); return ans; |
| } |
| } |
| } |
| else if (IS_SCALAR(arg1, INTSXP)) { |
| int i1 = SCALAR_IVAL(arg1); |
| if (IS_SCALAR(arg2, REALSXP)) { |
| double x1 = i1 != NA_INTEGER ? (double) i1 : NA_REAL; |
| double x2 = SCALAR_DVAL(arg2); |
| ans = ScalarValue1(arg2); |
| switch (PRIMVAL(op)) { |
| case PLUSOP: SET_SCALAR_DVAL(ans, x1 + x2); return ans; |
| case MINUSOP: SET_SCALAR_DVAL(ans, x1 - x2); return ans; |
| case TIMESOP: SET_SCALAR_DVAL(ans, x1 * x2); return ans; |
| case DIVOP: SET_SCALAR_DVAL(ans, x1 / x2); return ans; |
| } |
| } |
| else if (IS_SCALAR(arg2, INTSXP)) { |
| Rboolean naflag = FALSE; |
| int i2 = SCALAR_IVAL(arg2); |
| switch (PRIMVAL(op)) { |
| case PLUSOP: |
| ans = ScalarValue2(arg1, arg2); |
| SET_SCALAR_IVAL(ans, R_integer_plus(i1, i2, &naflag)); |
| CHECK_INTEGER_OVERFLOW(call, ans, naflag); |
| return ans; |
| case MINUSOP: |
| ans = ScalarValue2(arg1, arg2); |
| SET_SCALAR_IVAL(ans, R_integer_minus(i1, i2, &naflag)); |
| CHECK_INTEGER_OVERFLOW(call, ans, naflag); |
| return ans; |
| case TIMESOP: |
| ans = ScalarValue2(arg1, arg2); |
| SET_SCALAR_IVAL(ans, R_integer_times(i1, i2, &naflag)); |
| CHECK_INTEGER_OVERFLOW(call, ans, naflag); |
| return ans; |
| case DIVOP: |
| return ScalarReal(R_integer_divide(i1, i2)); |
| } |
| } |
| } |
| } |
| else if (argc == 1) { |
| if (IS_SCALAR(arg1, REALSXP)) { |
| switch(PRIMVAL(op)) { |
| case PLUSOP: return(arg1); |
| case MINUSOP: |
| ans = ScalarValue1(arg1); |
| SET_SCALAR_DVAL(ans, -SCALAR_DVAL(arg1)); |
| return ans; |
| } |
| } |
| else if (IS_SCALAR(arg1, INTSXP)) { |
| int ival; |
| switch(PRIMVAL(op)) { |
| case PLUSOP: return(arg1); |
| case MINUSOP: |
| ival = SCALAR_IVAL(arg1); |
| ans = ScalarValue1(arg1); |
| SET_SCALAR_IVAL(ans, ival == NA_INTEGER ? NA_INTEGER : -ival); |
| return ans; |
| } |
| } |
| } |
| |
| if (argc == 2) |
| return R_binary(call, op, arg1, arg2); |
| else if (argc == 1) |
| return R_unary(call, op, arg1); |
| else |
| errorcall(call,_("operator needs one or two arguments")); |
| return ans; /* never used; to keep -Wall happy */ |
| } |
| |
| #define COERCE_IF_NEEDED(v, tp, vpi) do { \ |
| if (TYPEOF(v) != (tp)) { \ |
| int __vo__ = OBJECT(v); \ |
| REPROTECT(v = coerceVector(v, (tp)), vpi); \ |
| if (__vo__) SET_OBJECT(v, 1); \ |
| } \ |
| } while (0) |
| |
| #define FIXUP_NULL_AND_CHECK_TYPES(v, vpi) do { \ |
| switch (TYPEOF(v)) { \ |
| case NILSXP: REPROTECT(v = allocVector(INTSXP,0), vpi); break; \ |
| case CPLXSXP: case REALSXP: case INTSXP: case LGLSXP: break; \ |
| default: errorcall(call, _("non-numeric argument to binary operator")); \ |
| } \ |
| } while (0) |
| |
| SEXP attribute_hidden R_binary(SEXP call, SEXP op, SEXP x, SEXP y) |
| { |
| Rboolean xattr, yattr, xarray, yarray, xts, yts, xS4, yS4; |
| PROTECT_INDEX xpi, ypi; |
| ARITHOP_TYPE oper = (ARITHOP_TYPE) PRIMVAL(op); |
| int nprotect = 2; /* x and y */ |
| |
| |
| PROTECT_WITH_INDEX(x, &xpi); |
| PROTECT_WITH_INDEX(y, &ypi); |
| |
| FIXUP_NULL_AND_CHECK_TYPES(x, xpi); |
| FIXUP_NULL_AND_CHECK_TYPES(y, ypi); |
| |
| R_xlen_t |
| nx = XLENGTH(x), |
| ny = XLENGTH(y); |
| if (ATTRIB(x) != R_NilValue) { |
| xattr = TRUE; |
| xarray = isArray(x); |
| xts = isTs(x); |
| xS4 = isS4(x); |
| } |
| else xattr = xarray = xts = xS4 = FALSE; |
| if (ATTRIB(y) != R_NilValue) { |
| yattr = TRUE; |
| yarray = isArray(y); |
| yts = isTs(y); |
| yS4 = isS4(y); |
| } |
| else yattr = yarray = yts = yS4 = FALSE; |
| |
| #define R_ARITHMETIC_ARRAY_1_SPECIAL |
| |
| #ifdef R_ARITHMETIC_ARRAY_1_SPECIAL |
| /* If either x or y is a matrix with length 1 and the other is a |
| vector of a different length, we want to coerce the matrix to be a vector. |
| Do we want to? We don't do it! BDR 2004-03-06 |
| |
| From 3.4.0 (Sep. 2016), this signals a warning, |
| and in the future we will disable these 2 clauses, |
| so it will give an error. |
| */ |
| |
| /* FIXME: Danger Will Robinson. |
| * ----- We might be trashing arguments here. |
| */ |
| if (xarray != yarray) { |
| if (xarray && nx==1 && ny!=1) { |
| if(ny != 0) |
| warningcall(call, _( |
| "Recycling array of length 1 in array-vector arithmetic is deprecated.\n\ |
| Use c() or as.vector() instead.\n")); |
| REPROTECT(x = duplicate(x), xpi); |
| setAttrib(x, R_DimSymbol, R_NilValue); |
| } |
| if (yarray && ny==1 && nx!=1) { |
| if(nx != 0) |
| warningcall(call, _( |
| "Recycling array of length 1 in vector-array arithmetic is deprecated.\n\ |
| Use c() or as.vector() instead.\n")); |
| REPROTECT(y = duplicate(y), ypi); |
| setAttrib(y, R_DimSymbol, R_NilValue); |
| } |
| } |
| #endif |
| |
| SEXP dims, xnames, ynames; |
| if (xarray || yarray) { |
| /* if one is a length-atleast-1-array and the |
| * other is a length-0 *non*array, then do not use array treatment */ |
| if (xarray && yarray) { |
| if (!conformable(x, y)) |
| errorcall(call, _("non-conformable arrays")); |
| PROTECT(dims = getAttrib(x, R_DimSymbol)); nprotect++; |
| } |
| else if (xarray && (ny != 0 || nx == 0)) { |
| PROTECT(dims = getAttrib(x, R_DimSymbol)); nprotect++; |
| } |
| else if (yarray && (nx != 0 || ny == 0)) { |
| PROTECT(dims = getAttrib(y, R_DimSymbol)); nprotect++; |
| } else |
| dims = R_NilValue; |
| if (xattr) { |
| PROTECT(xnames = getAttrib(x, R_DimNamesSymbol)); |
| nprotect++; |
| } |
| else xnames = R_NilValue; |
| if (yattr) { |
| PROTECT(ynames = getAttrib(y, R_DimNamesSymbol)); |
| nprotect++; |
| } |
| else ynames = R_NilValue; |
| } |
| else { |
| dims = R_NilValue; |
| if (xattr) { |
| PROTECT(xnames = getAttrib(x, R_NamesSymbol)); |
| nprotect++; |
| } |
| else xnames = R_NilValue; |
| if (yattr) { |
| PROTECT(ynames = getAttrib(y, R_NamesSymbol)); |
| nprotect++; |
| } |
| else ynames = R_NilValue; |
| } |
| |
| SEXP klass = NULL, tsp = NULL; // -Wall |
| if (xts || yts) { |
| if (xts && yts) { |
| /* could check ts conformance here */ |
| PROTECT(tsp = getAttrib(x, R_TspSymbol)); |
| PROTECT(klass = getAttrib(x, R_ClassSymbol)); |
| } |
| else if (xts) { |
| if (nx < ny) |
| ErrorMessage(call, ERROR_TSVEC_MISMATCH); |
| PROTECT(tsp = getAttrib(x, R_TspSymbol)); |
| PROTECT(klass = getAttrib(x, R_ClassSymbol)); |
| } |
| else { /* (yts) */ |
| if (ny < nx) |
| ErrorMessage(call, ERROR_TSVEC_MISMATCH); |
| PROTECT(tsp = getAttrib(y, R_TspSymbol)); |
| PROTECT(klass = getAttrib(y, R_ClassSymbol)); |
| } |
| nprotect += 2; |
| } |
| |
| if (nx > 0 && ny > 0 && |
| ((nx > ny) ? nx % ny : ny % nx) != 0) // mismatch |
| warningcall(call, |
| _("longer object length is not a multiple of shorter object length")); |
| |
| SEXP val; |
| /* need to preserve object here, as *_binary copies class attributes */ |
| if (TYPEOF(x) == CPLXSXP || TYPEOF(y) == CPLXSXP) { |
| COERCE_IF_NEEDED(x, CPLXSXP, xpi); |
| COERCE_IF_NEEDED(y, CPLXSXP, ypi); |
| val = complex_binary(oper, x, y); |
| } |
| else if (TYPEOF(x) == REALSXP || TYPEOF(y) == REALSXP) { |
| /* real_binary can handle REALSXP or INTSXP operand, but not LGLSXP. */ |
| /* Can get a LGLSXP. In base-Ex.R on 24 Oct '06, got 8 of these. */ |
| if (TYPEOF(x) != INTSXP) COERCE_IF_NEEDED(x, REALSXP, xpi); |
| if (TYPEOF(y) != INTSXP) COERCE_IF_NEEDED(y, REALSXP, ypi); |
| val = real_binary(oper, x, y); |
| } |
| else val = integer_binary(oper, x, y, call); |
| |
| /* quick return if there are no attributes */ |
| if (! xattr && ! yattr) { |
| UNPROTECT(nprotect); |
| return val; |
| } |
| |
| PROTECT(val); |
| nprotect++; |
| |
| if (dims != R_NilValue) { |
| setAttrib(val, R_DimSymbol, dims); |
| if (xnames != R_NilValue) |
| setAttrib(val, R_DimNamesSymbol, xnames); |
| else if (ynames != R_NilValue) |
| setAttrib(val, R_DimNamesSymbol, ynames); |
| } |
| else { |
| if (XLENGTH(val) == xlength(xnames)) |
| setAttrib(val, R_NamesSymbol, xnames); |
| else if (XLENGTH(val) == xlength(ynames)) |
| setAttrib(val, R_NamesSymbol, ynames); |
| } |
| |
| if (xts || yts) { /* must set *after* dims! */ |
| setAttrib(val, R_TspSymbol, tsp); |
| setAttrib(val, R_ClassSymbol, klass); |
| } |
| |
| if(xS4 || yS4) { /* Only set the bit: no method defined! */ |
| val = asS4(val, TRUE, TRUE); |
| } |
| UNPROTECT(nprotect); |
| return val; |
| } |
| |
| SEXP attribute_hidden R_unary(SEXP call, SEXP op, SEXP s1) |
| { |
| ARITHOP_TYPE operation = (ARITHOP_TYPE) PRIMVAL(op); |
| switch (TYPEOF(s1)) { |
| case LGLSXP: |
| return logical_unary(operation, s1, call); |
| case INTSXP: |
| return integer_unary(operation, s1, call); |
| case REALSXP: |
| return real_unary(operation, s1, call); |
| case CPLXSXP: |
| return complex_unary(operation, s1, call); |
| default: |
| errorcall(call, _("invalid argument to unary operator")); |
| } |
| return s1; /* never used; to keep -Wall happy */ |
| } |
| |
| static SEXP logical_unary(ARITHOP_TYPE code, SEXP s1, SEXP call) |
| { |
| R_xlen_t n = XLENGTH(s1); |
| SEXP ans = PROTECT(allocVector(INTSXP, n)); |
| SEXP names = PROTECT(getAttrib(s1, R_NamesSymbol)); |
| SEXP dim = PROTECT(getAttrib(s1, R_DimSymbol)); |
| SEXP dimnames = PROTECT(getAttrib(s1, R_DimNamesSymbol)); |
| if(names != R_NilValue) setAttrib(ans, R_NamesSymbol, names); |
| if(dim != R_NilValue) setAttrib(ans, R_DimSymbol, dim); |
| if(dimnames != R_NilValue) setAttrib(ans, R_DimNamesSymbol, dimnames); |
| UNPROTECT(3); |
| |
| int *pa = INTEGER(ans); |
| const int *px = LOGICAL_RO(s1); |
| |
| switch (code) { |
| case PLUSOP: |
| for (R_xlen_t i = 0; i < n; i++) pa[i] = px[i]; |
| break; |
| case MINUSOP: |
| for (R_xlen_t i = 0; i < n; i++) { |
| int x = px[i]; |
| pa[i] = (x == NA_INTEGER) ? |
| NA_INTEGER : ((x == 0.0) ? 0 : -x); |
| } |
| break; |
| default: |
| errorcall(call, _("invalid unary operator")); |
| } |
| UNPROTECT(1); |
| return ans; |
| } |
| |
| static SEXP integer_unary(ARITHOP_TYPE code, SEXP s1, SEXP call) |
| { |
| R_xlen_t i, n; |
| SEXP ans; |
| |
| switch (code) { |
| case PLUSOP: |
| return s1; |
| case MINUSOP: |
| ans = NO_REFERENCES(s1) ? s1 : duplicate(s1); |
| int *pa = INTEGER(ans); |
| const int *px = INTEGER_RO(s1); |
| n = XLENGTH(s1); |
| for (i = 0; i < n; i++) { |
| int x = px[i]; |
| pa[i] = (x == NA_INTEGER) ? |
| NA_INTEGER : ((x == 0.0) ? 0 : -x); |
| } |
| return ans; |
| default: |
| errorcall(call, _("invalid unary operator")); |
| } |
| return s1; /* never used; to keep -Wall happy */ |
| } |
| |
| static SEXP real_unary(ARITHOP_TYPE code, SEXP s1, SEXP lcall) |
| { |
| R_xlen_t i, n; |
| SEXP ans; |
| |
| switch (code) { |
| case PLUSOP: return s1; |
| case MINUSOP: |
| ans = NO_REFERENCES(s1) ? s1 : duplicate(s1); |
| double *pa = REAL(ans); |
| const double *px = REAL_RO(s1); |
| n = XLENGTH(s1); |
| for (i = 0; i < n; i++) |
| pa[i] = -px[i]; |
| return ans; |
| default: |
| errorcall(lcall, _("invalid unary operator")); |
| } |
| return s1; /* never used; to keep -Wall happy */ |
| } |
| |
| static SEXP integer_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2, SEXP lcall) |
| { |
| R_xlen_t i, i1, i2, n, n1, n2; |
| int x1, x2; |
| SEXP ans; |
| Rboolean naflag = FALSE; |
| |
| n1 = XLENGTH(s1); |
| n2 = XLENGTH(s2); |
| /* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */ |
| if (n1 == 0 || n2 == 0) n = 0; else n = (n1 > n2) ? n1 : n2; |
| |
| if (code == DIVOP || code == POWOP) |
| ans = allocVector(REALSXP, n); |
| else |
| ans = R_allocOrReuseVector(s1, s2, INTSXP, n); |
| if (n == 0) return(ans); |
| PROTECT(ans); |
| |
| switch (code) { |
| case PLUSOP: |
| { |
| int *pa = INTEGER(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| x1 = px1[i1]; |
| x2 = px2[i2]; |
| pa[i] = R_integer_plus(x1, x2, &naflag); |
| }); |
| if (naflag) |
| warningcall(lcall, INTEGER_OVERFLOW_WARNING); |
| } |
| break; |
| case MINUSOP: |
| { |
| int *pa = INTEGER(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| x1 = px1[i1]; |
| x2 = px2[i2]; |
| pa[i] = R_integer_minus(x1, x2, &naflag); |
| }); |
| if (naflag) |
| warningcall(lcall, INTEGER_OVERFLOW_WARNING); |
| } |
| break; |
| case TIMESOP: |
| { |
| int *pa = INTEGER(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| x1 = px1[i1]; |
| x2 = px2[i2]; |
| pa[i] = R_integer_times(x1, x2, &naflag); |
| }); |
| if (naflag) |
| warningcall(lcall, INTEGER_OVERFLOW_WARNING); |
| } |
| break; |
| case DIVOP: |
| { |
| double *pa = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| x1 = px1[i1]; |
| x2 = px2[i2]; |
| pa[i] = R_integer_divide(x1, x2); |
| }); |
| } |
| break; |
| case POWOP: |
| { |
| double *pa = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| if((x1 = px1[i1]) == 1 || (x2 = px2[i2]) == 0) |
| pa[i] = 1.; |
| else if (x1 == NA_INTEGER || x2 == NA_INTEGER) |
| pa[i] = NA_REAL; |
| else |
| pa[i] = R_POW((double) x1, (double) x2); |
| }); |
| } |
| break; |
| case MODOP: |
| { |
| int *pa = INTEGER(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| x1 = px1[i1]; |
| x2 = px2[i2]; |
| if (x1 == NA_INTEGER || x2 == NA_INTEGER || x2 == 0) |
| pa[i] = NA_INTEGER; |
| else { |
| pa[i] = /* till 0.63.2: x1 % x2 */ |
| (x1 >= 0 && x2 > 0) ? x1 % x2 : |
| (int)myfmod((double)x1,(double)x2); |
| } |
| }); |
| } |
| break; |
| case IDIVOP: |
| { |
| int *pa = INTEGER(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, { |
| x1 = px1[i1]; |
| x2 = px2[i2]; |
| /* This had x %/% 0 == 0 prior to 2.14.1, but |
| it seems conventionally to be undefined */ |
| if (x1 == NA_INTEGER || x2 == NA_INTEGER || x2 == 0) |
| pa[i] = NA_INTEGER; |
| else |
| pa[i] = (int) floor((double)x1 / (double)x2); |
| }); |
| } |
| break; |
| } |
| UNPROTECT(1); |
| |
| /* quick return if there are no attributes */ |
| if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue) |
| return ans; |
| |
| /* Copy attributes from longer argument. */ |
| |
| if (ans != s2 && n == n2 && ATTRIB(s2) != R_NilValue) |
| copyMostAttrib(s2, ans); |
| if (ans != s1 && n == n1 && ATTRIB(s1) != R_NilValue) |
| copyMostAttrib(s1, ans); /* Done 2nd so s1's attrs overwrite s2's */ |
| |
| return ans; |
| } |
| |
| #define R_INTEGER(x) (double) ((x) == NA_INTEGER ? NA_REAL : (x)) |
| |
| static SEXP real_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2) |
| { |
| R_xlen_t i, i1, i2, n, n1, n2; |
| SEXP ans; |
| |
| /* Note: "s1" and "s2" are protected above. */ |
| n1 = XLENGTH(s1); |
| n2 = XLENGTH(s2); |
| |
| /* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */ |
| if (n1 == 0 || n2 == 0) return(allocVector(REALSXP, 0)); |
| |
| n = (n1 > n2) ? n1 : n2; |
| PROTECT(ans = R_allocOrReuseVector(s1, s2, REALSXP, n)); |
| |
| switch (code) { |
| case PLUSOP: |
| if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) { |
| double *da = REAL(ans); |
| const double *dx = REAL_RO(s1); |
| const double *dy = REAL_RO(s2); |
| if (n2 == 1) { |
| double tmp = dy[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] + tmp;); |
| } |
| else if (n1 == 1) { |
| double tmp = dx[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp + dy[i];); |
| } |
| else if (n1 == n2) |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] + dy[i];); |
| else |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = dx[i1] + dy[i2];); |
| } |
| else if(TYPEOF(s1) == INTSXP ) { |
| double *da = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = R_INTEGER(px1[i1]) + px2[i2];); |
| } |
| else if(TYPEOF(s2) == INTSXP ) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = px1[i1] + R_INTEGER(px2[i2]);); |
| } |
| break; |
| case MINUSOP: |
| if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) { |
| double *da = REAL(ans); |
| const double *dx = REAL_RO(s1); |
| const double *dy = REAL_RO(s2); |
| if (n2 == 1) { |
| double tmp = dy[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] - tmp;); |
| } |
| else if (n1 == 1) { |
| double tmp = dx[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp - dy[i];); |
| } |
| else if (n1 == n2) |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] - dy[i];); |
| else |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = dx[i1] - dy[i2];); |
| } |
| else if(TYPEOF(s1) == INTSXP ) { |
| double *da = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = R_INTEGER(px1[i1]) - px2[i2];); |
| } |
| else if(TYPEOF(s2) == INTSXP ) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = px1[i1] - R_INTEGER(px2[i2]);); |
| } |
| break; |
| case TIMESOP: |
| if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) { |
| double *da = REAL(ans); |
| const double *dx = REAL_RO(s1); |
| const double *dy = REAL_RO(s2); |
| if (n2 == 1) { |
| double tmp = dy[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] * tmp;); |
| } |
| else if (n1 == 1) { |
| double tmp = dx[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp * dy[i];); |
| } |
| else if (n1 == n2) |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] * dy[i];); |
| else |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = dx[i1] * dy[i2];); |
| } |
| else if(TYPEOF(s1) == INTSXP ) { |
| double *da = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = R_INTEGER(px1[i1]) * px2[i2];); |
| } |
| else if(TYPEOF(s2) == INTSXP ) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = px1[i1] * R_INTEGER(px2[i2]);); |
| } |
| break; |
| case DIVOP: |
| if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) { |
| double *da = REAL(ans); |
| const double *dx = REAL_RO(s1); |
| const double *dy = REAL_RO(s2); |
| if (n2 == 1) { |
| double tmp = dy[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] / tmp;); |
| } |
| else if (n1 == 1) { |
| double tmp = dx[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp / dy[i];); |
| } |
| else if (n1 == n2) |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] / dy[i];); |
| else |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = dx[i1] / dy[i2];); |
| } |
| else if(TYPEOF(s1) == INTSXP ) { |
| double *da = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = R_INTEGER(px1[i1]) / px2[i2];); |
| } |
| else if(TYPEOF(s2) == INTSXP ) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = px1[i1] / R_INTEGER(px2[i2]);); |
| } |
| break; |
| case POWOP: |
| if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) { |
| double *da = REAL(ans); |
| const double *dx = REAL_RO(s1); |
| const double *dy = REAL_RO(s2); |
| if (n2 == 1) { |
| double tmp = dy[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = R_POW(dx[i], tmp);); |
| } |
| else if (n1 == 1) { |
| double tmp = dx[0]; |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = R_POW(tmp, dy[i]);); |
| } |
| else if (n1 == n2) |
| R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = R_POW(dx[i], dy[i]);); |
| else |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = R_POW(dx[i1], dy[i2]);); |
| } |
| else if(TYPEOF(s1) == INTSXP ) { |
| double *da = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = R_POW( R_INTEGER(px1[i1]), px2[i2]);); |
| } |
| else if(TYPEOF(s2) == INTSXP ) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = R_POW(px1[i1], R_INTEGER(px2[i2]));); |
| } |
| break; |
| case MODOP: |
| if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = myfmod(px1[i1], px2[i2]);); |
| } |
| else if(TYPEOF(s1) == INTSXP ) { |
| double *da = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = myfmod(R_INTEGER(px1[i1]), px2[i2]);); |
| } |
| else if(TYPEOF(s2) == INTSXP ) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = myfmod(px1[i1], R_INTEGER(px2[i2]));); |
| } |
| break; |
| case IDIVOP: |
| if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = myfloor(px1[i1], px2[i2]);); |
| } |
| else if(TYPEOF(s1) == INTSXP ) { |
| double *da = REAL(ans); |
| const int *px1 = INTEGER_RO(s1); |
| const double *px2 = REAL_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = myfloor(R_INTEGER(px1[i1]), px2[i2]);); |
| } |
| else if(TYPEOF(s2) == INTSXP ) { |
| double *da = REAL(ans); |
| const double *px1 = REAL_RO(s1); |
| const int *px2 = INTEGER_RO(s2); |
| MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, |
| da[i] = myfloor(px1[i1], R_INTEGER(px2[i2]));); |
| } |
| break; |
| } |
| UNPROTECT(1); |
| |
| /* quick return if there are no attributes */ |
| if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue) |
| return ans; |
| |
| /* Copy attributes from longer argument. */ |
| |
| if (ans != s2 && n == n2 && ATTRIB(s2) != R_NilValue) |
| copyMostAttrib(s2, ans); |
| if (ans != s1 && n == n1 && ATTRIB(s1) != R_NilValue) |
| copyMostAttrib(s1, ans); /* Done 2nd so s1's attrs overwrite s2's */ |
| |
| return ans; |
| } |
| |
| |
| /* Mathematical Functions of One Argument */ |
| |
| static SEXP math1(SEXP sa, double(*f)(double), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, n; |
| int naflag; |
| |
| if (!isNumeric(sa)) |
| errorcall(lcall, R_MSG_NONNUM_MATH); |
| |
| n = XLENGTH(sa); |
| /* coercion can lose the object bit */ |
| PROTECT(sa = coerceVector(sa, REALSXP)); |
| PROTECT(sy = NO_REFERENCES(sa) ? sa : allocVector(REALSXP, n)); |
| const double *a = REAL_RO(sa); |
| double *y = REAL(sy); |
| naflag = 0; |
| for (i = 0; i < n; i++) { |
| double x = a[i]; /* in case y == a */ |
| /* This code assumes that ISNAN(x) implies ISNAN(f(x)), so we |
| only need to check ISNAN(x) if ISNAN(f(x)) is true. */ |
| y[i] = f(x); |
| if (ISNAN(y[i])) { |
| if (ISNAN(x)) |
| y[i] = x; /* make sure the incoming NaN is preserved */ |
| else |
| naflag = 1; |
| } |
| } |
| /* These are primitives, so need to use the call */ |
| if(naflag) warningcall(lcall, R_MSG_NA); |
| |
| if (sa != sy && ATTRIB(sa) != R_NilValue) |
| SHALLOW_DUPLICATE_ATTRIB(sy, sa); |
| UNPROTECT(2); |
| return sy; |
| } |
| |
| |
| SEXP attribute_hidden do_math1(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP s; |
| |
| checkArity(op, args); |
| check1arg(args, call, "x"); |
| |
| if (DispatchGroup("Math", call, op, args, env, &s)) |
| return s; |
| |
| if (isComplex(CAR(args))) |
| return complex_math1(call, op, args, env); |
| |
| #define MATH1(x) math1(CAR(args), x, call); |
| switch (PRIMVAL(op)) { |
| case 1: return MATH1(floor); |
| case 2: return MATH1(ceil); |
| case 3: return MATH1(sqrt); |
| case 4: return MATH1(sign); |
| /* case 5: return MATH1(trunc); separate from 2.6.0 */ |
| |
| case 10: return MATH1(exp); |
| case 11: return MATH1(expm1); |
| case 12: return MATH1(log1p); |
| |
| case 20: return MATH1(cos); |
| case 21: return MATH1(sin); |
| case 22: return MATH1(tan); |
| case 23: return MATH1(acos); |
| case 24: return MATH1(asin); |
| case 25: return MATH1(atan); |
| |
| case 30: return MATH1(cosh); |
| case 31: return MATH1(sinh); |
| case 32: return MATH1(tanh); |
| case 33: return MATH1(acosh); |
| case 34: return MATH1(asinh); |
| case 35: return MATH1(atanh); |
| |
| case 40: return MATH1(lgammafn); |
| case 41: return MATH1(gammafn); |
| |
| case 42: return MATH1(digamma); |
| case 43: return MATH1(trigamma); |
| /* case 44: return MATH1(tetragamma); |
| case 45: return MATH1(pentagamma); |
| removed in 2.0.0 -- rather use Math2's psigamma() |
| |
| case 46: return MATH1(Rf_gamma_cody); removed in 2.8.0 |
| */ |
| case 47: return MATH1(cospi); |
| case 48: return MATH1(sinpi); |
| #if defined(HAVE_TANPI) || defined(HAVE___TANPI) |
| case 49: return MATH1(Rtanpi); |
| #else |
| case 49: return MATH1(tanpi); |
| #endif |
| |
| default: |
| errorcall(call, _("unimplemented real function of 1 argument")); |
| } |
| return s; /* never used; to keep -Wall happy */ |
| } |
| |
| /* methods are allowed to have more than one arg */ |
| SEXP attribute_hidden do_trunc(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP s; |
| if (DispatchGroup("Math", call, op, args, env, &s)) |
| return s; |
| // checkArity(op, args); /* is -1 in names.c */ |
| check1arg(args, call, "x"); |
| if (isComplex(CAR(args))) |
| errorcall(call, _("unimplemented complex function")); |
| return math1(CAR(args), trunc, call); |
| } |
| |
| /* |
| Note that this is slightly different from the do_math1 set, |
| both for integer/logical inputs and what it dispatches to for complex ones. |
| */ |
| |
| SEXP attribute_hidden do_abs(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP x, s = R_NilValue /* -Wall */; |
| |
| checkArity(op, args); |
| check1arg(args, call, "x"); |
| x = CAR(args); |
| |
| if (DispatchGroup("Math", call, op, args, env, &s)) |
| return s; |
| |
| if (isInteger(x) || isLogical(x)) { |
| /* integer or logical ==> return integer, |
| factor was covered by Math.factor. */ |
| R_xlen_t i, n = XLENGTH(x); |
| s = (NO_REFERENCES(x) && TYPEOF(x) == INTSXP) ? |
| x : allocVector(INTSXP, n); |
| PROTECT(s); |
| /* Note: relying on INTEGER(.) === LOGICAL(.) : */ |
| int *pa = INTEGER(s); |
| const int *px = INTEGER_RO(x); |
| for(i = 0 ; i < n ; i++) { |
| int xi = px[i]; |
| pa[i] = (xi == NA_INTEGER) ? xi : abs(xi); |
| } |
| } else if (TYPEOF(x) == REALSXP) { |
| R_xlen_t i, n = XLENGTH(x); |
| PROTECT(s = NO_REFERENCES(x) ? x : allocVector(REALSXP, n)); |
| double *pa = REAL(s); |
| const double *px = REAL_RO(x); |
| for(i = 0 ; i < n ; i++) |
| pa[i] = fabs(px[i]); |
| } else if (isComplex(x)) { |
| SET_TAG(args, R_NilValue); /* cmathfuns want "z"; we might have "x" PR#16047 */ |
| return do_cmathfuns(call, op, args, env); |
| } else |
| errorcall(call, R_MSG_NONNUM_MATH); |
| |
| if (x != s && ATTRIB(x) != R_NilValue) |
| SHALLOW_DUPLICATE_ATTRIB(s, x); |
| UNPROTECT(1); |
| return s; |
| } |
| |
| /* Mathematical Functions of Two Numeric Arguments (plus 1 int) */ |
| |
| /* math2_1 and math2_2 and related can be removed once the byte |
| compiler knows how to optimize to .External rather than |
| .Internal */ |
| |
| #define if_NA_Math2_set(y,a,b) \ |
| if (ISNA (a) || ISNA (b)) y = NA_REAL; \ |
| else if (ISNAN(a) || ISNAN(b)) y = R_NaN; |
| |
| static SEXP math2(SEXP sa, SEXP sb, double (*f)(double, double), |
| SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, n, na, nb; |
| double ai, bi, *y; |
| const double *a, *b; |
| int naflag; |
| |
| if (!isNumeric(sa) || !isNumeric(sb)) |
| errorcall(lcall, R_MSG_NONNUM_MATH); |
| |
| /* for 0-length a we want the attributes of a, not those of b |
| as no recycling will occur */ |
| #define SETUP_Math2 \ |
| na = XLENGTH(sa); \ |
| nb = XLENGTH(sb); \ |
| if ((na == 0) || (nb == 0)) { \ |
| PROTECT(sy = allocVector(REALSXP, 0)); \ |
| if (na == 0) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \ |
| UNPROTECT(1); \ |
| return(sy); \ |
| } \ |
| n = (na < nb) ? nb : na; \ |
| PROTECT(sa = coerceVector(sa, REALSXP)); \ |
| PROTECT(sb = coerceVector(sb, REALSXP)); \ |
| PROTECT(sy = allocVector(REALSXP, n)); \ |
| a = REAL_RO(sa); \ |
| b = REAL_RO(sb); \ |
| y = REAL(sy); \ |
| naflag = 0 |
| |
| SETUP_Math2; |
| |
| MOD_ITERATE2(n, na, nb, i, ia, ib, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| if_NA_Math2_set(y[i], ai, bi) |
| else { |
| y[i] = f(ai, bi); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| |
| #define FINISH_Math2 \ |
| if(naflag) warning(R_MSG_NA); \ |
| if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \ |
| else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \ |
| UNPROTECT(3) |
| |
| FINISH_Math2; |
| |
| return sy; |
| } /* math2() */ |
| |
| static SEXP math2_1(SEXP sa, SEXP sb, SEXP sI, |
| double (*f)(double, double, int), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, n, na, nb; |
| double ai, bi, *y; |
| const double *a, *b; |
| int m_opt; |
| int naflag; |
| |
| if (!isNumeric(sa) || !isNumeric(sb)) |
| errorcall(lcall, R_MSG_NONNUM_MATH); |
| |
| SETUP_Math2; |
| m_opt = asInteger(sI); |
| |
| MOD_ITERATE2(n, na, nb, i, ia, ib, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| if_NA_Math2_set(y[i], ai, bi) |
| else { |
| y[i] = f(ai, bi, m_opt); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| FINISH_Math2; |
| return sy; |
| } /* math2_1() */ |
| |
| static SEXP math2_2(SEXP sa, SEXP sb, SEXP sI1, SEXP sI2, |
| double (*f)(double, double, int, int), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, n, na, nb; |
| double ai, bi, *y; |
| const double *a, *b; |
| int i_1, i_2; |
| int naflag; |
| if (!isNumeric(sa) || !isNumeric(sb)) |
| errorcall(lcall, R_MSG_NONNUM_MATH); |
| |
| SETUP_Math2; |
| i_1 = asInteger(sI1); |
| i_2 = asInteger(sI2); |
| |
| MOD_ITERATE2(n, na, nb, i, ia, ib, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| if_NA_Math2_set(y[i], ai, bi) |
| else { |
| y[i] = f(ai, bi, i_1, i_2); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| FINISH_Math2; |
| return sy; |
| } /* math2_2() */ |
| |
| /* This is only used directly by .Internal for Bessel functions, |
| so managing R_alloc stack is only prudence */ |
| static SEXP math2B(SEXP sa, SEXP sb, double (*f)(double, double, double *), |
| SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, n, na, nb; |
| double ai, bi, *y; |
| const double *a, *b; |
| int naflag; |
| double amax, *work; |
| size_t nw; |
| |
| #define besselJY_max_nu 1e7 |
| |
| if (!isNumeric(sa) || !isNumeric(sb)) |
| errorcall(lcall, R_MSG_NONNUM_MATH); |
| |
| /* for 0-length a we want the attributes of a, not those of b |
| as no recycling will occur */ |
| SETUP_Math2; |
| |
| /* allocate work array for BesselJ, BesselY large enough for all |
| arguments */ |
| amax = 0.0; |
| for (i = 0; i < nb; i++) { |
| double av = b[i] < 0 ? -b[i] : b[i]; |
| if (amax < av) |
| amax = av; |
| } |
| if (amax > besselJY_max_nu) |
| amax = besselJY_max_nu; // and warning will happen in ../nmath/bessel_[jy].c |
| const void *vmax = vmaxget(); |
| nw = 1 + (size_t)floor(amax); |
| work = (double *) R_alloc(nw, sizeof(double)); |
| |
| MOD_ITERATE2(n, na, nb, i, ia, ib, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| if_NA_Math2_set(y[i], ai, bi) |
| else { |
| y[i] = f(ai, bi, work); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| |
| vmaxset(vmax); |
| FINISH_Math2; |
| |
| return sy; |
| } /* math2B() */ |
| |
| #define Math2(A, FUN) math2(CAR(A), CADR(A), FUN, call); |
| #define Math2_1(A, FUN) math2_1(CAR(A), CADR(A), CADDR(A), FUN, call); |
| #define Math2_2(A, FUN) math2_2(CAR(A), CADR(A), CADDR(A), CADDDR(A), FUN, call) |
| #define Math2B(A, FUN) math2B(CAR(A), CADR(A), FUN, call); |
| |
| SEXP attribute_hidden do_math2(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| /* For .Internals, fix up call so errorcall() behaves like error(). */ |
| if (TYPEOF(CAR(call)) == SYMSXP && INTERNAL(CAR(call)) == op) |
| call = R_CurrentExpression; |
| |
| checkArity(op, args); |
| |
| if (isComplex(CAR(args)) || |
| (PRIMVAL(op) == 0 && isComplex(CADR(args)))) |
| return complex_math2(call, op, args, env); |
| |
| |
| switch (PRIMVAL(op)) { |
| |
| case 0: return Math2(args, atan2); |
| case 10001: return Math2(args, fround);// round(), ../nmath/fround.c |
| case 10004: return Math2(args, fprec); // signif(), ../nmath/fprec.c |
| |
| case 2: return Math2(args, lbeta); |
| case 3: return Math2(args, beta); |
| case 4: return Math2(args, lchoose); |
| case 5: return Math2(args, choose); |
| |
| case 6: return Math2_1(args, dchisq); |
| case 7: return Math2_2(args, pchisq); |
| case 8: return Math2_2(args, qchisq); |
| |
| case 9: return Math2_1(args, dexp); |
| case 10: return Math2_2(args, pexp); |
| case 11: return Math2_2(args, qexp); |
| |
| case 12: return Math2_1(args, dgeom); |
| case 13: return Math2_2(args, pgeom); |
| case 14: return Math2_2(args, qgeom); |
| |
| case 15: return Math2_1(args, dpois); |
| case 16: return Math2_2(args, ppois); |
| case 17: return Math2_2(args, qpois); |
| |
| case 18: return Math2_1(args, dt); |
| case 19: return Math2_2(args, pt); |
| case 20: return Math2_2(args, qt); |
| |
| case 21: return Math2_1(args, dsignrank); |
| case 22: return Math2_2(args, psignrank); |
| case 23: return Math2_2(args, qsignrank); |
| |
| case 24: return Math2B(args, bessel_j_ex); |
| case 25: return Math2B(args, bessel_y_ex); |
| case 26: return Math2(args, psigamma); |
| |
| default: |
| error(_("unimplemented real function of %d numeric arguments"), 2); |
| } |
| return op; /* never used; to keep -Wall happy */ |
| } |
| |
| |
| /* The S4 Math2 group, round and signif */ |
| /* This is a primitive SPECIALSXP with internal argument matching */ |
| SEXP attribute_hidden do_Math2(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP res, call2; |
| int n, nprotect = 2; |
| static SEXP do_Math2_formals = NULL; |
| |
| if (length(args) >= 2 && |
| isSymbol(CADR(args)) && R_isMissing(CADR(args), env)) { |
| double digits = 0; |
| if(PRIMVAL(op) == 10004) digits = 6.0; // for signif() |
| PROTECT(args = list2(CAR(args), ScalarReal(digits))); nprotect++; |
| } |
| |
| PROTECT(args = evalListKeepMissing(args, env)); |
| PROTECT(call2 = lang2(CAR(call), R_NilValue)); |
| SETCDR(call2, args); |
| |
| n = length(args); |
| if (n != 1 && n != 2) |
| error(ngettext("%d argument passed to '%s' which requires 1 or 2 arguments", |
| "%d arguments passed to '%s'which requires 1 or 2 arguments", n), |
| n, PRIMNAME(op)); |
| |
| if (! DispatchGroup("Math", call2, op, args, env, &res)) { |
| if(n == 1) { |
| double digits = 0.0; |
| if(PRIMVAL(op) == 10004) digits = 6.0; |
| SETCDR(args, CONS(ScalarReal(digits), R_NilValue)); |
| } else { |
| /* If named, do argument matching by name */ |
| if (TAG(args) != R_NilValue || TAG(CDR(args)) != R_NilValue) { |
| if (do_Math2_formals == NULL) |
| do_Math2_formals = allocFormalsList2(install("x"), |
| install("digits")); |
| PROTECT(args = matchArgs(do_Math2_formals, args, call)); |
| nprotect++; |
| } |
| if (length(CADR(args)) == 0) |
| errorcall(call, _("invalid second argument of length 0")); |
| } |
| res = do_math2(call, op, args, env); |
| } |
| UNPROTECT(nprotect); |
| return res; |
| } |
| |
| /* log{2,10}() builtins : */ |
| SEXP attribute_hidden do_log1arg(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| SEXP res, call2, args2, tmp = R_NilValue /* -Wall */; |
| |
| checkArity(op, args); |
| check1arg(args, call, "x"); |
| |
| if (DispatchGroup("Math", call, op, args, env, &res)) return res; |
| |
| SEXP sLog = install("log"); |
| if(PRIMVAL(op) == 10) tmp = ScalarReal(10.0); |
| if(PRIMVAL(op) == 2) tmp = ScalarReal(2.0); |
| |
| PROTECT(call2 = lang3(sLog, CAR(args), tmp)); |
| PROTECT(args2 = lang2(CAR(args), tmp)); |
| if (! DispatchGroup("Math", call2, op, args2, env, &res)) { |
| if (isComplex(CAR(args))) |
| res = complex_math2(call2, op, args2, env); |
| else |
| res = math2(CAR(args), tmp, logbase, call); |
| } |
| UNPROTECT(2); |
| return res; |
| } |
| |
| #ifdef M_E |
| # define DFLT_LOG_BASE M_E |
| #else |
| # define DFLT_LOG_BASE exp(1.) |
| #endif |
| |
| /* do_log is a primitive SPECIALSXP with internal argument |
| matching. do_log_builtin is the BUILTIN version that expects |
| evaluated arguments to be passed as 'args', expect that these may |
| contain missing arguments. */ |
| SEXP attribute_hidden do_log(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| args = evalListKeepMissing(args, env); |
| return do_log_builtin(call, op, args, env); |
| } |
| |
| SEXP attribute_hidden do_log_builtin(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| PROTECT(args); |
| int n = length(args); |
| SEXP res; |
| |
| if (n == 1 && TAG(args) == R_NilValue) { |
| /* log(x) is handled here */ |
| SEXP x = CAR(args); |
| if (x != R_MissingArg && ! OBJECT(x)) { |
| if (isComplex(x)) |
| res = complex_math1(call, op, args, env); |
| else |
| res = math1(x, R_log, call); |
| UNPROTECT(1); |
| return res; |
| } |
| } |
| else if (n == 2 && |
| TAG(args) == R_NilValue && |
| (TAG(CDR(args)) == R_NilValue || TAG(CDR(args)) == R_BaseSymbol)) { |
| /* log(x, y) or log(x, base = y) are handled here */ |
| SEXP x = CAR(args); |
| SEXP y = CADR(args); |
| if (x != R_MissingArg && y != R_MissingArg && |
| ! OBJECT(x) && ! OBJECT(y)) { |
| if (isComplex(x) || isComplex(y)) |
| res = complex_math2(call, op, args, env); |
| else |
| res = math2(x, y, logbase, call); |
| UNPROTECT(1); |
| return res; |
| } |
| } |
| |
| static SEXP do_log_formals = NULL; |
| static SEXP R_x_Symbol = NULL; |
| if (do_log_formals == NULL) { |
| R_x_Symbol = install("x"); |
| do_log_formals = allocFormalsList2(R_x_Symbol, R_BaseSymbol); |
| } |
| |
| if (n == 1) { |
| if (CAR(args) == R_MissingArg || |
| (TAG(args) != R_NilValue && TAG(args) != R_x_Symbol)) |
| error(_("argument \"%s\" is missing, with no default"), "x"); |
| |
| if (! DispatchGroup("Math", call, op, args, env, &res)) { |
| if (isComplex(CAR(args))) |
| res = complex_math1(call, op, args, env); |
| else |
| res = math1(CAR(args), R_log, call); |
| } |
| UNPROTECT(1); |
| return res; |
| } |
| else { |
| /* match argument names if supplied */ |
| /* will signal an error unless there are one or two arguments */ |
| /* after the match, length(args) will be 2 */ |
| PROTECT(args = matchArgs(do_log_formals, args, call)); |
| |
| if(CAR(args) == R_MissingArg) |
| error(_("argument \"%s\" is missing, with no default"), "x"); |
| if (CADR(args) == R_MissingArg) |
| SETCADR(args, ScalarReal(DFLT_LOG_BASE)); |
| |
| if (! DispatchGroup("Math", call, op, args, env, &res)) { |
| if (length(CADR(args)) == 0) |
| errorcall(call, _("invalid argument 'base' of length 0")); |
| if (isComplex(CAR(args)) || isComplex(CADR(args))) |
| res = complex_math2(call, op, args, env); |
| else |
| res = math2(CAR(args), CADR(args), logbase, call); |
| } |
| UNPROTECT(2); |
| return res; |
| } |
| } |
| |
| |
| /* Mathematical Functions of Three (Real) Arguments */ |
| |
| /* math3_1 and math3_2 and related can be removed once the byte |
| compiler knows how to optimize to .External rather than |
| .Internal */ |
| |
| |
| #define if_NA_Math3_set(y,a,b,c) \ |
| if (ISNA (a) || ISNA (b)|| ISNA (c)) y = NA_REAL; \ |
| else if (ISNAN(a) || ISNAN(b)|| ISNAN(c)) y = R_NaN; |
| |
| #define SETUP_Math3 \ |
| if (!isNumeric(sa) || !isNumeric(sb) || !isNumeric(sc)) \ |
| error(R_MSG_NONNUM_MATH); \ |
| \ |
| na = XLENGTH(sa); \ |
| nb = XLENGTH(sb); \ |
| nc = XLENGTH(sc); \ |
| if ((na == 0) || (nb == 0) || (nc == 0)) \ |
| return(allocVector(REALSXP, 0)); \ |
| n = na; \ |
| if (n < nb) n = nb; \ |
| if (n < nc) n = nc; \ |
| PROTECT(sa = coerceVector(sa, REALSXP)); \ |
| PROTECT(sb = coerceVector(sb, REALSXP)); \ |
| PROTECT(sc = coerceVector(sc, REALSXP)); \ |
| PROTECT(sy = allocVector(REALSXP, n)); \ |
| a = REAL_RO(sa); \ |
| b = REAL_RO(sb); \ |
| c = REAL_RO(sc); \ |
| y = REAL(sy); \ |
| naflag = 0 |
| |
| #define FINISH_Math3 \ |
| if(naflag) warning(R_MSG_NA); \ |
| \ |
| if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \ |
| else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \ |
| else if (n == nc) SHALLOW_DUPLICATE_ATTRIB(sy, sc); \ |
| UNPROTECT(4) |
| |
| static SEXP math3_1(SEXP sa, SEXP sb, SEXP sc, SEXP sI, |
| double (*f)(double, double, double, int), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, ic, n, na, nb, nc; |
| double ai, bi, ci, *y; |
| const double *a, *b, *c; |
| int i_1; |
| int naflag; |
| |
| SETUP_Math3; |
| i_1 = asInteger(sI); |
| |
| MOD_ITERATE3(n, na, nb, nc, i, ia, ib, ic, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| ci = c[ic]; |
| if_NA_Math3_set(y[i], ai,bi,ci) |
| else { |
| y[i] = f(ai, bi, ci, i_1); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| |
| FINISH_Math3; |
| return sy; |
| } /* math3_1 */ |
| |
| static SEXP math3_2(SEXP sa, SEXP sb, SEXP sc, SEXP sI, SEXP sJ, |
| double (*f)(double, double, double, int, int), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, ic, n, na, nb, nc; |
| double ai, bi, ci, *y; |
| const double *a, *b, *c; |
| int i_1,i_2; |
| int naflag; |
| |
| SETUP_Math3; |
| i_1 = asInteger(sI); |
| i_2 = asInteger(sJ); |
| |
| MOD_ITERATE3 (n, na, nb, nc, i, ia, ib, ic, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| ci = c[ic]; |
| if_NA_Math3_set(y[i], ai,bi,ci) |
| else { |
| y[i] = f(ai, bi, ci, i_1, i_2); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| |
| FINISH_Math3; |
| return sy; |
| } /* math3_2 */ |
| |
| /* This is only used directly by .Internal for Bessel functions, |
| so managing R_alloc stack is only prudence */ |
| static SEXP math3B(SEXP sa, SEXP sb, SEXP sc, |
| double (*f)(double, double, double, double *), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, ic, n, na, nb, nc; |
| double ai, bi, ci, *y; |
| const double *a, *b, *c; |
| int naflag; |
| double amax, *work; |
| size_t nw; |
| |
| SETUP_Math3; |
| |
| /* allocate work array for BesselI, BesselK large enough for all |
| arguments */ |
| amax = 0.0; |
| for (i = 0; i < nb; i++) { |
| double av = b[i] < 0 ? -b[i] : b[i]; |
| if (av > amax) amax = av; |
| } |
| const void *vmax = vmaxget(); |
| nw = 1 + (size_t)floor(amax); |
| work = (double *) R_alloc(nw, sizeof(double)); |
| |
| MOD_ITERATE3 (n, na, nb, nc, i, ia, ib, ic, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| ci = c[ic]; |
| if_NA_Math3_set(y[i], ai,bi,ci) |
| else { |
| y[i] = f(ai, bi, ci, work); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| |
| FINISH_Math3; |
| vmaxset(vmax); |
| |
| return sy; |
| } /* math3B */ |
| |
| #define Math3_1(A, FUN) math3_1(CAR(A), CADR(A), CADDR(A), CADDDR(A), FUN, call); |
| #define Math3_2(A, FUN) math3_2(CAR(A), CADR(A), CADDR(A), CADDDR(A), CAD4R(A), FUN, call) |
| #define Math3B(A, FUN) math3B (CAR(A), CADR(A), CADDR(A), FUN, call); |
| |
| SEXP attribute_hidden do_math3(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| checkArity(op, args); |
| |
| switch (PRIMVAL(op)) { |
| |
| case 1: return Math3_1(args, dbeta); |
| case 2: return Math3_2(args, pbeta); |
| case 3: return Math3_2(args, qbeta); |
| |
| case 4: return Math3_1(args, dbinom); |
| case 5: return Math3_2(args, pbinom); |
| case 6: return Math3_2(args, qbinom); |
| |
| case 7: return Math3_1(args, dcauchy); |
| case 8: return Math3_2(args, pcauchy); |
| case 9: return Math3_2(args, qcauchy); |
| |
| case 10: return Math3_1(args, df); |
| case 11: return Math3_2(args, pf); |
| case 12: return Math3_2(args, qf); |
| |
| case 13: return Math3_1(args, dgamma); |
| case 14: return Math3_2(args, pgamma); |
| case 15: return Math3_2(args, qgamma); |
| |
| case 16: return Math3_1(args, dlnorm); |
| case 17: return Math3_2(args, plnorm); |
| case 18: return Math3_2(args, qlnorm); |
| |
| case 19: return Math3_1(args, dlogis); |
| case 20: return Math3_2(args, plogis); |
| case 21: return Math3_2(args, qlogis); |
| |
| case 22: return Math3_1(args, dnbinom); |
| case 23: return Math3_2(args, pnbinom); |
| case 24: return Math3_2(args, qnbinom); |
| |
| case 25: return Math3_1(args, dnorm); |
| case 26: return Math3_2(args, pnorm); |
| case 27: return Math3_2(args, qnorm); |
| |
| case 28: return Math3_1(args, dunif); |
| case 29: return Math3_2(args, punif); |
| case 30: return Math3_2(args, qunif); |
| |
| case 31: return Math3_1(args, dweibull); |
| case 32: return Math3_2(args, pweibull); |
| case 33: return Math3_2(args, qweibull); |
| |
| case 34: return Math3_1(args, dnchisq); |
| case 35: return Math3_2(args, pnchisq); |
| case 36: return Math3_2(args, qnchisq); |
| |
| case 37: return Math3_1(args, dnt); |
| case 38: return Math3_2(args, pnt); |
| case 39: return Math3_2(args, qnt); |
| |
| case 40: return Math3_1(args, dwilcox); |
| case 41: return Math3_2(args, pwilcox); |
| case 42: return Math3_2(args, qwilcox); |
| |
| case 43: return Math3B(args, bessel_i_ex); |
| case 44: return Math3B(args, bessel_k_ex); |
| |
| case 45: return Math3_1(args, dnbinom_mu); |
| case 46: return Math3_2(args, pnbinom_mu); |
| case 47: return Math3_2(args, qnbinom_mu); |
| |
| default: |
| error(_("unimplemented real function of %d numeric arguments"), 3); |
| } |
| return op; /* never used; to keep -Wall happy */ |
| } /* do_math3() */ |
| |
| /* Mathematical Functions of Four (Real) Arguments */ |
| |
| /* This can be removed completely once the byte compiler knows how to |
| optimize to .External rather than .Internal */ |
| |
| #define if_NA_Math4_set(y,a,b,c,d) \ |
| if (ISNA (a)|| ISNA (b)|| ISNA (c)|| ISNA (d)) y = NA_REAL;\ |
| else if (ISNAN(a)|| ISNAN(b)|| ISNAN(c)|| ISNAN(d)) y = R_NaN; |
| |
| static SEXP math4(SEXP sa, SEXP sb, SEXP sc, SEXP sd, |
| double (*f)(double, double, double, double), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, ic, id, n, na, nb, nc, nd; |
| double ai, bi, ci, di, *y; |
| const double *a, *b, *c, *d; |
| int naflag; |
| |
| #define SETUP_Math4 \ |
| if(!isNumeric(sa)|| !isNumeric(sb)|| !isNumeric(sc)|| !isNumeric(sd))\ |
| error(R_MSG_NONNUM_MATH); \ |
| \ |
| na = XLENGTH(sa); \ |
| nb = XLENGTH(sb); \ |
| nc = XLENGTH(sc); \ |
| nd = XLENGTH(sd); \ |
| if ((na == 0) || (nb == 0) || (nc == 0) || (nd == 0)) \ |
| return(allocVector(REALSXP, 0)); \ |
| n = na; \ |
| if (n < nb) n = nb; \ |
| if (n < nc) n = nc; \ |
| if (n < nd) n = nd; \ |
| PROTECT(sa = coerceVector(sa, REALSXP)); \ |
| PROTECT(sb = coerceVector(sb, REALSXP)); \ |
| PROTECT(sc = coerceVector(sc, REALSXP)); \ |
| PROTECT(sd = coerceVector(sd, REALSXP)); \ |
| PROTECT(sy = allocVector(REALSXP, n)); \ |
| a = REAL_RO(sa); \ |
| b = REAL_RO(sb); \ |
| c = REAL_RO(sc); \ |
| d = REAL_RO(sd); \ |
| y = REAL(sy); \ |
| naflag = 0 |
| |
| SETUP_Math4; |
| |
| MOD_ITERATE4 (n, na, nb, nc, nd, i, ia, ib, ic, id, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| ci = c[ic]; |
| di = d[id]; |
| if_NA_Math4_set(y[i], ai,bi,ci,di) |
| else { |
| y[i] = f(ai, bi, ci, di); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| |
| #define FINISH_Math4 \ |
| if(naflag) warning(R_MSG_NA); \ |
| \ |
| if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \ |
| else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \ |
| else if (n == nc) SHALLOW_DUPLICATE_ATTRIB(sy, sc); \ |
| else if (n == nd) SHALLOW_DUPLICATE_ATTRIB(sy, sd); \ |
| UNPROTECT(5) |
| |
| FINISH_Math4; |
| |
| return sy; |
| } /* math4() */ |
| |
| static SEXP math4_1(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP sI, double (*f)(double, double, double, double, int), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, ic, id, n, na, nb, nc, nd; |
| double ai, bi, ci, di, *y; |
| const double *a, *b, *c, *d; |
| int i_1; |
| int naflag; |
| |
| SETUP_Math4; |
| i_1 = asInteger(sI); |
| |
| MOD_ITERATE4 (n, na, nb, nc, nd, i, ia, ib, ic, id, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| ci = c[ic]; |
| di = d[id]; |
| if_NA_Math4_set(y[i], ai,bi,ci,di) |
| else { |
| y[i] = f(ai, bi, ci, di, i_1); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| FINISH_Math4; |
| return sy; |
| } /* math4_1() */ |
| |
| static SEXP math4_2(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP sI, SEXP sJ, |
| double (*f)(double, double, double, double, int, int), SEXP lcall) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, ic, id, n, na, nb, nc, nd; |
| double ai, bi, ci, di, *y; |
| const double *a, *b, *c, *d; |
| int i_1, i_2; |
| int naflag; |
| |
| SETUP_Math4; |
| i_1 = asInteger(sI); |
| i_2 = asInteger(sJ); |
| |
| MOD_ITERATE4 (n, na, nb, nc, nd, i, ia, ib, ic, id, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| ci = c[ic]; |
| di = d[id]; |
| if_NA_Math4_set(y[i], ai,bi,ci,di) |
| else { |
| y[i] = f(ai, bi, ci, di, i_1, i_2); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| FINISH_Math4; |
| return sy; |
| } /* math4_2() */ |
| |
| |
| #define CAD3R CADDDR |
| /* This is not (yet) in Rinternals.h : */ |
| #define CAD5R(e) CAR(CDR(CDR(CDR(CDR(CDR(e)))))) |
| |
| #define Math4(A, FUN) math4 (CAR(A), CADR(A), CADDR(A), CAD3R(A), FUN, call) |
| #define Math4_1(A, FUN) math4_1(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), \ |
| FUN, call) |
| #define Math4_2(A, FUN) math4_2(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), \ |
| CAD5R(A), FUN, call) |
| |
| |
| SEXP attribute_hidden do_math4(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| checkArity(op, args); |
| |
| |
| switch (PRIMVAL(op)) { |
| |
| /* Completely dummy for -Wall -- math4() at all! : */ |
| case -99: return Math4(args, (double (*)(double, double, double, double))NULL); |
| |
| case 1: return Math4_1(args, dhyper); |
| case 2: return Math4_2(args, phyper); |
| case 3: return Math4_2(args, qhyper); |
| |
| case 4: return Math4_1(args, dnbeta); |
| case 5: return Math4_2(args, pnbeta); |
| case 6: return Math4_2(args, qnbeta); |
| case 7: return Math4_1(args, dnf); |
| case 8: return Math4_2(args, pnf); |
| case 9: return Math4_2(args, qnf); |
| #ifdef UNIMP |
| case 10: return Math4_1(args, dtukey); |
| #endif |
| case 11: return Math4_2(args, ptukey); |
| case 12: return Math4_2(args, qtukey); |
| default: |
| error(_("unimplemented real function of %d numeric arguments"), 4); |
| } |
| return op; /* never used; to keep -Wall happy */ |
| } |
| |
| |
| #ifdef WHEN_MATH5_IS_THERE/* as in ./arithmetic.h */ |
| |
| /* Mathematical Functions of Five (Real) Arguments */ |
| |
| #define if_NA_Math5_set(y,a,b,c,d,e) \ |
| if (ISNA (a)|| ISNA (b)|| ISNA (c)|| ISNA (d)|| ISNA (e)) \ |
| y = NA_REAL; \ |
| else if(ISNAN(a)|| ISNAN(b)|| ISNAN(c)|| ISNAN(d)|| ISNAN(e)) \ |
| y = R_NaN; |
| |
| static SEXP math5(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP se, double (*f)()) |
| { |
| SEXP sy; |
| R_xlen_t i, ia, ib, ic, id, ie, n, na, nb, nc, nd, ne; |
| double ai, bi, ci, di, ei, *y; |
| const double *a, *b, *c, *d, *e; |
| |
| #define SETUP_Math5 \ |
| if (!isNumeric(sa) || !isNumeric(sb) || !isNumeric(sc) || \ |
| !isNumeric(sd) || !isNumeric(se)) \ |
| error(R_MSG_NONNUM_MATH); \ |
| \ |
| na = XLENGTH(sa); \ |
| nb = XLENGTH(sb); \ |
| nc = XLENGTH(sc); \ |
| nd = XLENGTH(sd); \ |
| ne = XLENGTH(se); \ |
| if ((na == 0) || (nb == 0) || (nc == 0) || (nd == 0) || (ne == 0)) \ |
| return(allocVector(REALSXP, 0)); \ |
| n = na; \ |
| if (n < nb) n = nb; \ |
| if (n < nc) n = nc; \ |
| if (n < nd) n = nd; \ |
| if (n < ne) n = ne; /* n = max(na,nb,nc,nd,ne) */ \ |
| PROTECT(sa = coerceVector(sa, REALSXP)); \ |
| PROTECT(sb = coerceVector(sb, REALSXP)); \ |
| PROTECT(sc = coerceVector(sc, REALSXP)); \ |
| PROTECT(sd = coerceVector(sd, REALSXP)); \ |
| PROTECT(se = coerceVector(se, REALSXP)); \ |
| PROTECT(sy = allocVector(REALSXP, n)); \ |
| a = REAL_RO(sa); \ |
| b = REAL_RO(sb); \ |
| c = REAL_RO(sc); \ |
| d = REAL_RO(sd); \ |
| e = REAL_RO(se); \ |
| y = REAL(sy); \ |
| naflag = 0 |
| |
| SETUP_Math5; |
| |
| MOD_ITERATE5 (n, na, nb, nc, nd, ne, |
| i, ia, ib, ic, id, ie, { |
| // if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt(); |
| ai = a[ia]; |
| bi = b[ib]; |
| ci = c[ic]; |
| di = d[id]; |
| ei = e[ie]; |
| if_NA_Math5_set(y[i], ai,bi,ci,di,ei) |
| else { |
| y[i] = f(ai, bi, ci, di, ei); |
| if (ISNAN(y[i])) naflag = 1; |
| } |
| }); |
| |
| #define FINISH_Math5 \ |
| if(naflag) warning(R_MSG_NA); \ |
| \ |
| if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \ |
| else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \ |
| else if (n == nc) SHALLOW_DUPLICATE_ATTRIB(sy, sc); \ |
| else if (n == nd) SHALLOW_DUPLICATE_ATTRIB(sy, sd); \ |
| else if (n == ne) SHALLOW_DUPLICATE_ATTRIB(sy, se); \ |
| UNPROTECT(6) |
| |
| FINISH_Math5; |
| |
| return sy; |
| } /* math5() */ |
| |
| #define Math5(A, FUN) \ |
| math5(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), FUN); |
| |
| SEXP attribute_hidden do_math5(SEXP call, SEXP op, SEXP args, SEXP env) |
| { |
| checkArity(op, args); |
| lcall = call; |
| |
| switch (PRIMVAL(op)) { |
| |
| /* Completely dummy for -Wall -- use math5() at all! : */ |
| case -99: return Math5(args, dhyper); |
| #ifdef UNIMP |
| case 2: return Math5(args, p...); |
| case 3: return Math5(args, q...); |
| #endif |
| default: |
| error(_("unimplemented real function of %d numeric arguments"), 5); |
| } |
| return op; /* never used; to keep -Wall happy */ |
| } /* do_math5() */ |
| |
| #endif /* Math5 is there */ |
| |
| /* This is used for experimenting with parallelized nmath functions -- LT */ |
| CCODE R_get_arith_function(int which) |
| { |
| switch (which) { |
| case 1: return do_math1; |
| case 2: return do_math2; |
| case 3: return do_math3; |
| case 4: return do_math4; |
| case 11: return complex_math1; |
| case 12: return complex_math2; |
| default: error("bad arith function index"); return NULL; |
| } |
| } |