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



static int naflag;
 
SEXP complex_unary(int code, SEXP s1)
{
    int i, n;
    Rcomplex x;
    SEXP ans;
 
    switch(code) {
    case PLUSOP:
	return s1;

	error("illegal complex unary operator");
	return R_NilValue;	/* -Wall*/
    }
}
 
static void complex_div(Rcomplex *c, Rcomplex *a, Rcomplex *b)
{
    double ratio, den;
    double abr, abi;
 
    if( (abr = b->r) < 0)

	c->r = (a->r + a->i*ratio) / den;
	c->i = (a->i - a->r*ratio) / den;
    }
}
 
static void complex_pow(Rcomplex *r, Rcomplex *a, Rcomplex *b)
{
    double logr, logi, x, y;
    int ib;
 
    if(b->i == 0.) {		/* be fast (and more accurate)*/

/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */
 
SEXP complex_binary(int code, SEXP s1, SEXP s2)
{
    int i, n, n1, n2;
    Rcomplex x1, x2;
    SEXP ans;
    n1 = LENGTH(s1);
    n2 = LENGTH(s2);
    n = (n1 > n2) ? n1 : n2;
    /* Note: "s1" and "s1" are protected in the calling code. */

    OBJECT(y) = OBJECT(x);
    UNPROTECT(2);
    return y;
}
 
static void z_rround(Rcomplex *r, Rcomplex *x, Rcomplex *p)
{
    r->r = rround(x->r, p->r);
    r->i = rround(x->i, p->r);
}
 
/* Question:  This treats real and imaginary parts separately.  Should
   it do them jointly? */
 
static void z_prec(Rcomplex *r, Rcomplex *x, Rcomplex *p)
{
    r->r = prec(x->r, p->r);
    r->i = prec(x->i, p->r);
}
 
static void z_log(Rcomplex *r, Rcomplex *z)
{
    r->i = atan2(z->i, z->r);
    r->r = log(hypot( z->r, z->i ));
}
 
static void z_logbase(Rcomplex *r, Rcomplex *z, Rcomplex *base)
{
    Rcomplex t1, t2;
    z_log(&t1, z);
    z_log(&t2, base);
    complex_div(r, &t1, &t2);
}
 
static void z_exp(Rcomplex *r, Rcomplex *z)
{
    double expx;
    expx = exp(z->r);
    r->r = expx * cos(z->i);
    r->i = expx * sin(z->i);
}
 
static void z_sqrt(Rcomplex *r, Rcomplex *z)
{
    double mag;
 
    if( (mag = hypot(z->r, z->i)) == 0.0)
	r->r = r->i = 0.0;

	    r->i = - r->i;
	r->r = z->i / r->i / 2;
    }
}
 
static void z_cos(Rcomplex *r, Rcomplex *z)
{
    r->r = cos(z->r) * cosh(z->i);
    r->i = - sin(z->r) * sinh(z->i);
}
 
static void z_sin(Rcomplex *r, Rcomplex *z)
{
    r->r = sin(z->r) * cosh(z->i);
    r->i = cos(z->r) * sinh(z->i);
}
 
static void z_tan(Rcomplex *r, Rcomplex *z)
{
    double x2, y2, den;
    x2 = 2.0 * z->r;
    y2 = 2.0 * z->i;
    den = cos(x2) + cosh(y2);

}
 
	/* Complex Arcsin and Arccos Functions */
	/* Equation (4.4.37) Abramowitz and Stegun */
 
static void z_asin(Rcomplex *r, Rcomplex *z)
{
    double alpha, beta, t1, t2, x, y;
    x = z->r;
    y = z->i;
    t1 = 0.5 * sqrt((x + 1) * (x + 1) + y * y);

    beta = t1 - t2;
    r->r = asin(beta);
    r->i = log(alpha + sqrt(alpha*alpha - 1));
}
 
static void z_acos(Rcomplex *r, Rcomplex *z)
{
    Rcomplex asin;
    z_asin(&asin, z);
    r->r = M_PI_half - asin.r;
    r->i = - asin.i;
}
 
	/* Complex Arctangent Function */
	/* Equation (4.4.39) Abramowitz and Stegun */
 
static void z_atan(Rcomplex *r, Rcomplex *z)
{
    double x, y;
    x = z->r;
    y = z->i;
    r->r = 0.5 * atan(2 * x / ( 1 - x * x - y * y));
    r->i = 0.25 * log((x * x + (y + 1) * (y + 1)) /
		      (x * x + (y - 1) * (y - 1)));
}
 
static void z_atan2(Rcomplex *r, Rcomplex *csn, Rcomplex *ccs)
{
    Rcomplex tmp;
    if (ccs->r == 0 && ccs->i == 0) {
	if(csn->r == 0 && csn->r == 0) {
	    r->r = NA_REAL;
	    r->i = NA_REAL;
	}

	if(ccs->r < 0) r->r += M_PI;
	if(r->r > M_PI) r->r -= 2 * M_PI;
    }
}
 
static void z_acosh(Rcomplex *r, Rcomplex *z)
{
    Rcomplex a;
    z_acos(&a, z);
    r->r = -a.i;
    r->i = a.r;
}
 
static void z_asinh(Rcomplex *r, Rcomplex *z)
{
    Rcomplex a, b;
    b.r = -z->i;
    b.i =  z->r;
    z_asin(&a, &b);
    r->r =  a.i;
    r->i = -a.r;
}
 
static void z_atanh(Rcomplex *r, Rcomplex *z)
{
    Rcomplex a, b;
    b.r = -z->i;
    b.i =  z->r;
    z_atan(&a, &b);
    r->r =  a.i;
    r->i = -a.r;
}
 
static void z_cosh(Rcomplex *r, Rcomplex *z)
{
    Rcomplex a;
    a.r = -z->i;
    a.i =  z->r;
    z_cos(r, &a);
}
 
static void z_sinh(Rcomplex *r, Rcomplex *z)
{
    Rcomplex a, b;
    b.r = -z->i;
    b.i =  z->r;
    z_sin(&a, &b);
    r->r =  a.i;
    r->i = -a.r;
}
 
static void z_tanh(Rcomplex *r, Rcomplex *z)
{
    Rcomplex a, b;
    b.r = -z->i;
    b.i =  z->r;
    z_tan(&a, &b);
    r->r =  a.i;
    r->i = -a.r;
}
 
static void cmath1(void (*f)(), Rcomplex *x, Rcomplex *y, int n)
{
    int i;
    for (i = 0 ; i < n ; i++) {
	if (ISNA(x[i].r) || ISNA(x[i].i)) {
	    y[i].r = NA_REAL;

/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */
 
static SEXP cmath2(SEXP op, SEXP sa, SEXP sb, void (*f)())
{
    int i, n, na, nb;
    Rcomplex ai, bi, *a, *b, *y;
    SEXP sy;
 
    na = length(sa);
    nb = length(sb);
    if ((na == 0) || (nb == 0))

Rev 6191	Rev 6994
Line 31...	Line 31...
31	`static int naflag;`	31	`static int naflag;`
32		32
33	`SEXP complex_unary(int code, SEXP s1)`	33	`SEXP complex_unary(int code, SEXP s1)`
34	`{`	34	`{`
35	`int i, n;`	35	`int i, n;`
36	`complex x;`	36	`Rcomplex x;`
37	`SEXP ans;`	37	`SEXP ans;`
38		38
39	`switch(code) {`	39	`switch(code) {`
40	`case PLUSOP:`	40	`case PLUSOP:`
41	`return s1;`	41	`return s1;`
Line 63...	Line 63...
63	`error("illegal complex unary operator");`	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(Rcomplex c, Rcomplex a, Rcomplex *b)`
69	`{`	69	`{`
70	`double ratio, den;`	70	`double ratio, den;`
71	`double abr, abi;`	71	`double abr, abi;`
72		72
73	`if( (abr = b->r) < 0)`	73	`if( (abr = b->r) < 0)`
Line 93...	Line 93...
93	`c->r = (a->r + a->i*ratio) / den;`	93	`c->r = (a->r + a->i*ratio) / den;`
94	`c->i = (a->i - a->r*ratio) / den;`	94	`c->i = (a->i - a->r*ratio) / den;`
95	`}`	95	`}`
96	`}`	96	`}`
97		97
98	`static void complex_pow(complex r, complex a, complex *b)`	98	`static void complex_pow(Rcomplex r, Rcomplex a, Rcomplex *b)`
99	`{`	99	`{`
100	`double logr, logi, x, y;`	100	`double logr, logi, x, y;`
101	`int ib;`	101	`int ib;`
102		102
103	`if(b->i == 0.) { /* be fast (and more accurate)*/`	103	`if(b->i == 0.) { /* be fast (and more accurate)*/`
Line 129...	Line 129...
129	`/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */`	129	`/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */`
130		130
131	`SEXP complex_binary(int code, SEXP s1, SEXP s2)`	131	`SEXP complex_binary(int code, SEXP s1, SEXP s2)`
132	`{`	132	`{`
133	`int i, n, n1, n2;`	133	`int i, n, n1, n2;`
134	`complex x1, x2;`	134	`Rcomplex x1, x2;`
135	`SEXP ans;`	135	`SEXP ans;`
136	`n1 = LENGTH(s1);`	136	`n1 = LENGTH(s1);`
137	`n2 = LENGTH(s2);`	137	`n2 = LENGTH(s2);`
138	`n = (n1 > n2) ? n1 : n2;`	138	`n = (n1 > n2) ? n1 : n2;`
139	`/* Note: "s1" and "s1" are protected in the calling code. */`	139	`/* Note: "s1" and "s1" are protected in the calling code. */`
Line 364...	Line 364...
364	`OBJECT(y) = OBJECT(x);`	364	`OBJECT(y) = OBJECT(x);`
365	`UNPROTECT(2);`	365	`UNPROTECT(2);`
366	`return y;`	366	`return y;`
367	`}`	367	`}`
368		368
369	`static void z_rround(complex r, complex x, complex *p)`	369	`static void z_rround(Rcomplex r, Rcomplex x, Rcomplex *p)`
370	`{`	370	`{`
371	`r->r = rround(x->r, p->r);`	371	`r->r = rround(x->r, p->r);`
372	`r->i = rround(x->i, p->r);`	372	`r->i = rround(x->i, p->r);`
373	`}`	373	`}`
374		374
375	`/* Question: This treats real and imaginary parts separately. Should`	375	`/* Question: This treats real and imaginary parts separately. Should`
376	`it do them jointly? */`	376	`it do them jointly? */`
377		377
378	`static void z_prec(complex r, complex x, complex *p)`	378	`static void z_prec(Rcomplex r, Rcomplex x, Rcomplex *p)`
379	`{`	379	`{`
380	`r->r = prec(x->r, p->r);`	380	`r->r = prec(x->r, p->r);`
381	`r->i = prec(x->i, p->r);`	381	`r->i = prec(x->i, p->r);`
382	`}`	382	`}`
383		383
384	`static void z_log(complex r, complex z)`	384	`static void z_log(Rcomplex r, Rcomplex z)`
385	`{`	385	`{`
386	`r->i = atan2(z->i, z->r);`	386	`r->i = atan2(z->i, z->r);`
387	`r->r = log(hypot( z->r, z->i ));`	387	`r->r = log(hypot( z->r, z->i ));`
388	`}`	388	`}`
389		389
390	`static void z_logbase(complex r, complex z, complex *base)`	390	`static void z_logbase(Rcomplex r, Rcomplex z, Rcomplex *base)`
391	`{`	391	`{`
392	`complex t1, t2;`	392	`Rcomplex t1, t2;`
393	`z_log(&t1, z);`	393	`z_log(&t1, z);`
394	`z_log(&t2, base);`	394	`z_log(&t2, base);`
395	`complex_div(r, &t1, &t2);`	395	`complex_div(r, &t1, &t2);`
396	`}`	396	`}`
397		397
398	`static void z_exp(complex r, complex z)`	398	`static void z_exp(Rcomplex r, Rcomplex z)`
399	`{`	399	`{`
400	`double expx;`	400	`double expx;`
401	`expx = exp(z->r);`	401	`expx = exp(z->r);`
402	`r->r = expx * cos(z->i);`	402	`r->r = expx * cos(z->i);`
403	`r->i = expx * sin(z->i);`	403	`r->i = expx * sin(z->i);`
404	`}`	404	`}`
405		405
406	`static void z_sqrt(complex r, complex z)`	406	`static void z_sqrt(Rcomplex r, Rcomplex z)`
407	`{`	407	`{`
408	`double mag;`	408	`double mag;`
409		409
410	`if( (mag = hypot(z->r, z->i)) == 0.0)`	410	`if( (mag = hypot(z->r, z->i)) == 0.0)`
411	`r->r = r->i = 0.0;`	411	`r->r = r->i = 0.0;`
Line 419...	Line 419...
419	`r->i = - r->i;`	419	`r->i = - r->i;`
420	`r->r = z->i / r->i / 2;`	420	`r->r = z->i / r->i / 2;`
421	`}`	421	`}`
422	`}`	422	`}`
423		423
424	`static void z_cos(complex r, complex z)`	424	`static void z_cos(Rcomplex r, Rcomplex z)`
425	`{`	425	`{`
426	`r->r = cos(z->r) * cosh(z->i);`	426	`r->r = cos(z->r) * cosh(z->i);`
427	`r->i = - sin(z->r) * sinh(z->i);`	427	`r->i = - sin(z->r) * sinh(z->i);`
428	`}`	428	`}`
429		429
430	`static void z_sin(complex r, complex z)`	430	`static void z_sin(Rcomplex r, Rcomplex z)`
431	`{`	431	`{`
432	`r->r = sin(z->r) * cosh(z->i);`	432	`r->r = sin(z->r) * cosh(z->i);`
433	`r->i = cos(z->r) * sinh(z->i);`	433	`r->i = cos(z->r) * sinh(z->i);`
434	`}`	434	`}`
435		435
436	`static void z_tan(complex r, complex z)`	436	`static void z_tan(Rcomplex r, Rcomplex z)`
437	`{`	437	`{`
438	`double x2, y2, den;`	438	`double x2, y2, den;`
439	`x2 = 2.0 * z->r;`	439	`x2 = 2.0 * z->r;`
440	`y2 = 2.0 * z->i;`	440	`y2 = 2.0 * z->i;`
441	`den = cos(x2) + cosh(y2);`	441	`den = cos(x2) + cosh(y2);`
Line 444...	Line 444...
444	`}`	444	`}`
445		445
446	`/* Complex Arcsin and Arccos Functions */`	446	`/* Complex Arcsin and Arccos Functions */`
447	`/* Equation (4.4.37) Abramowitz and Stegun */`	447	`/* Equation (4.4.37) Abramowitz and Stegun */`
448		448
449	`static void z_asin(complex r, complex z)`	449	`static void z_asin(Rcomplex r, Rcomplex z)`
450	`{`	450	`{`
451	`double alpha, beta, t1, t2, x, y;`	451	`double alpha, beta, t1, t2, x, y;`
452	`x = z->r;`	452	`x = z->r;`
453	`y = z->i;`	453	`y = z->i;`
454	`t1 = 0.5 * sqrt((x + 1) * (x + 1) + y * y);`	454	`t1 = 0.5 * sqrt((x + 1) * (x + 1) + y * y);`
Line 457...	Line 457...
457	`beta = t1 - t2;`	457	`beta = t1 - t2;`
458	`r->r = asin(beta);`	458	`r->r = asin(beta);`
459	`r->i = log(alpha + sqrt(alpha*alpha - 1));`	459	`r->i = log(alpha + sqrt(alpha*alpha - 1));`
460	`}`	460	`}`
461		461
462	`static void z_acos(complex r, complex z)`	462	`static void z_acos(Rcomplex r, Rcomplex z)`
463	`{`	463	`{`
464	`complex asin;`	464	`Rcomplex asin;`
465	`z_asin(&asin, z);`	465	`z_asin(&asin, z);`
466	`r->r = M_PI_half - asin.r;`	466	`r->r = M_PI_half - asin.r;`
467	`r->i = - asin.i;`	467	`r->i = - asin.i;`
468	`}`	468	`}`
469		469
470	`/* Complex Arctangent Function */`	470	`/* Complex Arctangent Function */`
471	`/* Equation (4.4.39) Abramowitz and Stegun */`	471	`/* Equation (4.4.39) Abramowitz and Stegun */`
472		472
473	`static void z_atan(complex r, complex z)`	473	`static void z_atan(Rcomplex r, Rcomplex z)`
474	`{`	474	`{`
475	`double x, y;`	475	`double x, y;`
476	`x = z->r;`	476	`x = z->r;`
477	`y = z->i;`	477	`y = z->i;`
478	`r->r = 0.5 * atan(2 * x / ( 1 - x * x - y * y));`	478	`r->r = 0.5 * atan(2 * x / ( 1 - x * x - y * y));`
479	`r->i = 0.25 * log((x * x + (y + 1) * (y + 1)) /`	479	`r->i = 0.25 * log((x * x + (y + 1) * (y + 1)) /`
480	`(x * x + (y - 1) * (y - 1)));`	480	`(x * x + (y - 1) * (y - 1)));`
481	`}`	481	`}`
482		482
483	`static void z_atan2(complex r, complex csn, complex *ccs)`	483	`static void z_atan2(Rcomplex r, Rcomplex csn, Rcomplex *ccs)`
484	`{`	484	`{`
485	`complex tmp;`	485	`Rcomplex tmp;`
486	`if (ccs->r == 0 && ccs->i == 0) {`	486	`if (ccs->r == 0 && ccs->i == 0) {`
487	`if(csn->r == 0 && csn->r == 0) {`	487	`if(csn->r == 0 && csn->r == 0) {`
488	`r->r = NA_REAL;`	488	`r->r = NA_REAL;`
489	`r->i = NA_REAL;`	489	`r->i = NA_REAL;`
490	`}`	490	`}`
Line 499...	Line 499...
499	`if(ccs->r < 0) r->r += M_PI;`	499	`if(ccs->r < 0) r->r += M_PI;`
500	`if(r->r > M_PI) r->r -= 2 * M_PI;`	500	`if(r->r > M_PI) r->r -= 2 * M_PI;`
501	`}`	501	`}`
502	`}`	502	`}`
503		503
504	`static void z_acosh(complex r, complex z)`	504	`static void z_acosh(Rcomplex r, Rcomplex z)`
505	`{`	505	`{`
506	`complex a;`	506	`Rcomplex a;`
507	`z_acos(&a, z);`	507	`z_acos(&a, z);`
508	`r->r = -a.i;`	508	`r->r = -a.i;`
509	`r->i = a.r;`	509	`r->i = a.r;`
510	`}`	510	`}`
511		511
512	`static void z_asinh(complex r, complex z)`	512	`static void z_asinh(Rcomplex r, Rcomplex z)`
513	`{`	513	`{`
514	`complex a, b;`	514	`Rcomplex a, b;`
515	`b.r = -z->i;`	515	`b.r = -z->i;`
516	`b.i = z->r;`	516	`b.i = z->r;`
517	`z_asin(&a, &b);`	517	`z_asin(&a, &b);`
518	`r->r = a.i;`	518	`r->r = a.i;`
519	`r->i = -a.r;`	519	`r->i = -a.r;`
520	`}`	520	`}`
521		521
522	`static void z_atanh(complex r, complex z)`	522	`static void z_atanh(Rcomplex r, Rcomplex z)`
523	`{`	523	`{`
524	`complex a, b;`	524	`Rcomplex a, b;`
525	`b.r = -z->i;`	525	`b.r = -z->i;`
526	`b.i = z->r;`	526	`b.i = z->r;`
527	`z_atan(&a, &b);`	527	`z_atan(&a, &b);`
528	`r->r = a.i;`	528	`r->r = a.i;`
529	`r->i = -a.r;`	529	`r->i = -a.r;`
530	`}`	530	`}`
531		531
532	`static void z_cosh(complex r, complex z)`	532	`static void z_cosh(Rcomplex r, Rcomplex z)`
533	`{`	533	`{`
534	`complex a;`	534	`Rcomplex a;`
535	`a.r = -z->i;`	535	`a.r = -z->i;`
536	`a.i = z->r;`	536	`a.i = z->r;`
537	`z_cos(r, &a);`	537	`z_cos(r, &a);`
538	`}`	538	`}`
539		539
540	`static void z_sinh(complex r, complex z)`	540	`static void z_sinh(Rcomplex r, Rcomplex z)`
541	`{`	541	`{`
542	`complex a, b;`	542	`Rcomplex a, b;`
543	`b.r = -z->i;`	543	`b.r = -z->i;`
544	`b.i = z->r;`	544	`b.i = z->r;`
545	`z_sin(&a, &b);`	545	`z_sin(&a, &b);`
546	`r->r = a.i;`	546	`r->r = a.i;`
547	`r->i = -a.r;`	547	`r->i = -a.r;`
548	`}`	548	`}`
549		549
550	`static void z_tanh(complex r, complex z)`	550	`static void z_tanh(Rcomplex r, Rcomplex z)`
551	`{`	551	`{`
552	`complex a, b;`	552	`Rcomplex a, b;`
553	`b.r = -z->i;`	553	`b.r = -z->i;`
554	`b.i = z->r;`	554	`b.i = z->r;`
555	`z_tan(&a, &b);`	555	`z_tan(&a, &b);`
556	`r->r = a.i;`	556	`r->r = a.i;`
557	`r->i = -a.r;`	557	`r->i = -a.r;`
558	`}`	558	`}`
559		559
560	`static void cmath1(void (f)(), complex x, complex *y, int n)`	560	`static void cmath1(void (f)(), Rcomplex x, Rcomplex *y, int n)`
561	`{`	561	`{`
562	`int i;`	562	`int i;`
563	`for (i = 0 ; i < n ; i++) {`	563	`for (i = 0 ; i < n ; i++) {`
564	`if (ISNA(x[i].r) \|\| ISNA(x[i].i)) {`	564	`if (ISNA(x[i].r) \|\| ISNA(x[i].i)) {`
565	`y[i].r = NA_REAL;`	565	`y[i].r = NA_REAL;`
Line 627...	Line 627...
627	`/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */`	627	`/* FIXME : Use the trick in arithmetic.c to eliminate "modulo" ops */`
628		628
629	`static SEXP cmath2(SEXP op, SEXP sa, SEXP sb, void (*f)())`	629	`static SEXP cmath2(SEXP op, SEXP sa, SEXP sb, void (*f)())`
630	`{`	630	`{`
631	`int i, n, na, nb;`	631	`int i, n, na, nb;`
632	`complex ai, bi, a, b, *y;`	632	`Rcomplex ai, bi, a, b, *y;`
633	`SEXP sy;`	633	`SEXP sy;`
634		634
635	`na = length(sa);`	635	`na = length(sa);`
636	`nb = length(sb);`	636	`nb = length(sb);`
637	`if ((na == 0) \|\| (nb == 0))`	637	`if ((na == 0) \|\| (nb == 0))`

The R Project SVN R

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