The R Project SVN R

Rev

Rev 24539 | Rev 31908 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 r 1 
/*
831 maechler 2 
 *  R : A Computer Language for Statistical Data Analysis
2 r 3 
 *  Copyright (C) 1995, 1996, 1997  Robert Gentleman and Ross Ihaka
7609 maechler 4 
 *  Copyright (C) 2000		    The R Development Core Team.
2 r 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
10 
 *
11 
 *  This program is distributed in the hope that it will be useful,
12 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 
 *  GNU General Public License for more details.
15 
 *
16 
 *  You should have received a copy of the GNU General Public License
17 
 *  along with this program; if not, write to the Free Software
5458 ripley 18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
2 r 19 
 */
20 
 
5187 hornik 21 
#ifdef HAVE_CONFIG_H
7701 hornik 22 
#include <config.h>
5187 hornik 23 
#endif
24 
 
15508 pd 25 
#include <Defn.h>		/* -> ../include/R_ext/Complex.h */
11499 ripley 26 
#include <Rmath.h>
15508 pd 27 
#include <R_ext/Applic.h>	/* R_cpoly */
2 r 28 
 
9417 ripley 29 
#include "arithmetic.h"		/* complex_*  */
2632 maechler 30 
 
9417 ripley 31 
#ifndef HAVE_HYPOT
32 
# define hypot pythag
33 
#endif
34 
 
10718 maechler 35 
SEXP complex_unary(ARITHOP_TYPE code, SEXP s1)
2 r 36 
{
580 ihaka 37 
    int i, n;
6994 pd 38 
    Rcomplex x;
580 ihaka 39 
    SEXP ans;
2 r 40 
 
580 ihaka 41 
    switch(code) {
42 
    case PLUSOP:
43 
	return s1;
44 
    case MINUSOP:
9147 maechler 45 
 
580 ihaka 46 
	ans = duplicate(s1);
47 
	n = LENGTH(s1);
48 
	for (i = 0; i < n; i++) {
49 
	    x = COMPLEX(s1)[i];
15508 pd 50 
#ifdef IEEE_754
580 ihaka 51 
	    COMPLEX(ans)[i].r = -x.r;
52 
	    COMPLEX(ans)[i].i = -x.i;
53 
#else
54 
	    if(ISNAN(x.r) || ISNAN(x.i)) {
55 
		COMPLEX(ans)[i].r = NA_REAL;
56 
		COMPLEX(ans)[i].i = NA_REAL;
57 
	    }
58 
	    else {
59 
		COMPLEX(ans)[i].r = -x.r;
60 
		COMPLEX(ans)[i].i = -x.i;
61 
	    }
62 
#endif
2 r 63 
	}
580 ihaka 64 
	return ans;
65 
    default:
9147 maechler 66 
	error_return("illegal complex unary operator");
580 ihaka 67 
    }
2 r 68 
}
69 
 
6994 pd 70 
static void complex_div(Rcomplex *c, Rcomplex *a, Rcomplex *b)
2 r 71 
{
580 ihaka 72 
    double ratio, den;
73 
    double abr, abi;
2 r 74 
 
580 ihaka 75 
    if( (abr = b->r) < 0)
76 
	abr = - abr;
77 
    if( (abi = b->i) < 0)
78 
	abi = - abi;
79 
    if( abr <= abi ) {
80 
#ifndef IEEE_754
81 
	if(abi == 0) {
82 
	    c->r = NA_REAL;
83 
	    c->i = NA_REAL;
84 
	    return;
2 r 85 
	}
580 ihaka 86 
#endif
87 
	ratio = b->r / b->i ;
88 
	den = b->i * (1 + ratio*ratio);
89 
	c->r = (a->r*ratio + a->i) / den;
90 
	c->i = (a->i*ratio - a->r) / den;
91 
    }
92 
    else {
93 
	ratio = b->i / b->r ;
94 
	den = b->r * (1 + ratio*ratio);
95 
	c->r = (a->r + a->i*ratio) / den;
96 
	c->i = (a->i - a->r*ratio) / den;
97 
    }
2 r 98 
}
99 
 
6994 pd 100 
static void complex_pow(Rcomplex *r, Rcomplex *a, Rcomplex *b)
2 r 101 
{
7609 maechler 102 
/* r := a^b */
580 ihaka 103 
    double logr, logi, x, y;
2212 maechler 104 
    int ib;
2 r 105 
 
7609 maechler 106 
    if(b->i == 0.) {		/* ^ "real" : be fast (and more accurate)*/
6098 pd 107 
	if(b->r == 1.) {	/* a^1 */
108 
	    r->r = a->r; r->i = a->i; return;
109 
	}
7609 maechler 110 
	if(a->i == 0. && a->r >= 0.) {
111 
	    r->r = R_pow(a->r, b->r); r->i = 0.; return;
6098 pd 112 
	}
2212 maechler 113 
	if(a->r == 0. && b->r == (ib = (int)b->r)) {/* (|a|*i)^b */
7867 ripley 114 
	    x = R_pow_di(a->i, ib);
6098 pd 115 
	    if(ib % 2) {	/* ib is odd ==> imaginary r */
2212 maechler 116 
		r->r = 0.;
117 
		r->i = ((ib>0 && ib %4 == 3)||(ib<0 && (-ib)%4 == 1))? -x : x;
6098 pd 118 
	    } else {		/* even exponent b : real r */
2212 maechler 119 
		r->r = (ib %4)? -x : x; r->i = 0.;
120 
	    }
121 
	    return;
122 
	}
123 
    }
580 ihaka 124 
    logr = log(hypot(a->r, a->i) );
125 
    logi = atan2(a->i, a->r);
126 
    x = exp( logr * b->r - logi * b->i );
127 
    y = logr * b->i + logi * b->r;
128 
    r->r = x * cos(y);
129 
    r->i = x * sin(y);
2 r 130 
}
131 
 
1839 ihaka 132 
/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */
133 
 
10718 maechler 134 
SEXP complex_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2)
2 r 135 
{
580 ihaka 136 
    int i, n, n1, n2;
6994 pd 137 
    Rcomplex x1, x2;
580 ihaka 138 
    SEXP ans;
16613 ripley 139 
 
140 
    /* Note: "s1" and "s1" are protected in the calling code. */
580 ihaka 141 
    n1 = LENGTH(s1);
142 
    n2 = LENGTH(s2);
16613 ripley 143 
     /* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */
144 
    if (n1 == 0 || n2 == 0) return(allocVector(CPLXSXP, 0));
145 
 
580 ihaka 146 
    n = (n1 > n2) ? n1 : n2;
147 
    ans = allocVector(CPLXSXP, n);
16613 ripley 148 
/*    if (n1 < 1 || n2 < 1) {
580 ihaka 149 
	for (i = 0; i < n; i++) {
150 
	    COMPLEX(ans)[i].r = NA_REAL;
151 
	    COMPLEX(ans)[i].i = NA_REAL;
2 r 152 
	}
580 ihaka 153 
	return ans;
154 
    }
16613 ripley 155 
*/
156 
 
580 ihaka 157 
    switch (code) {
158 
    case PLUSOP:
159 
	for (i = 0; i < n; i++) {
160 
	    x1 = COMPLEX(s1)[i % n1];
161 
	    x2 = COMPLEX(s2)[i % n2];
162 
#ifdef IEEE_754
163 
	    COMPLEX(ans)[i].r = x1.r + x2.r;
164 
	    COMPLEX(ans)[i].i = x1.i + x2.i;
165 
#else
166 
	    if (ISNAN(x1.r) || ISNAN(x1.i) || ISNAN(x2.r) || ISNAN(x2.i)) {
167 
		COMPLEX(ans)[i].r = NA_REAL;
168 
		COMPLEX(ans)[i].i = NA_REAL;
169 
	    }
170 
	    else {
171 
		COMPLEX(ans)[i].r = MATH_CHECK(x1.r + x2.r);
172 
		COMPLEX(ans)[i].i = MATH_CHECK(x1.i + x2.i);
173 
	    }
174 
#endif
2 r 175 
	}
580 ihaka 176 
	break;
177 
    case MINUSOP:
178 
	for (i = 0; i < n; i++) {
179 
	    x1 = COMPLEX(s1)[i % n1];
180 
	    x2 = COMPLEX(s2)[i % n2];
181 
#ifdef IEEE_754
182 
	    COMPLEX(ans)[i].r = x1.r - x2.r;
183 
	    COMPLEX(ans)[i].i = x1.i - x2.i;
184 
#else
185 
	    if (ISNAN(x1.r) || ISNAN(x1.i) || ISNAN(x2.r) || ISNAN(x2.i)) {
186 
		COMPLEX(ans)[i].r = NA_REAL;
187 
		COMPLEX(ans)[i].i = NA_REAL;
188 
	    }
189 
	    else {
190 
		COMPLEX(ans)[i].r = MATH_CHECK(x1.r - x2.r);
191 
		COMPLEX(ans)[i].i = MATH_CHECK(x1.i - x2.i);
192 
	    }
193 
#endif
194 
	}
195 
	break;
196 
    case TIMESOP:
197 
	for (i = 0; i < n; i++) {
198 
	    x1 = COMPLEX(s1)[i % n1];
199 
	    x2 = COMPLEX(s2)[i % n2];
200 
#ifdef IEEE_754
201 
	    COMPLEX(ans)[i].r = x1.r * x2.r - x1.i * x2.i;
202 
	    COMPLEX(ans)[i].i = x1.r * x2.i + x1.i * x2.r;
203 
#else
204 
	    if (ISNAN(x1.r) || ISNAN(x1.i) || ISNAN(x2.r) || ISNAN(x2.i)) {
205 
		COMPLEX(ans)[i].r = NA_REAL;
206 
		COMPLEX(ans)[i].i = NA_REAL;
207 
	    }
208 
	    else {
209 
		COMPLEX(ans)[i].r = MATH_CHECK(x1.r * x2.r - x1.i * x2.i);
210 
		COMPLEX(ans)[i].i = MATH_CHECK(x1.r * x2.i + x1.i * x2.r);
211 
	    }
212 
#endif
213 
	}
214 
	break;
215 
    case DIVOP:
216 
	for (i = 0; i < n; i++) {
217 
	    x1 = COMPLEX(s1)[i % n1];
218 
	    x2 = COMPLEX(s2)[i % n2];
219 
#ifdef IEEE_754
220 
	    complex_div(&COMPLEX(ans)[i], &x1, &x2);
221 
#else
222 
	    if (ISNAN(x1.r) || ISNAN(x1.i) || ISNAN(x2.r) || ISNAN(x2.i)) {
223 
		COMPLEX(ans)[i].r = NA_REAL;
224 
		COMPLEX(ans)[i].i = NA_REAL;
225 
	    }
226 
	    else complex_div(&COMPLEX(ans)[i], &x1, &x2);
227 
#endif
228 
	}
229 
	break;
230 
    case POWOP:
231 
	for (i = 0; i < n; i++) {
232 
	    x1 = COMPLEX(s1)[i % n1];
233 
	    x2 = COMPLEX(s2)[i % n2];
234 
#ifdef IEEE_754
235 
	    complex_pow(&COMPLEX(ans)[i], &x1, &x2);
236 
#else
237 
	    if (ISNAN(x1.r) || ISNAN(x1.i) || ISNAN(x2.r) || ISNAN(x2.i)) {
238 
		COMPLEX(ans)[i].r = NA_REAL;
239 
		COMPLEX(ans)[i].i = NA_REAL;
240 
	    }
241 
	    else complex_pow(&COMPLEX(ans)[i], &x1, &x2);
242 
#endif
243 
	}
244 
	break;
245 
    default:
5731 ripley 246 
	error("unimplemented complex operation");
580 ihaka 247 
    }
16613 ripley 248 
 
24539 luke 249 
    /* quick return if there are no attributes */
250 
    if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue)
251 
	return ans;
252 
 
16613 ripley 253 
    /* Copy attributes from longer argument. */
1010 ihaka 254 
    if (n1 > n2)
255 
	copyMostAttrib(s1, ans);
256 
    else if (n1 == n2) {
257 
	copyMostAttrib(s2, ans);
258 
	copyMostAttrib(s1, ans);
259 
    }
260 
    else
261 
	copyMostAttrib(s2, ans);
580 ihaka 262 
    return ans;
2 r 263 
}
264 
 
1839 ihaka 265 
/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */
266 
 
2 r 267 
SEXP do_cmathfuns(SEXP call, SEXP op, SEXP args, SEXP env)
831 maechler 268 
{
6098 pd 269 
    SEXP x, y = R_NilValue;	/* -Wall*/
580 ihaka 270 
    int i, n;
2 r 271 
 
580 ihaka 272 
    checkArity(op, args);
26610 ripley 273 
    if (DispatchGroup("Complex", call, op, args, env, &x))
274 
        return x;
580 ihaka 275 
    x = CAR(args);
276 
    n = length(x);
277 
    if (isComplex(x)) {
278 
	switch(PRIMVAL(op)) {
279 
	case 1:	/* Re */
280 
	    y = allocVector(REALSXP, n);
281 
	    for(i=0 ; i<n ; i++)
282 
		REAL(y)[i] = COMPLEX(x)[i].r;
283 
	    break;
284 
	case 2:	/* Im */
285 
	    y = allocVector(REALSXP, n);
286 
	    for(i=0 ; i<n ; i++)
287 
		REAL(y)[i] = COMPLEX(x)[i].i;
288 
	    break;
289 
	case 3:	/* Mod */
6098 pd 290 
	case 6:	/* abs */
580 ihaka 291 
	    y = allocVector(REALSXP, n);
292 
	    for(i=0 ; i<n ; i++) {
293 
#ifdef IEEE_754
294 
		REAL(y)[i] = hypot(COMPLEX(x)[i].r, COMPLEX(x)[i].i);
295 
#else
296 
		if(ISNAN(COMPLEX(x)[i].r) || ISNAN(COMPLEX(x)[i].i)) {
297 
		    REAL(y)[i] = NA_REAL;
2 r 298 
		}
580 ihaka 299 
		else {
300 
		    REAL(y)[i] = hypot(COMPLEX(x)[i].r, COMPLEX(x)[i].i);
2 r 301 
		}
580 ihaka 302 
#endif
303 
	    }
304 
	    break;
305 
	case 4:	/* Arg */
306 
	    y = allocVector(REALSXP, n);
307 
	    for(i=0 ; i<n ; i++) {
308 
#ifdef IEEE_754
309 
		REAL(y)[i] = atan2(COMPLEX(x)[i].i, COMPLEX(x)[i].r);
310 
#else
311 
		if(ISNAN(COMPLEX(x)[i].r) || ISNAN(COMPLEX(x)[i].i)) {
312 
		    REAL(y)[i] = NA_REAL;
313 
		}
314 
		else {
315 
		    REAL(y)[i] = atan2(COMPLEX(x)[i].i, COMPLEX(x)[i].r);
316 
		}
317 
#endif
318 
	    }
319 
	    break;
320 
	case 5:	/* Conj */
321 
	    y = allocVector(CPLXSXP, n);
322 
	    for(i=0 ; i<n ; i++) {
323 
#ifdef IEEE_754
324 
		COMPLEX(y)[i].r = COMPLEX(x)[i].r;
325 
		COMPLEX(y)[i].i = -COMPLEX(x)[i].i;
326 
#else
327 
		if(ISNAN(COMPLEX(x)[i].r) || ISNAN(COMPLEX(x)[i].i)) {
328 
		    COMPLEX(y)[i].r = NA_REAL;
329 
		    COMPLEX(y)[i].i = NA_REAL;
330 
		}
331 
		else {
332 
		    COMPLEX(y)[i].r = COMPLEX(x)[i].r;
333 
		    COMPLEX(y)[i].i = -COMPLEX(x)[i].i;
334 
		}
335 
#endif
336 
	    }
337 
	    break;
2 r 338 
	}
580 ihaka 339 
    }
340 
    else if(isNumeric(x)) {
341 
	if(isReal(x)) PROTECT(x);
342 
	else PROTECT(x = coerceVector(x, REALSXP));
343 
	switch(PRIMVAL(op)) {
344 
	case 1:	/* Re */
345 
	case 5:	/* Conj */
346 
	    y = allocVector(REALSXP, n);
347 
	    for(i=0 ; i<n ; i++)
348 
		REAL(y)[i] = REAL(x)[i];
349 
	    break;
350 
	case 2:	/* Im */
351 
	case 4:	/* Arg */
352 
	    y = allocVector(REALSXP, n);
353 
	    for(i=0 ; i<n ; i++)
354 
		if(ISNAN(REAL(x)[i]))
355 
		    REAL(y)[i] = REAL(x)[i];
356 
		else
357 
		    REAL(y)[i] = 0;
358 
	    break;
359 
	case 3:	/* Mod */
6098 pd 360 
	case 6:	/* abs */
580 ihaka 361 
	    y = allocVector(REALSXP, n);
362 
	    for(i=0 ; i<n ; i++) {
363 
#ifdef IEEE_754
364 
		REAL(y)[i] = fabs(REAL(x)[i]);
365 
#else
366 
		if(ISNAN(REAL(x)[i]))
367 
		    REAL(y)[i] = REAL(x)[i];
368 
		else
369 
		    REAL(y)[i] = fabs(REAL(x)[i]);
370 
#endif
371 
	    }
372 
	    break;
373 
	}
374 
	UNPROTECT(1);
375 
    }
5731 ripley 376 
    else errorcall(call, "non-numeric argument to function");
580 ihaka 377 
    PROTECT(x);
378 
    PROTECT(y);
10172 luke 379 
    SET_ATTRIB(y, duplicate(ATTRIB(x)));
380 
    SET_OBJECT(y, OBJECT(x));
580 ihaka 381 
    UNPROTECT(2);
382 
    return y;
2 r 383 
}
384 
 
6994 pd 385 
static void z_rround(Rcomplex *r, Rcomplex *x, Rcomplex *p)
2 r 386 
{
580 ihaka 387 
    r->r = rround(x->r, p->r);
388 
    r->i = rround(x->i, p->r);
2 r 389 
}
390 
 
6098 pd 391 
/* Question:  This treats real and imaginary parts separately.  Should
392 
   it do them jointly? */
2 r 393 
 
6994 pd 394 
static void z_prec(Rcomplex *r, Rcomplex *x, Rcomplex *p)
2 r 395 
{
580 ihaka 396 
    r->r = prec(x->r, p->r);
397 
    r->i = prec(x->i, p->r);
2 r 398 
}
399 
 
6994 pd 400 
static void z_log(Rcomplex *r, Rcomplex *z)
2 r 401 
{
580 ihaka 402 
    r->i = atan2(z->i, z->r);
403 
    r->r = log(hypot( z->r, z->i ));
2 r 404 
}
405 
 
6994 pd 406 
static void z_logbase(Rcomplex *r, Rcomplex *z, Rcomplex *base)
2 r 407 
{
6994 pd 408 
    Rcomplex t1, t2;
580 ihaka 409 
    z_log(&t1, z);
410 
    z_log(&t2, base);
411 
    complex_div(r, &t1, &t2);
2 r 412 
}
413 
 
6994 pd 414 
static void z_exp(Rcomplex *r, Rcomplex *z)
2 r 415 
{
580 ihaka 416 
    double expx;
417 
    expx = exp(z->r);
831 maechler 418 
    r->r = expx * cos(z->i);
580 ihaka 419 
    r->i = expx * sin(z->i);
2 r 420 
}
421 
 
6994 pd 422 
static void z_sqrt(Rcomplex *r, Rcomplex *z)
2 r 423 
{
580 ihaka 424 
    double mag;
2 r 425 
 
580 ihaka 426 
    if( (mag = hypot(z->r, z->i)) == 0.0)
427 
	r->r = r->i = 0.0;
428 
    else if(z->r > 0) {
429 
	r->r = sqrt(0.5 * (mag + z->r) );
430 
	r->i = z->i / r->r / 2;
431 
    }
432 
    else {
433 
	r->i = sqrt(0.5 * (mag - z->r) );
434 
	if(z->i < 0)
435 
	    r->i = - r->i;
436 
	r->r = z->i / r->i / 2;
437 
    }
2 r 438 
}
439 
 
6994 pd 440 
static void z_cos(Rcomplex *r, Rcomplex *z)
2 r 441 
{
580 ihaka 442 
    r->r = cos(z->r) * cosh(z->i);
443 
    r->i = - sin(z->r) * sinh(z->i);
2 r 444 
}
445 
 
6994 pd 446 
static void z_sin(Rcomplex *r, Rcomplex *z)
2 r 447 
{
580 ihaka 448 
    r->r = sin(z->r) * cosh(z->i);
831 maechler 449 
    r->i = cos(z->r) * sinh(z->i);
2 r 450 
}
451 
 
6994 pd 452 
static void z_tan(Rcomplex *r, Rcomplex *z)
2 r 453 
{
580 ihaka 454 
    double x2, y2, den;
455 
    x2 = 2.0 * z->r;
456 
    y2 = 2.0 * z->i;
457 
    den = cos(x2) + cosh(y2);
458 
    r->r = sin(x2)/den;
459 
    r->i = sinh(y2)/den;
2 r 460 
}
461 
 
1839 ihaka 462 
	/* Complex Arcsin and Arccos Functions */
2 r 463 
	/* Equation (4.4.37) Abramowitz and Stegun */
831 maechler 464 
 
6994 pd 465 
static void z_asin(Rcomplex *r, Rcomplex *z)
2 r 466 
{
10866 maechler 467 
    double alpha, bet, t1, t2, x, y;
580 ihaka 468 
    x = z->r;
469 
    y = z->i;
9307 maechler 470 
    t1 = 0.5 * hypot(x + 1, y);
471 
    t2 = 0.5 * hypot(x - 1, y);
580 ihaka 472 
    alpha = t1 + t2;
10866 maechler 473 
    bet = t1 - t2;
474 
    r->r = asin(bet);
580 ihaka 475 
    r->i = log(alpha + sqrt(alpha*alpha - 1));
2 r 476 
}
477 
 
6994 pd 478 
static void z_acos(Rcomplex *r, Rcomplex *z)
2 r 479 
{
10866 maechler 480 
    Rcomplex Asin;
481 
    z_asin(&Asin, z);
482 
    r->r = M_PI_2 - Asin.r;
483 
    r->i = - Asin.i;
2 r 484 
}
485 
 
486 
	/* Complex Arctangent Function */
487 
	/* Equation (4.4.39) Abramowitz and Stegun */
488 
 
6994 pd 489 
static void z_atan(Rcomplex *r, Rcomplex *z)
2 r 490 
{
580 ihaka 491 
    double x, y;
492 
    x = z->r;
493 
    y = z->i;
494 
    r->r = 0.5 * atan(2 * x / ( 1 - x * x - y * y));
495 
    r->i = 0.25 * log((x * x + (y + 1) * (y + 1)) /
496 
		      (x * x + (y - 1) * (y - 1)));
2 r 497 
}
498 
 
6994 pd 499 
static void z_atan2(Rcomplex *r, Rcomplex *csn, Rcomplex *ccs)
2 r 500 
{
6994 pd 501 
    Rcomplex tmp;
580 ihaka 502 
    if (ccs->r == 0 && ccs->i == 0) {
503 
	if(csn->r == 0 && csn->r == 0) {
504 
	    r->r = NA_REAL;
505 
	    r->i = NA_REAL;
2 r 506 
	}
507 
	else {
7527 maechler 508 
	    r->r = fsign(M_PI_2, csn->r);
831 maechler 509 
	    r->i = 0;
580 ihaka 510 
	}
511 
    }
512 
    else {
513 
	complex_div(&tmp, csn, ccs);
514 
	z_atan(r, &tmp);
2487 maechler 515 
	if(ccs->r < 0) r->r += M_PI;
516 
	if(r->r > M_PI) r->r -= 2 * M_PI;
831 maechler 517 
    }
2 r 518 
}
519 
 
6994 pd 520 
static void z_acosh(Rcomplex *r, Rcomplex *z)
2 r 521 
{
6994 pd 522 
    Rcomplex a;
580 ihaka 523 
    z_acos(&a, z);
524 
    r->r = -a.i;
525 
    r->i = a.r;
2 r 526 
}
527 
 
6994 pd 528 
static void z_asinh(Rcomplex *r, Rcomplex *z)
2 r 529 
{
6994 pd 530 
    Rcomplex a, b;
580 ihaka 531 
    b.r = -z->i;
532 
    b.i =  z->r;
533 
    z_asin(&a, &b);
534 
    r->r =  a.i;
535 
    r->i = -a.r;
2 r 536 
}
537 
 
6994 pd 538 
static void z_atanh(Rcomplex *r, Rcomplex *z)
2 r 539 
{
6994 pd 540 
    Rcomplex a, b;
580 ihaka 541 
    b.r = -z->i;
542 
    b.i =  z->r;
543 
    z_atan(&a, &b);
544 
    r->r =  a.i;
545 
    r->i = -a.r;
2 r 546 
}
547 
 
6994 pd 548 
static void z_cosh(Rcomplex *r, Rcomplex *z)
2 r 549 
{
6994 pd 550 
    Rcomplex a;
580 ihaka 551 
    a.r = -z->i;
552 
    a.i =  z->r;
553 
    z_cos(r, &a);
2 r 554 
}
555 
 
6994 pd 556 
static void z_sinh(Rcomplex *r, Rcomplex *z)
2 r 557 
{
6994 pd 558 
    Rcomplex a, b;
580 ihaka 559 
    b.r = -z->i;
560 
    b.i =  z->r;
561 
    z_sin(&a, &b);
562 
    r->r =  a.i;
563 
    r->i = -a.r;
2 r 564 
}
565 
 
6994 pd 566 
static void z_tanh(Rcomplex *r, Rcomplex *z)
2 r 567 
{
6994 pd 568 
    Rcomplex a, b;
580 ihaka 569 
    b.r = -z->i;
570 
    b.i =  z->r;
571 
    z_tan(&a, &b);
572 
    r->r =  a.i;
573 
    r->i = -a.r;
2 r 574 
}
575 
 
19912 duncan 576 
static Rboolean cmath1(void (*f)(), Rcomplex *x, Rcomplex *y, int n)
2 r 577 
{
831 maechler 578 
    int i;
19912 duncan 579 
    Rboolean naflag = FALSE;
1839 ihaka 580 
    for (i = 0 ; i < n ; i++) {
580 ihaka 581 
	if (ISNA(x[i].r) || ISNA(x[i].i)) {
582 
	    y[i].r = NA_REAL;
583 
	    y[i].i = NA_REAL;
2 r 584 
	}
580 ihaka 585 
	else {
586 
	    f(&y[i], &x[i]);
587 
#ifndef IEEE_754
588 
	    if(ISNA(y[i].r) || ISNA(y[i].i)) {
589 
		y[i].r = NA_REAL;
590 
		y[i].i = NA_REAL;
19912 duncan 591 
		naflag = TRUE;
580 ihaka 592 
	    }
593 
#endif
594 
	}
595 
    }
19912 duncan 596 
 
597 
    return(naflag);
2 r 598 
}
599 
 
600 
SEXP complex_math1(SEXP call, SEXP op, SEXP args, SEXP env)
601 
{
580 ihaka 602 
    SEXP x, y;
603 
    int n;
19912 duncan 604 
    Rboolean naflag = FALSE;
6098 pd 605 
    PROTECT(x = CAR(args));
580 ihaka 606 
    n = length(x);
6098 pd 607 
    PROTECT(y = allocVector(CPLXSXP, n));
2 r 608 
 
580 ihaka 609 
    switch (PRIMVAL(op)) {
19912 duncan 610 
    case 10002: naflag = cmath1(z_atan, COMPLEX(x), COMPLEX(y), n); break;
611 
    case 10003: naflag = cmath1(z_log, COMPLEX(x), COMPLEX(y), n); break;
2 r 612 
 
19912 duncan 613 
    case 3: naflag = cmath1(z_sqrt, COMPLEX(x), COMPLEX(y), n); break;
2 r 614 
 
19912 duncan 615 
    case 10: naflag = cmath1(z_exp, COMPLEX(x), COMPLEX(y), n); break;
2 r 616 
 
19912 duncan 617 
    case 20: naflag = cmath1(z_cos, COMPLEX(x), COMPLEX(y), n); break;
618 
    case 21: naflag = cmath1(z_sin, COMPLEX(x), COMPLEX(y), n); break;
619 
    case 22: naflag = cmath1(z_tan, COMPLEX(x), COMPLEX(y), n); break;
620 
    case 23: naflag = cmath1(z_acos, COMPLEX(x), COMPLEX(y), n); break;
621 
    case 24: naflag = cmath1(z_asin, COMPLEX(x), COMPLEX(y), n); break;
2 r 622 
 
19912 duncan 623 
    case 30: naflag = cmath1(z_cosh, COMPLEX(x), COMPLEX(y), n); break;
624 
    case 31: naflag = cmath1(z_sinh, COMPLEX(x), COMPLEX(y), n); break;
625 
    case 32: naflag = cmath1(z_tanh, COMPLEX(x), COMPLEX(y), n); break;
626 
    case 33: naflag = cmath1(z_acosh, COMPLEX(x), COMPLEX(y), n); break;
627 
    case 34: naflag = cmath1(z_asinh, COMPLEX(x), COMPLEX(y), n); break;
628 
    case 35: naflag = cmath1(z_atanh, COMPLEX(x), COMPLEX(y), n); break;
2 r 629 
 
630 
#ifdef NOTYET
2278 maechler 631 
	MATH1(40, lgammafn);
632 
	MATH1(41, gammafn);
2 r 633 
#endif
634 
 
580 ihaka 635 
    default:
5731 ripley 636 
	errorcall(call, "unimplemented complex function");
580 ihaka 637 
    }
6098 pd 638 
    if (naflag)
639 
	warning("NAs produced in function \"%s\"", PRIMNAME(op));
10172 luke 640 
    SET_ATTRIB(y, duplicate(ATTRIB(x)));
641 
    SET_OBJECT(y, OBJECT(x));
6098 pd 642 
    UNPROTECT(2);
580 ihaka 643 
    return y;
2 r 644 
}
645 
 
1839 ihaka 646 
/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */
647 
 
2 r 648 
static SEXP cmath2(SEXP op, SEXP sa, SEXP sb, void (*f)())
649 
{
580 ihaka 650 
    int i, n, na, nb;
6994 pd 651 
    Rcomplex ai, bi, *a, *b, *y;
580 ihaka 652 
    SEXP sy;
19912 duncan 653 
    int naflag = 0;
580 ihaka 654 
    na = length(sa);
655 
    nb = length(sb);
6098 pd 656 
    if ((na == 0) || (nb == 0))
657 
	return(allocVector(CPLXSXP, 0));
580 ihaka 658 
    n = (na < nb) ? nb : na;
659 
    PROTECT(sa = coerceVector(sa, CPLXSXP));
660 
    PROTECT(sb = coerceVector(sb, CPLXSXP));
661 
    PROTECT(sy = allocVector(CPLXSXP, n));
662 
    a = COMPLEX(sa);
663 
    b = COMPLEX(sb);
664 
    y = COMPLEX(sy);
6098 pd 665 
    naflag = 0;
666 
    for (i = 0; i < n; i++) {
667 
	ai = a[i % na];
668 
	bi = b[i % nb];
669 
	if(ISNA(ai.r) && ISNA(ai.i) &&
670 
	   ISNA(bi.r) && ISNA(bi.i)) {
580 ihaka 671 
	    y[i].r = NA_REAL;
672 
	    y[i].i = NA_REAL;
2 r 673 
	}
6098 pd 674 
	else {
675 
	    f(&y[i], &ai, &bi);
676 
#ifndef IEEE_754
677 
	    if(ISNA(y[i].r) || ISNA(y[i].i)) {
580 ihaka 678 
		y[i].r = NA_REAL;
679 
		y[i].i = NA_REAL;
6098 pd 680 
		naflag = 1;
580 ihaka 681 
	    }
682 
#endif
2 r 683 
	}
580 ihaka 684 
    }
6098 pd 685 
    if (naflag)
686 
	warning("NAs produced in function \"%s\"", PRIMNAME(op));
580 ihaka 687 
    if(n == na) {
10172 luke 688 
	SET_ATTRIB(sy, duplicate(ATTRIB(sa)));
689 
	SET_OBJECT(sy, OBJECT(sa));
580 ihaka 690 
    }
691 
    else if(n == nb) {
10172 luke 692 
	SET_ATTRIB(sy, duplicate(ATTRIB(sb)));
693 
	SET_OBJECT(sy, OBJECT(sb));
580 ihaka 694 
    }
695 
    UNPROTECT(3);
696 
    return sy;
2 r 697 
}
698 
 
831 maechler 699 
	/* Complex Functions of Two Arguments */
580 ihaka 700 
 
2 r 701 
SEXP complex_math2(SEXP call, SEXP op, SEXP args, SEXP env)
702 
{
1839 ihaka 703 
    switch (PRIMVAL(op)) {
704 
    case 10001:
705 
	return cmath2(op, CAR(args), CADR(args), z_rround);
706 
    case 10002:
707 
	return cmath2(op, CAR(args), CADR(args), z_atan2);
708 
    case 10003:
709 
	return cmath2(op, CAR(args), CADR(args), z_logbase);
710 
    case 10004:
711 
	return cmath2(op, CAR(args), CADR(args), z_prec);
712 
    case 0:
713 
	return cmath2(op, CAR(args), CADR(args), z_atan2);
714 
    default:
9147 maechler 715 
	errorcall_return(call, "unimplemented complex function");
1839 ihaka 716 
    }
2 r 717 
}
718 
 
719 
SEXP do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)
720 
{
1839 ihaka 721 
    /* complex(length, real, imaginary) */
580 ihaka 722 
    SEXP ans, re, im;
723 
    int i, na, nr, ni;
724 
    na = asInteger(CAR(args));
725 
    if(na == NA_INTEGER || na < 0)
5731 ripley 726 
	errorcall(call, "invalid length");
580 ihaka 727 
    PROTECT(re = coerceVector(CADR(args), REALSXP));
728 
    PROTECT(im = coerceVector(CADDR(args), REALSXP));
729 
    nr = length(re);
730 
    ni = length(im);
1839 ihaka 731 
    /* is always true: if (na >= 0) {*/
732 
    na = (nr > na) ? nr : na;
733 
    na = (ni > na) ? ni : na;
734 
    /* }*/
580 ihaka 735 
    ans = allocVector(CPLXSXP, na);
736 
    for(i=0 ; i<na ; i++) {
737 
	COMPLEX(ans)[i].r = 0;
738 
	COMPLEX(ans)[i].i = 0;
739 
    }
740 
    UNPROTECT(2);
741 
    if(na > 0 && nr > 0) {
742 
	for(i=0 ; i<na ; i++)
743 
	    COMPLEX(ans)[i].r = REAL(re)[i%nr];
744 
    }
745 
    if(na > 0 && ni > 0) {
746 
	for(i=0 ; i<na ; i++)
747 
	    COMPLEX(ans)[i].i = REAL(im)[i%ni];
748 
    }
749 
    return ans;
2 r 750 
}
751 
 
580 ihaka 752 
 
2 r 753 
SEXP do_polyroot(SEXP call, SEXP op, SEXP args, SEXP rho)
754 
{
580 ihaka 755 
    SEXP z, zr, zi, r, rr, ri;
10866 maechler 756 
    Rboolean fail;
757 
    int degree, i, n;
758 
 
580 ihaka 759 
    checkArity(op, args);
760 
    z = CAR(args);
761 
    switch(TYPEOF(z)) {
762 
    case CPLXSXP:
763 
	PROTECT(z);
764 
	break;
765 
    case REALSXP:
766 
    case INTSXP:
767 
    case LGLSXP:
768 
	PROTECT(z = coerceVector(z, CPLXSXP));
769 
	break;
770 
    default:
5731 ripley 771 
	errorcall(call, "invalid argument type");
580 ihaka 772 
    }
773 
    n = length(z);
18490 ripley 774 
    degree = 0;
775 
    for(i = 0; i < n; i++) {
776 
	if(COMPLEX(z)[i].r!= 0.0 || COMPLEX(z)[i].i != 0.0) degree = i;
777 
    }
778 
    n = degree + 1; /* omit trailing zeroes */
580 ihaka 779 
    if(degree >= 1) {
5731 ripley 780 
	if(n > 49) errorcall(call, "polynomial degree too high (49 max)");
2212 maechler 781 
	/* <==>	 #define NMAX 50  in  ../appl/cpoly.c */
782 
 
18490 ripley 783 
	/* if(COMPLEX(z)[n-1].r == 0.0 && COMPLEX(z)[n-1].i == 0.0)
784 
	   errorcall(call, "highest power has coefficient 0");*/
2 r 785 
 
580 ihaka 786 
	PROTECT(rr = allocVector(REALSXP, n));
787 
	PROTECT(ri = allocVector(REALSXP, n));
788 
	PROTECT(zr = allocVector(REALSXP, n));
789 
	PROTECT(zi = allocVector(REALSXP, n));
2 r 790 
 
580 ihaka 791 
	for(i=0 ; i<n ; i++) {
5107 maechler 792 
	    if(!R_FINITE(COMPLEX(z)[i].r) || !R_FINITE(COMPLEX(z)[i].i))
6191 maechler 793 
		errorcall(call, "invalid polynomial coefficient");
580 ihaka 794 
	    REAL(zr)[degree-i] = COMPLEX(z)[i].r;
795 
	    REAL(zi)[degree-i] = COMPLEX(z)[i].i;
2 r 796 
	}
10885 maechler 797 
	R_cpolyroot(REAL(zr), REAL(zi), &degree, REAL(rr), REAL(ri), &fail);
5731 ripley 798 
	if(fail) errorcall(call, "root finding code failed");
580 ihaka 799 
	UNPROTECT(2);
800 
	r = allocVector(CPLXSXP, degree);
8987 luke 801 
	for(i=0 ; i<degree ; i++) {
580 ihaka 802 
	    COMPLEX(r)[i].r = REAL(rr)[i];
803 
	    COMPLEX(r)[i].i = REAL(ri)[i];
2 r 804 
	}
580 ihaka 805 
	UNPROTECT(3);
806 
    }
807 
    else {
808 
	UNPROTECT(1);
809 
	r = allocVector(CPLXSXP, 0);
810 
    }
811 
    return r;
2 r 812 
}