The R Project SVN – R – Diff – /branches/R-exp-API/src/main/complex.c

 #include <Rconfig.h>
 #endif
 #include "Defn.h"
 #include "Mathlib.h"
-#include "Fortran.h"/* POW_DI */
+#include "Fortran.h"		/* POW_DI */
-#include "Applic.h"/* cpoly */
+#include "Applic.h"		/* cpoly */
-#include "arithmetic.h"/* complex_ */
+#include "arithmetic.h"		/* complex_ */
 static int naflag;
 SEXP complex_unary(int code, SEXP s1)
+{
+	    }
 #endif
+	}
 	return ans;
     default:
-	error("illegal complex unary operator");
+	error("illegal complex unary operator\n");
-	return R_NilValue;/* -Wall*/
+	return R_NilValue;	/* -Wall*/
+    }
+}
 static void complex_div(complex *c, complex *a, complex *b)
+{
 static void complex_pow(complex *r, complex *a, complex *b)
+{
     double logr, logi, x, y;
     int ib;
-    if(b->i == 0.) {/* be fast (and more accurate)*/
+    if(b->i == 0.) {		/* be fast (and more accurate)*/
+	if(b->r == 1.) {	/* a^1 */
-	if(b->r == 1.) { /* a^1 */ r->r = a->r; r->i = a->i; return;}
+	    r->r = a->r; r->i = a->i; return;
+	}
+	if(a->i == 0.) {
-	if(a->i == 0.) { r->r = pow(a->r, b->r); r->i = 0.; return;}
+	    r->r = pow(a->r, b->r); r->i = 0.; return;
+	}
 	if(a->r == 0. && b->r == (ib = (int)b->r)) {/* (|a|*i)^b */
 	    x = POW_DI(&(a->i), &ib);
-	    if(ib % 2) { /* ib is odd ==> imaginary r */
+	    if(ib % 2) {	/* ib is odd ==> imaginary r */
 		r->r = 0.;
 		r->i = ((ib>0 && ib %4 == 3)||(ib<0 && (-ib)%4 == 1))? -x : x;
-	    } else { /* even exponent b : real r */
+	    } else {		/* even exponent b : real r */
 		r->r = (ib %4)? -x : x; r->i = 0.;
+	    }
 	    return;
+	}
+    }
 /* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */
 SEXP do_cmathfuns(SEXP call, SEXP op, SEXP args, SEXP env)
+{
-    SEXP x, y = R_NilValue;/* -Wall*/
+    SEXP x, y = R_NilValue;	/* -Wall*/
     int i, n;
     checkArity(op, args);
     x = CAR(args);
     n = length(x);
 	    y = allocVector(REALSXP, n);
 	    for(i=0 ; i<n ; i++)
 		REAL(y)[i] = COMPLEX(x)[i].i;
 	    break;
 	case 3:	/* Mod */
+	case 6:	/* abs */
 	    y = allocVector(REALSXP, n);
 	    for(i=0 ; i<n ; i++) {
 #ifdef IEEE_754
 		REAL(y)[i] = hypot(COMPLEX(x)[i].r, COMPLEX(x)[i].i);
 #else
 		    REAL(y)[i] = REAL(x)[i];
 		else
 		    REAL(y)[i] = 0;
 	    break;
 	case 3:	/* Mod */
+	case 6:	/* abs */
 	    y = allocVector(REALSXP, n);
 	    for(i=0 ; i<n ; i++) {
 #ifdef IEEE_754
 		REAL(y)[i] = fabs(REAL(x)[i]);
 #else
+{
     r->r = rround(x->r, p->r);
     r->i = rround(x->i, p->r);
+}
-	/* Question:  This treats real and imaginary parts */
+/* Question:  This treats real and imaginary parts separately.  Should
-	/* separately.	Should it do them jointly? */
+   it do them jointly? */
 static void z_prec(complex *r, complex *x, complex *p)
+{
     r->r = prec(x->r, p->r);
     r->i = prec(x->i, p->r);
 SEXP complex_math1(SEXP call, SEXP op, SEXP args, SEXP env)
+{
     SEXP x, y;
     int n;
-    x = CAR(args);
+    PROTECT(x = CAR(args));
     n = length(x);
-    y = allocVector(CPLXSXP, n);
+    PROTECT(y = allocVector(CPLXSXP, n));
     naflag = 0;
     switch (PRIMVAL(op)) {
     case 10002: cmath1(z_atan, COMPLEX(x), COMPLEX(y), n); break;
     case 10003: cmath1(z_log, COMPLEX(x), COMPLEX(y), n); break;
-    case 0:
-	errorcall(call,
-		  "abs() unimplemented for complex; use Mod()");
-	break;
     case 3:  cmath1(z_sqrt, COMPLEX(x), COMPLEX(y), n); break;
     case 10: cmath1(z_exp, COMPLEX(x), COMPLEX(y), n); break;
     case 20: cmath1(z_cos, COMPLEX(x), COMPLEX(y), n); break;
 #endif
     default:
 	errorcall(call, "unimplemented complex function");
+    }
+    if (naflag)
-    if (naflag) warning("NAs produced in function \"%s\"", PRIMNAME(op));
+	warning("NAs produced in function \"%s\"", PRIMNAME(op));
+    ATTRIB(y) = duplicate(ATTRIB(x));
+    OBJECT(y) = OBJECT(x);
+    UNPROTECT(2);
     return y;
+}
 /* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */
     complex ai, bi, *a, *b, *y;
     SEXP sy;
     na = length(sa);
     nb = length(sb);
+    if ((na == 0) || (nb == 0))
+	return(allocVector(CPLXSXP, 0));
     n = (na < nb) ? nb : na;
     PROTECT(sa = coerceVector(sa, CPLXSXP));
     PROTECT(sb = coerceVector(sb, CPLXSXP));
     PROTECT(sy = allocVector(CPLXSXP, n));
     a = COMPLEX(sa);
     b = COMPLEX(sb);
     y = COMPLEX(sy);
-    if (na < 1 || na < 1) {
+    naflag = 0;
-	for (i = 0; i < n; i++) {
+    for (i = 0; i < n; i++) {
+	ai = a[i % na];
+	bi = b[i % nb];
+	if(ISNA(ai.r) && ISNA(ai.i) &&
+	   ISNA(bi.r) && ISNA(bi.i)) {
 	    y[i].r = NA_REAL;
 	    y[i].i = NA_REAL;
+	}
-    }
-    else {
+	else {
-	naflag = 0;
-	for (i = 0; i < n; i++) {
-	    ai = a[i % na];
+	    f(&y[i], &ai, &bi);
-	    bi = b[i % nb];
+#ifndef IEEE_754
-	    if(ISNA(ai.r) && ISNA(ai.i) &&
-	       ISNA(bi.r) && ISNA(bi.i)) {
+	    if(ISNA(y[i].r) || ISNA(y[i].i)) {
 		y[i].r = NA_REAL;
 		y[i].i = NA_REAL;
+		naflag = 1;
+	    }
-	    else {
-		f(&y[i], &ai, &bi);
-#ifndef IEEE_754
-		if(ISNA(y[i].r) || ISNA(y[i].i)) {
-		    y[i].r = NA_REAL;
-		    y[i].i = NA_REAL;
-		    naflag = 1;
-		}
 #endif
-	    }
+	}
+    }
+    if (naflag)
-    if (naflag) warning("NAs produced in function \"%s\"", PRIMNAME(op));
+	warning("NAs produced in function \"%s\"", PRIMNAME(op));
     if(n == na) {
 	ATTRIB(sy) = duplicate(ATTRIB(sa));
 	OBJECT(sy) = OBJECT(sa);
+    }
     else if(n == nb) {
     case 10004:
 	return cmath2(op, CAR(args), CADR(args), z_prec);
     case 0:
 	return cmath2(op, CAR(args), CADR(args), z_atan2);
     default:
-	errorcall(call, "unimplemented complex function");
+	errorcall(call, "unimplemented complex function\n");
-	return call;/* just for -Wall */
+	return call;		/* just for -Wall */
+    }
+}
 SEXP do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)
+{
 	PROTECT(zr = allocVector(REALSXP, n));
 	PROTECT(zi = allocVector(REALSXP, n));
 	for(i=0 ; i<n ; i++) {
 	    if(!R_FINITE(COMPLEX(z)[i].r) || !R_FINITE(COMPLEX(z)[i].i))
-		errorcall(call, "invalid polynomial coefficient");
+		errorcall(call, "invalid polynomial coefficient\n");
 	    REAL(zr)[degree-i] = COMPLEX(z)[i].r;
 	    REAL(zi)[degree-i] = COMPLEX(z)[i].i;
+	}
 	F77_SYMBOL(cpoly)(REAL(zr), REAL(zi), &degree,
 			  REAL(rr), REAL(ri), &fail);

The R Project SVN R

(root)/branches/R-exp-API/src/main/complex.c @ 28362 – Rev 5731 → 6098