The R Project SVN – R – Diff – /trunk/src/main/complex.c

 #include <Defn.h>		/* -> ../include/R_ext/Complex.h */
 #include <Internal.h>
 #include <Rmath.h>
 #include "arithmetic.h"		/* complex_*  */
-#include <complex.h>
+#include <complex.h>            /* incl "complex math.h" */
 #include "Rcomplex.h"		/* I, SET_C99_COMPLEX, toC99 */
 #include <R_ext/Itermacros.h>
 /* interval at which to check interrupts, a guess */
 #ifndef HAVE_CASIN
 #define casin R_casin
 static double complex casin(double complex z)
+{
     /* A&S 4.4.37 */
-    double alpha, t1, t2, x = creal(z), y = cimag(z), ri;
+    double x = creal(z), y = cimag(z),
-    t1 = 0.5 * hypot(x + 1, y);
+	t1 = 0.5 * hypot(x + 1, y),
-    t2 = 0.5 * hypot(x - 1, y);
+	t2 = 0.5 * hypot(x - 1, y),
-    alpha = t1 + t2;
+	alpha = t1 + t2,
-    ri = log(alpha + sqrt(alpha*alpha - 1));
+	ri = log(alpha + sqrt(alpha*alpha - 1));
     /* This comes from
        'z_asin() is continuous from below if x >= 1
 	and continuous from above if x <= -1.'
     */
     if(y < 0 || (y == 0 && x > 1)) ri *= -1;
 static double complex z_atanh(double complex z)
+{
     return -I * z_atan(z * I);
+}
+#ifdef HAVE_CABS
+# define R_CABS(Z) cabs(Z)
+#else
+# define R_CABS(Z) hypot(creal(Z), cimag(Z))
+#endif
 static bool cmath1(double complex (*f)(double complex),
 		       const Rcomplex *x, Rcomplex *y, R_xlen_t n)
+{
     R_xlen_t i;
     bool naflag = false;
 	/* Complex Functions of Two Arguments */
 typedef void (*cm2_fun)(Rcomplex *, Rcomplex *, Rcomplex *);
 attribute_hidden SEXP complex_math2(SEXP call, SEXP op, SEXP args, SEXP env)
+{
-    R_xlen_t i, n, na, nb, ia, ib;
-    Rcomplex ai, bi, *y;
-    const Rcomplex *a, *b;
-    SEXP sa, sb, sy;
-    bool naflag = false;
     cm2_fun f;
     switch (PRIMVAL(op)) {
     case 0: /* atan2 */
 	f = z_atan2; break;
     case 10001: /* round */
 	f = z_rround; break;
     case 10002: /* passed from do_log1arg */
-    case 10010:
+    case 10010: /*   "     "    " */
     case 10003: /* passed from do_log */
 	f = z_logbase; break;
     case 10004: /* signif */
 	f = z_prec; break;
     default:
 	error_return(_("unimplemented complex function"));
+    }
+    SEXP
-    PROTECT(sa = coerceVector(CAR(args), CPLXSXP));
+	sa = PROTECT(coerceVector(CAR(args), CPLXSXP)),
-    PROTECT(sb = coerceVector(CADR(args), CPLXSXP));
+	sb = PROTECT(coerceVector(CADR(args), CPLXSXP));
-    na = XLENGTH(sa); nb = XLENGTH(sb);
+    R_xlen_t na = XLENGTH(sa), nb = XLENGTH(sb), n;
     if ((na == 0) || (nb == 0)) {
 	UNPROTECT(2);
 	return(allocVector(CPLXSXP, 0));
+    }
     n = (na < nb) ? nb : na;
-    PROTECT(sy = allocVector(CPLXSXP, n));
+    SEXP sy = PROTECT(allocVector(CPLXSXP, n));
+    const Rcomplex
-    a = COMPLEX_RO(sa); b = COMPLEX_RO(sb);
+	*a = COMPLEX_RO(sa),
-    y = COMPLEX(sy);
+	*b = COMPLEX_RO(sb);
+    Rcomplex ai, bi, *y = COMPLEX(sy);
+    R_xlen_t i, ia, ib;
+    bool naflag = false;
     MOD_ITERATE2(n, na, nb, i, ia, ib, {
 	ai = a[ia]; bi = b[ib];
 	if(ISNA(ai.r) && ISNA(ai.i) &&
 	   ISNA(bi.r) && ISNA(bi.i)) {
 	    y[i].r = NA_REAL; y[i].i = NA_REAL;
 static void R_cpolyroot(double *opr, double *opi, int *degree,
 			double *zeror, double *zeroi, bool *fail);
 attribute_hidden SEXP do_polyroot(SEXP call, SEXP op, SEXP args, SEXP rho)
+{
-    SEXP z, zr, zi, r, rr, ri;
-    bool fail;
-    int degree, i, n;
     checkArity(op, args);
-    z = CAR(args);
+    SEXP z = CAR(args);
     switch(TYPEOF(z)) {
     case CPLXSXP:
 	PROTECT(z);
 	break;
     case REALSXP:
 	PROTECT(z = coerceVector(z, CPLXSXP));
 	break;
     default:
 	UNIMPLEMENTED_TYPE("polyroot", z);
+    }
+    int i, n;
 #ifdef LONG_VECTOR_SUPPORT
     R_xlen_t nn = XLENGTH(z);
     if (nn > R_SHORT_LEN_MAX) error("long vectors are not supported");
     n = (int) nn;
 #else
     n = LENGTH(z);
 #endif
     const Rcomplex *pz = COMPLEX_RO(z);
-    degree = 0;
+    int degree = 0; // := max{i; z[i] != 0}
     for(i = 0; i < n; i++) {
 	if(pz[i].r!= 0.0 || pz[i].i != 0.0) degree = i;
+    }
     n = degree + 1; /* omit trailing zeroes */
+    SEXP r;
     if(degree >= 1) {
+	SEXP zr, zi, rr, ri;
 	PROTECT(rr = allocVector(REALSXP, n));
 	PROTECT(ri = allocVector(REALSXP, n));
 	PROTECT(zr = allocVector(REALSXP, n));
 	PROTECT(zi = allocVector(REALSXP, n));
 	    if(!R_FINITE(pz[i].r) || !R_FINITE(pz[i].i))
 		error(_("invalid polynomial coefficient"));
 	    p_zr[degree-i] = pz[i].r;
 	    p_zi[degree-i] = pz[i].i;
+	}
+	bool fail;
 	R_cpolyroot(p_zr, p_zi, &degree, p_rr, p_ri, &fail);
 	if(fail) error(_("root finding code failed"));
 	UNPROTECT(2);
 	r = allocVector(CPLXSXP, degree);
 	Rcomplex *pr = COMPLEX(r);

The R Project SVN R

(root)/trunk/src/main/complex.c @ 63243 – Rev 88828 → 88866