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



	    }
#endif
	}
	return ans;
    default:
	error("illegal complex unary operator");
	return R_NilValue;/* -Wall*/
    }
}
 
static void complex_div(complex *c, complex *a, complex *b)

	    else complex_pow(&COMPLEX(ans)[i], &x1, &x2);
#endif
	}
	break;
    default:
	error("unimplemented complex operation");
    }
    /* Copy attributes from longest argument. */
    if (n1 > n2)
	copyMostAttrib(s1, ans);
    else if (n1 == n2) {

	    }
	    break;
	}
	UNPROTECT(1);
    }
    else errorcall(call, "non-numeric argument to function");
    PROTECT(x);
    PROTECT(y);
    ATTRIB(y) = duplicate(ATTRIB(x));
    OBJECT(y) = OBJECT(x);
    UNPROTECT(2);

    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;
 

	MATH1(40, lgammafn);
	MATH1(41, gammafn);
#endif
 
    default:
	errorcall(call, "unimplemented complex function");
    }
    if (naflag) warning("NAs produced in function \"%s\"", PRIMNAME(op));
    return y;
}
 

    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");
	return call;/* just for -Wall */
    }
}
 
SEXP do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)

    /* complex(length, real, imaginary) */
    SEXP ans, re, im;
    int i, na, nr, ni;
    na = asInteger(CAR(args));
    if(na == NA_INTEGER || na < 0)
	errorcall(call, "invalid length");
    PROTECT(re = coerceVector(CADR(args), REALSXP));
    PROTECT(im = coerceVector(CADDR(args), REALSXP));
    nr = length(re);
    ni = length(im);
    /* is always true: if (na >= 0) {*/

    case INTSXP:
    case LGLSXP:
	PROTECT(z = coerceVector(z, CPLXSXP));
	break;
    default:
	errorcall(call, "invalid argument type");
    }
    n = length(z);
    degree = n - 1;
    if(degree >= 1) {
	if(n > 49) errorcall(call, "polynomial degree too high (49 max)");
	/* <==>	 #define NMAX 50  in  ../appl/cpoly.c */
 
	if(COMPLEX(z)[n-1].r == 0.0 && COMPLEX(z)[n-1].i == 0.0)
	    errorcall(call, "highest power has coefficient 0");
 
	PROTECT(rr = allocVector(REALSXP, n));
	PROTECT(ri = allocVector(REALSXP, n));
	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");
	    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);
	if(fail) errorcall(call, "root finding code failed");
	UNPROTECT(2);
	r = allocVector(CPLXSXP, degree);
	for(i=0 ; i<n ; i++) {
	    COMPLEX(r)[i].r = REAL(rr)[i];
	    COMPLEX(r)[i].i = REAL(ri)[i];

Rev 5458	Rev 5731
Line 58...	Line 58...
58	`}`	58	`}`
59	`#endif`	59	`#endif`
60	`}`	60	`}`
61	`return ans;`	61	`return ans;`
62	`default:`	62	`default:`
63	`error("illegal complex unary operator\n");`	63	`error("illegal complex unary operator");`
64	`return R_NilValue;/* -Wall*/`	64	`return R_NilValue;/* -Wall*/`
65	`}`	65	`}`
66	`}`	66	`}`
67		67
68	`static void complex_div(complex c, complex a, complex *b)`	68	`static void complex_div(complex c, complex a, complex *b)`
Line 228...	Line 228...
228	`else complex_pow(&COMPLEX(ans)[i], &x1, &x2);`	228	`else complex_pow(&COMPLEX(ans)[i], &x1, &x2);`
229	`#endif`	229	`#endif`
230	`}`	230	`}`
231	`break;`	231	`break;`
232	`default:`	232	`default:`
233	`error("unimplemented complex operation\n");`	233	`error("unimplemented complex operation");`
234	`}`	234	`}`
235	`/* Copy attributes from longest argument. */`	235	`/* Copy attributes from longest argument. */`
236	`if (n1 > n2)`	236	`if (n1 > n2)`
237	`copyMostAttrib(s1, ans);`	237	`copyMostAttrib(s1, ans);`
238	`else if (n1 == n2) {`	238	`else if (n1 == n2) {`
Line 349...	Line 349...
349	`}`	349	`}`
350	`break;`	350	`break;`
351	`}`	351	`}`
352	`UNPROTECT(1);`	352	`UNPROTECT(1);`
353	`}`	353	`}`
354	`else errorcall(call, "non-numeric argument to function\n");`	354	`else errorcall(call, "non-numeric argument to function");`
355	`PROTECT(x);`	355	`PROTECT(x);`
356	`PROTECT(y);`	356	`PROTECT(y);`
357	`ATTRIB(y) = duplicate(ATTRIB(x));`	357	`ATTRIB(y) = duplicate(ATTRIB(x));`
358	`OBJECT(y) = OBJECT(x);`	358	`OBJECT(y) = OBJECT(x);`
359	`UNPROTECT(2);`	359	`UNPROTECT(2);`
Line 585...	Line 585...
585	`case 10002: cmath1(z_atan, COMPLEX(x), COMPLEX(y), n); break;`	585	`case 10002: cmath1(z_atan, COMPLEX(x), COMPLEX(y), n); break;`
586	`case 10003: cmath1(z_log, COMPLEX(x), COMPLEX(y), n); break;`	586	`case 10003: cmath1(z_log, COMPLEX(x), COMPLEX(y), n); break;`
587		587
588	`case 0:`	588	`case 0:`
589	`errorcall(call,`	589	`errorcall(call,`
590	`"abs() unimplemented for complex; use Mod()\n");`	590	`"abs() unimplemented for complex; use Mod()");`
591	`break;`	591	`break;`
592	`case 3: cmath1(z_sqrt, COMPLEX(x), COMPLEX(y), n); break;`	592	`case 3: cmath1(z_sqrt, COMPLEX(x), COMPLEX(y), n); break;`
593		593
594	`case 10: cmath1(z_exp, COMPLEX(x), COMPLEX(y), n); break;`	594	`case 10: cmath1(z_exp, COMPLEX(x), COMPLEX(y), n); break;`
595		595
Line 610...	Line 610...
610	`MATH1(40, lgammafn);`	610	`MATH1(40, lgammafn);`
611	`MATH1(41, gammafn);`	611	`MATH1(41, gammafn);`
612	`#endif`	612	`#endif`
613		613
614	`default:`	614	`default:`
615	`errorcall(call, "unimplemented complex function\n");`	615	`errorcall(call, "unimplemented complex function");`
616	`}`	616	`}`
617	`if (naflag) warning("NAs produced in function \"%s\"", PRIMNAME(op));`	617	`if (naflag) warning("NAs produced in function \"%s\"", PRIMNAME(op));`
618	`return y;`	618	`return y;`
619	`}`	619	`}`
620		620
Line 690...	Line 690...
690	`case 10004:`	690	`case 10004:`
691	`return cmath2(op, CAR(args), CADR(args), z_prec);`	691	`return cmath2(op, CAR(args), CADR(args), z_prec);`
692	`case 0:`	692	`case 0:`
693	`return cmath2(op, CAR(args), CADR(args), z_atan2);`	693	`return cmath2(op, CAR(args), CADR(args), z_atan2);`
694	`default:`	694	`default:`
695	`errorcall(call, "unimplemented complex function\n");`	695	`errorcall(call, "unimplemented complex function");`
696	`return call;/* just for -Wall */`	696	`return call;/* just for -Wall */`
697	`}`	697	`}`
698	`}`	698	`}`
699		699
700	`SEXP do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)`	700	`SEXP do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)`
Line 702...	Line 702...
702	`/* complex(length, real, imaginary) */`	702	`/* complex(length, real, imaginary) */`
703	`SEXP ans, re, im;`	703	`SEXP ans, re, im;`
704	`int i, na, nr, ni;`	704	`int i, na, nr, ni;`
705	`na = asInteger(CAR(args));`	705	`na = asInteger(CAR(args));`
706	`if(na == NA_INTEGER \|\| na < 0)`	706	`if(na == NA_INTEGER \|\| na < 0)`
707	`errorcall(call, "invalid length\n");`	707	`errorcall(call, "invalid length");`
708	`PROTECT(re = coerceVector(CADR(args), REALSXP));`	708	`PROTECT(re = coerceVector(CADR(args), REALSXP));`
709	`PROTECT(im = coerceVector(CADDR(args), REALSXP));`	709	`PROTECT(im = coerceVector(CADDR(args), REALSXP));`
710	`nr = length(re);`	710	`nr = length(re);`
711	`ni = length(im);`	711	`ni = length(im);`
712	`/* is always true: if (na >= 0) {*/`	712	`/* is always true: if (na >= 0) {*/`
Line 745...	Line 745...
745	`case INTSXP:`	745	`case INTSXP:`
746	`case LGLSXP:`	746	`case LGLSXP:`
747	`PROTECT(z = coerceVector(z, CPLXSXP));`	747	`PROTECT(z = coerceVector(z, CPLXSXP));`
748	`break;`	748	`break;`
749	`default:`	749	`default:`
750	`errorcall(call, "invalid argument type\n");`	750	`errorcall(call, "invalid argument type");`
751	`}`	751	`}`
752	`n = length(z);`	752	`n = length(z);`
753	`degree = n - 1;`	753	`degree = n - 1;`
754	`if(degree >= 1) {`	754	`if(degree >= 1) {`
755	`if(n > 49) errorcall(call, "polynomial degree too high (49 max)\n");`	755	`if(n > 49) errorcall(call, "polynomial degree too high (49 max)");`
756	`/* <==> #define NMAX 50 in ../appl/cpoly.c */`	756	`/* <==> #define NMAX 50 in ../appl/cpoly.c */`
757		757
758	`if(COMPLEX(z)[n-1].r == 0.0 && COMPLEX(z)[n-1].i == 0.0)`	758	`if(COMPLEX(z)[n-1].r == 0.0 && COMPLEX(z)[n-1].i == 0.0)`
759	`errorcall(call, "highest power has coefficient 0\n");`	759	`errorcall(call, "highest power has coefficient 0");`
760		760
761	`PROTECT(rr = allocVector(REALSXP, n));`	761	`PROTECT(rr = allocVector(REALSXP, n));`
762	`PROTECT(ri = allocVector(REALSXP, n));`	762	`PROTECT(ri = allocVector(REALSXP, n));`
763	`PROTECT(zr = allocVector(REALSXP, n));`	763	`PROTECT(zr = allocVector(REALSXP, n));`
764	`PROTECT(zi = allocVector(REALSXP, n));`	764	`PROTECT(zi = allocVector(REALSXP, n));`
765		765
766	`for(i=0 ; i<n ; i++) {`	766	`for(i=0 ; i<n ; i++) {`
767	`if(!R_FINITE(COMPLEX(z)[i].r) \|\| !R_FINITE(COMPLEX(z)[i].i))`	767	`if(!R_FINITE(COMPLEX(z)[i].r) \|\| !R_FINITE(COMPLEX(z)[i].i))`
768	`errorcall(call, "invalid polynomial coefficient\n");`	768	`errorcall(call, "invalid polynomial coefficient");`
769	`REAL(zr)[degree-i] = COMPLEX(z)[i].r;`	769	`REAL(zr)[degree-i] = COMPLEX(z)[i].r;`
770	`REAL(zi)[degree-i] = COMPLEX(z)[i].i;`	770	`REAL(zi)[degree-i] = COMPLEX(z)[i].i;`
771	`}`	771	`}`
772	`F77_SYMBOL(cpoly)(REAL(zr), REAL(zi), &degree,`	772	`F77_SYMBOL(cpoly)(REAL(zr), REAL(zi), &degree,`
773	`REAL(rr), REAL(ri), &fail);`	773	`REAL(rr), REAL(ri), &fail);`
774	`if(fail) errorcall(call, "root finding code failed\n");`	774	`if(fail) errorcall(call, "root finding code failed");`
775	`UNPROTECT(2);`	775	`UNPROTECT(2);`
776	`r = allocVector(CPLXSXP, degree);`	776	`r = allocVector(CPLXSXP, degree);`
777	`for(i=0 ; i<n ; i++) {`	777	`for(i=0 ; i<n ; i++) {`
778	`COMPLEX(r)[i].r = REAL(rr)[i];`	778	`COMPLEX(r)[i].r = REAL(rr)[i];`
779	`COMPLEX(r)[i].i = REAL(ri)[i];`	779	`COMPLEX(r)[i].i = REAL(ri)[i];`

The R Project SVN R

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