The R Project SVN R

    if (x >= -11. && x <= 11.) {

	if(y == 4.)

	    return x * x * x * x;

	if(y == 3.)

	    return x * x * x;

    }

/** compute f(sa)  in "ari"thmetic - with "extras", notably ensuring

 *  <correct f>(arg) = res *and* correct treatment of NAs in sa */

static SEXP math1_ari(SEXP sa, double(*f)(double), double arg, double res, 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 (when sy = sa) */

	if (x == arg)

	    y[i] = res;

	else

	/* 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;

}

/* Make the result precise for squares of integers up to 11 */

// needed?? I see exact sqrt() for all of {0:1e8} with glibc R (Linux Fedora 42):

//     chkRt <- function(n) {zN <- 0:n; stopifnot(sqrt(zN^2) == zN) } ; chkRt(1e8)

static double Rsqrt(double x)

{

    if (x == 0.) return x; // => sqrt(-0.) = -0.

    for(int i=1; i < 12; i++) {

	if (x == i*i) return i;

    }

    return sqrt(x);

}

static double Rexp(double x)

{

    /* exp(x) = 1 + x + x^2/2! + ...

     * should return 1+x for very small x i.e.,  x^2/2 < D_EPS / 2

     * <==> x^2 < D_EPS  <==> |x| < sqrt(D_EPS) : */

    return (fabs(x) <= sqrt(DBL_EPSILON)) ? 1. + x : exp(x);

}

/* should return x for very small x (e.g. expm1(x))*/

static double f_x_x(double x, double(*f)(double), double m)

{

    return (fabs(x) <= m) ? x : f(x);

}

static double Rexpm1(double x)

{

    /* expm1(x) = exp(x) - 1 =  x + x^2/2 + O(x^3)

     *                       =. x  when  x^2/2 < |x| * D_EPS/2

     *                             <==>  |x|   <  D_EPS */

    return f_x_x(x, expm1, DBL_EPSILON);

}

static double Rlog1p(double x)

{

    /* log1p(x) = log(1 + x) =  x - x^2/2  + O(x^3)

     *                       =. x  when  |x| < D_EPS, see Rexpm1() */

    return f_x_x(x, log1p, DBL_EPSILON);

}

static double Rsin(double x)

{

    /* sin(x) = x - x^3/6 + O(x^5) = x*(1 - x^2/6) + O(.) =!= x  iff

       (1 - x^2/6) .= 1  <==>  |x| < sqrt(6 eps); eps = DBL_Eps/2 */

    return f_x_x(x, sin, sqrt(3. * DBL_EPSILON));

}

static double Rtan(double x)

{

    /* tan(x) =  x + x^3/3 + O(x^5)

              =. x  when |x|^3/3 < |x| EPS/2  <==>

	                  x^2    < EPS * 3/2  <==>

 	                   |x|   < sqrt(3/2 * EPS) */

    return f_x_x(x, tan, sqrt(1.5 * DBL_EPSILON));

}

static double Rcos(double x)

{

    /* cos(x) =  1 - x^2/2! + x^4/4!

              =. 1 - x^2/2!  iff  x^4/24 < (1 - x^2/2) * EPS/2  ~= EPS/2 <==>

 	                          x^4 < 12*EPS  */

    if (fabs(x) < sqrt(sqrt(12. * DBL_EPSILON)))

	return (1. - x*x*0.5);

    else

	return cos(x);

}

static double Rasin(double x)

{

    /* asin(x) =  x + x^3/6 + 3/40 * x^5...   ==>

               =. x when |x| < sqrt(3 * EPS)  -- see Rsin() above */

    return f_x_x(x, asin, sqrt(3. * DBL_EPSILON));

}

static double Ratan(double x)

{

    /* atan(x) =  x - x^3/3 + x^5/5 ... ==>

               =. x  iff |x| < sqrt(3/2 * EPS) -- see Rtan() above*/

    return f_x_x(x, atan, sqrt(1.5 * DBL_EPSILON));

}

#define MATH1_SQRT(x) math1(CAR(args), Rsqrt, call);

#define MATH1_EXP(x)  math1_ari(CAR(args), Rexp,  0., 1., call);

#define MATH1_EXPM1(x)math1_ari(CAR(args), Rexpm1,0., 0., call);

#define MATH1_LOG1P(x)math1_ari(CAR(args), Rlog1p,0., 0., call);

#define MATH1_SIN(x)  math1_ari(CAR(args), Rsin,  0., 0., call);

#define MATH1_ASIN(x) math1_ari(CAR(args), Rasin, 0., 0., call);

#define MATH1_TAN(x)  math1_ari(CAR(args), Rtan,  0., 0., call);

#define MATH1_ATAN(x) math1_ari(CAR(args), Ratan, 0., 0., call);

#define MATH1_COS(x)  math1_ari(CAR(args), Rcos,  0., 1., call);

#define MATH1_00(x)   math1_ari(CAR(args),  x,    0., 0., call);

#define MATH1_01(x)   math1_ari(CAR(args),  x,    0., 1., call);

#define MATH1_10(x)   math1_ari(CAR(args),  x,    1., 0., call);