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

Rev	Author	Line No.	Line
73733	luke	1	`/*`
831	maechler	2	`* R : A Computer Language for Statistical Data Analysis`
77395	maechler	3	`* Copyright (C) 2000-2019 The R Core Team`
		4	`* Copyright (C) 2005 The R Foundation`
		5	`* Copyright (C) 1995-1997 Robert Gentleman and Ross Ihaka`
2	r	6	`*`
		7	`* This program is free software; you can redistribute it and/or modify`
		8	`* it under the terms of the GNU General Public License as published by`
		9	`* the Free Software Foundation; either version 2 of the License, or`
		10	`* (at your option) any later version.`
		11	`*`
		12	`* This program is distributed in the hope that it will be useful,`
		13	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
		14	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
		15	`* GNU General Public License for more details.`
		16	`*`
		17	`* You should have received a copy of the GNU General Public License`
42307	ripley	18	`* along with this program; if not, a copy is available at`
68947	ripley	19	`* https://www.R-project.org/Licenses/`
2	r	20	`*/`
		21
5187	hornik	22	`#ifdef HAVE_CONFIG_H`
7701	hornik	23	`#include <config.h>`
5187	hornik	24	`#endif`
		25
77395	maechler	26	`/* Note: gcc -pedantic may warn in several places about C99 features`
71380	ripley	27	`as extensions.`
		28	`This was a very-long-standing GCC bug, http://gcc.gnu.org/PR7263`
54279	ripley	29	`The system <complex.h> header can work around it: some do.`
71380	ripley	30	`It should have been resolved (after a decade) in 2012.`
54277	ripley	31	`*/`
54275	ripley	32
72661	kalibera	33	`#if defined(HAVE_CTANH) && !defined(HAVE_WORKING_CTANH)`
		34	`#undef HAVE_CTANH`
		35	`#endif`
		36
54256	ripley	37	`#if 0`
54775	ripley	38	`/* For testing substitute fns */`
54293	ripley	39	`#undef HAVE_CARG`
		40	`#undef HAVE_CABS`
54328	ripley	41	`#undef HAVE_CPOW`
54256	ripley	42	`#undef HAVE_CEXP`
		43	`#undef HAVE_CLOG`
		44	`#undef HAVE_CSQRT`
		45	`#undef HAVE_CSIN`
		46	`#undef HAVE_CCOS`
		47	`#undef HAVE_CTAN`
		48	`#undef HAVE_CASIN`
		49	`#undef HAVE_CACOS`
		50	`#undef HAVE_CATAN`
		51	`#undef HAVE_CSINH`
		52	`#undef HAVE_CCOSH`
		53	`#undef HAVE_CTANH`
		54	`#endif`
		55
58276	ripley	56	`#ifdef __SUNPRO_C`
58274	ripley	57	`/* segfaults in Solaris Studio 12.3 */`
		58	`#undef HAVE_CPOW`
		59	`#endif`
		60
15508	pd	61	`#include <Defn.h> /* -> ../include/R_ext/Complex.h */`
60667	ripley	62	`#include <Internal.h>`
11499	ripley	63	`#include <Rmath.h>`
2	r	64
9417	ripley	65	`#include "arithmetic.h" /* complex_* */`
54170	ripley	66	`#include <complex.h>`
73732	luke	67	`#include "Rcomplex.h" /* I, SET_C99_COMPLEX, toC99 */`
68663	luke	68	`#include <R_ext/Itermacros.h>`
2632	maechler	69
68663	luke	70
60166	ripley	71	`/* interval at which to check interrupts, a guess */`
		72	`#define NINTERRUPT 10000000`
		73
		74
41713	ripley	75	`SEXP attribute_hidden complex_unary(ARITHOP_TYPE code, SEXP s1, SEXP call)`
2	r	76	`{`
59048	ripley	77	`R_xlen_t i, n;`
580	ihaka	78	`SEXP ans;`
2	r	79
580	ihaka	80	`switch(code) {`
		81	`case PLUSOP:`
		82	`return s1;`
		83	`case MINUSOP:`
65020	luke	84	`ans = NO_REFERENCES(s1) ? s1 : duplicate(s1);`
73732	luke	85	`Rcomplex *pans = COMPLEX(ans);`
		86	`const Rcomplex *ps1 = COMPLEX_RO(s1);`
59048	ripley	87	`n = XLENGTH(s1);`
580	ihaka	88	`for (i = 0; i < n; i++) {`
73732	luke	89	`Rcomplex x = ps1[i];`
		90	`pans[i].r = -x.r;`
		91	`pans[i].i = -x.i;`
2	r	92	`}`
580	ihaka	93	`return ans;`
		94	`default:`
41713	ripley	95	`errorcall(call, _("invalid complex unary operator"));`
580	ihaka	96	`}`
41713	ripley	97	`return R_NilValue; /* -Wall */`
2	r	98	`}`
		99
54329	ripley	100	`static R_INLINE double complex R_cpow_n(double complex X, int k)`
54256	ripley	101	`{`
54163	ripley	102	`if(k == 0) return (double complex) 1.;`
		103	`else if(k == 1) return X;`
		104	`else if(k < 0) return 1. / R_cpow_n(X, -k);`
49062	maechler	105	`else {/* k > 0 */`
54163	ripley	106	`double complex z = (double complex) 1.;;`
49062	maechler	107	`while (k > 0) {`
54163	ripley	108	`if (k & 1) z = z * X;`
54256	ripley	109	`if (k == 1) break;`
49062	maechler	110	`k >>= 1; /* efficient division by 2; now have k >= 1 */`
		111	`X = X * X;`
		112	`}`
		113	`return z;`
		114	`}`
		115	`}`
		116
55513	ripley	117	`#if defined(Win32)`
55510	ripley	118	`# undef HAVE_CPOW`
		119	`#endif`
49062	maechler	120	`/* reason for this:`
55512	ripley	121	`1) X^n (e.g. for n = +/- 2, 3) is unnecessarily inaccurate in glibc;`
55510	ripley	122	`cut-off 65536 : guided from empirical speed measurements`
54329	ripley	123
55512	ripley	124	`2) On Mingw (but not Mingw-w64) the system cpow is explicitly linked`
		125	`against the (slow) MSVCRT pow, and gets (0+0i)^Y as 0+0i for all Y.`
55611	ripley	126
70779	ripley	127	`3) PPC macOS crashed on powers of 0+0i (at least under Rosetta).`
55618	ripley	128	`Really 0i^-1 should by Inf+NaNi, but getting that portably seems too hard.`
57740	ripley	129	`(C1x's CMPLX will eventually be possible.)`
55510	ripley	130	`*/`
49062	maechler	131
38705	ripley	132	`static double complex mycpow (double complex X, double complex Y)`
35225	ripley	133	`{`
55611	ripley	134	`double complex Z;`
68923	ripley	135	`double yr = creal(Y), yi = cimag(Y);`
55611	ripley	136	`int k;`
		137	`if (X == 0.0) {`
55618	ripley	138	`if (yi == 0.0) Z = R_pow(0.0, yr); else Z = R_NaN + R_NaN*I;`
55611	ripley	139	`} else if (yi == 0.0 && yr == (k = (int) yr) && abs(k) <= 65536)`
54328	ripley	140	`Z = R_cpow_n(X, k);`
55512	ripley	141	`else`
55510	ripley	142	`#ifdef HAVE_CPOW`
54328	ripley	143	`Z = cpow(X, Y);`
		144	`#else`
		145	`{`
55512	ripley	146	`/* Used for FreeBSD and MingGW, hence mainly with gcc */`
54328	ripley	147	`double rho, r, i, theta;`
		148	`r = hypot(creal(X), cimag(X));`
		149	`i = atan2(cimag(X), creal(X));`
		150	`theta = i * yr;`
		151	`if (yi == 0.0)`
		152	`rho = pow(r, yr);`
		153	`else {`
		154	`/* rearrangement of cexp(X * clog(Y)) */`
54329	ripley	155	`r = log(r);`
54328	ripley	156	`theta += r * yi;`
		157	`rho = exp(r * yr - i * yi);`
49062	maechler	158	`}`
55512	ripley	159	`#ifdef __GNUC__`
55611	ripley	160	`__real__ Z = rho * cos(theta);`
		161	`__imag__ Z = rho * sin(theta);`
55512	ripley	162	`#else`
55611	ripley	163	`Z = rho * cos(theta) + (rho * sin(theta)) * I;`
55512	ripley	164	`#endif`
54328	ripley	165	`}`
		166	`#endif`
		167	`return Z;`
35225	ripley	168	`}`
35227	ripley	169
1839	ihaka	170
68663	luke	171
38505	ripley	172	`SEXP attribute_hidden complex_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2)`
2	r	173	`{`
68663	luke	174	`R_xlen_t i, i1, i2, n, n1, n2;`
580	ihaka	175	`SEXP ans;`
16613	ripley	176
54256	ripley	177	`/* Note: "s1" and "s2" are protected in the calling code. */`
59048	ripley	178	`n1 = XLENGTH(s1);`
		179	`n2 = XLENGTH(s2);`
16613	ripley	180	`/* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */`
35069	maechler	181	`if (n1 == 0 \|\| n2 == 0) return(allocVector(CPLXSXP, 0));`
16613	ripley	182
580	ihaka	183	`n = (n1 > n2) ? n1 : n2;`
63462	luke	184	`ans = R_allocOrReuseVector(s1, s2, CPLXSXP, n);`
63461	luke	185	`PROTECT(ans);`
35069	maechler	186
73732	luke	187	`Rcomplex *pans = COMPLEX(ans);`
		188	`const Rcomplex *ps1 = COMPLEX_RO(s1);`
		189	`const Rcomplex *ps2 = COMPLEX_RO(s2);`
		190
580	ihaka	191	`switch (code) {`
		192	`case PLUSOP:`
68663	luke	193	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
73732	luke	194	`Rcomplex x1 = ps1[i1];`
		195	`Rcomplex x2 = ps2[i2];`
		196	`pans[i].r = x1.r + x2.r;`
		197	`pans[i].i = x1.i + x2.i;`
68663	luke	198	`});`
580	ihaka	199	`break;`
		200	`case MINUSOP:`
68663	luke	201	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
73732	luke	202	`Rcomplex x1 = ps1[i1];`
		203	`Rcomplex x2 = ps2[i2];`
		204	`pans[i].r = x1.r - x2.r;`
		205	`pans[i].i = x1.i - x2.i;`
68663	luke	206	`});`
580	ihaka	207	`break;`
		208	`case TIMESOP:`
68663	luke	209	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
73732	luke	210	`SET_C99_COMPLEX(pans, i,`
		211	`toC99(&ps1[i1]) * toC99(&ps2[i2]));`
68663	luke	212	`});`
580	ihaka	213	`break;`
		214	`case DIVOP:`
68663	luke	215	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
73732	luke	216	`SET_C99_COMPLEX(pans, i,`
		217	`toC99(&ps1[i1]) / toC99(&ps2[i2]));`
68663	luke	218	`});`
580	ihaka	219	`break;`
		220	`case POWOP:`
68663	luke	221	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
73732	luke	222	`SET_C99_COMPLEX(pans, i,`
		223	`mycpow(toC99(&ps1[i1]), toC99(&ps2[i2])));`
68663	luke	224	`});`
580	ihaka	225	`break;`
		226	`default:`
32861	ripley	227	`error(_("unimplemented complex operation"));`
580	ihaka	228	`}`
63461	luke	229	`UNPROTECT(1);`
16613	ripley	230
24539	luke	231	`/* quick return if there are no attributes */`
		232	`if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue)`
		233	`return ans;`
		234
16613	ripley	235	`/* Copy attributes from longer argument. */`
63461	luke	236
63462	luke	237	`if (ans != s2 && n == n2 && ATTRIB(s2) != R_NilValue)`
68923	ripley	238	`copyMostAttrib(s2, ans);`
63462	luke	239	`if (ans != s1 && n == n1 && ATTRIB(s1) != R_NilValue)`
68923	ripley	240	`copyMostAttrib(s1, ans); /* Done 2nd so s1's attrs overwrite s2's */`
63461	luke	241
580	ihaka	242	`return ans;`
2	r	243	`}`
		244
36990	ripley	245	`SEXP attribute_hidden do_cmathfuns(SEXP call, SEXP op, SEXP args, SEXP env)`
831	maechler	246	`{`
6098	pd	247	`SEXP x, y = R_NilValue; /* -Wall*/`
59048	ripley	248	`R_xlen_t i, n;`
2	r	249
580	ihaka	250	`checkArity(op, args);`
51316	ripley	251	`check1arg(args, call, "z");`
26610	ripley	252	`if (DispatchGroup("Complex", call, op, args, env, &x))`
45446	ripley	253	`return x;`
580	ihaka	254	`x = CAR(args);`
		255	`if (isComplex(x)) {`
63461	luke	256	`n = XLENGTH(x);`
73732	luke	257	`const Rcomplex *px = COMPLEX_RO(x);`
580	ihaka	258	`switch(PRIMVAL(op)) {`
		259	`case 1: /* Re */`
73732	luke	260	`{`
		261	`y = allocVector(REALSXP, n);`
		262	`double *py = REAL(y);`
		263	`for(i = 0 ; i < n ; i++)`
		264	`py[i] = px[i].r;`
		265	`}`
580	ihaka	266	`break;`
		267	`case 2: /* Im */`
73732	luke	268	`{`
		269	`y = allocVector(REALSXP, n);`
		270	`double *py = REAL(y);`
		271	`for(i = 0 ; i < n ; i++)`
		272	`py[i] = px[i].i;`
		273	`}`
580	ihaka	274	`break;`
		275	`case 3: /* Mod */`
6098	pd	276	`case 6: /* abs */`
73732	luke	277	`{`
		278	`y = allocVector(REALSXP, n);`
		279	`double *py = REAL(y);`
		280	`for(i = 0 ; i < n ; i++)`
54249	ripley	281	`#if HAVE_CABS`
73732	luke	282	`py[i] = cabs(toC99(&px[i]));`
54249	ripley	283	`#else`
73732	luke	284	`py[i] = hypot(px[i].r, px[i].i);`
54249	ripley	285	`#endif`
73732	luke	286	`}`
580	ihaka	287	`break;`
		288	`case 4: /* Arg */`
73732	luke	289	`{`
		290	`y = allocVector(REALSXP, n);`
		291	`double *py = REAL(y);`
		292	`for(i = 0 ; i < n ; i++)`
54249	ripley	293	`#if HAVE_CARG`
73732	luke	294	`py[i] = carg(toC99(&px[i]));`
54249	ripley	295	`#else`
73732	luke	296	`py[i] = atan2(px[i].i, px[i].r);`
54249	ripley	297	`#endif`
73732	luke	298	`}`
580	ihaka	299	`break;`
		300	`case 5: /* Conj */`
73732	luke	301	`{`
		302	`y = NO_REFERENCES(x) ? x : allocVector(CPLXSXP, n);`
		303	`Rcomplex *py = COMPLEX(y);`
		304	`for(i = 0 ; i < n ; i++) {`
		305	`py[i].r = px[i].r;`
		306	`py[i].i = -px[i].i;`
		307	`}`
580	ihaka	308	`}`
		309	`break;`
2	r	310	`}`
580	ihaka	311	`}`
35225	ripley	312	`else if(isNumeric(x)) { /* so no complex numbers involved */`
63461	luke	313	`n = XLENGTH(x);`
580	ihaka	314	`if(isReal(x)) PROTECT(x);`
		315	`else PROTECT(x = coerceVector(x, REALSXP));`
68923	ripley	316	`y = NO_REFERENCES(x) ? x : allocVector(REALSXP, n);`
73732	luke	317	`double *py = REAL(y);`
		318	`const double *px = REAL_RO(x);`
63461	luke	319
580	ihaka	320	`switch(PRIMVAL(op)) {`
		321	`case 1: /* Re */`
		322	`case 5: /* Conj */`
35225	ripley	323	`for(i = 0 ; i < n ; i++)`
73732	luke	324	`py[i] = px[i];`
580	ihaka	325	`break;`
		326	`case 2: /* Im */`
36129	ripley	327	`for(i = 0 ; i < n ; i++)`
73732	luke	328	`py[i] = 0.0;`
36129	ripley	329	`break;`
580	ihaka	330	`case 4: /* Arg */`
35225	ripley	331	`for(i = 0 ; i < n ; i++)`
73732	luke	332	`if(ISNAN(px[i]))`
		333	`py[i] = px[i];`
		334	`else if (px[i] >= 0)`
		335	`py[i] = 0;`
580	ihaka	336	`else`
73732	luke	337	`py[i] = M_PI;`
580	ihaka	338	`break;`
		339	`case 3: /* Mod */`
6098	pd	340	`case 6: /* abs */`
35225	ripley	341	`for(i = 0 ; i < n ; i++)`
73732	luke	342	`py[i] = fabs(px[i]);`
580	ihaka	343	`break;`
		344	`}`
		345	`UNPROTECT(1);`
		346	`}`
32861	ripley	347	`else errorcall(call, _("non-numeric argument to function"));`
63461	luke	348
63462	luke	349	`if (x != y && ATTRIB(x) != R_NilValue) {`
68923	ripley	350	`PROTECT(x);`
		351	`PROTECT(y);`
69108	luke	352	`SHALLOW_DUPLICATE_ATTRIB(y, x);`
68923	ripley	353	`UNPROTECT(2);`
63461	luke	354	`}`
580	ihaka	355	`return y;`
2	r	356	`}`
		357
77395	maechler	358	`/* Implementing signif(<complex>) and used in format.c and printutils.c */`
35250	ripley	359	`#define MAX_DIGITS 22`
73727	luke	360	`void attribute_hidden z_prec_r(Rcomplex r, const Rcomplex x, double digits)`
2	r	361	`{`
77395	maechler	362	`// Implement r <- signif(x, digits)`
35250	ripley	363
		364	`r->r = x->r; r->i = x->i;`
77395	maechler	365	`double m = 0.0,`
		366	`m1 = fabs(x->r),`
		367	`m2 = fabs(x->i);`
35253	ripley	368	`if(R_FINITE(m1)) m = m1;`
		369	`if(R_FINITE(m2) && m2 > m) m = m2;`
35250	ripley	370	`if (m == 0.0) return;`
		371	`if (!R_FINITE(digits)) {`
		372	`if(digits > 0) return; else {r->r = r->i = 0.0; return ;}`
		373	`}`
77395	maechler	374	`int dig = (int)floor(digits+0.5);`
35250	ripley	375	`if (dig > MAX_DIGITS) return; else if (dig < 1) dig = 1;`
77395	maechler	376	`int mag = (int)floor(log10(m));`
35250	ripley	377	`dig = dig - mag - 1;`
		378	`if (dig > 306) {`
		379	`double pow10 = 1.0e4;`
		380	`digits = (double)(dig - 4);`
54776	ripley	381	`r->r = fround(pow10 * x->r, digits)/pow10;`
		382	`r->i = fround(pow10 * x->i, digits)/pow10;`
35250	ripley	383	`} else {`
		384	`digits = (double)(dig);`
54776	ripley	385	`r->r = fround(x->r, digits);`
		386	`r->i = fround(x->i, digits);`
35250	ripley	387	`}`
2	r	388	`}`
		389
54275	ripley	390	`/* These substitute functions are rarely used, and not extensively`
		391	`tested, e.g. over CRAN. Please do not change without very good`
		392	`reason!`
54253	ripley	393
54275	ripley	394	`Currently (Feb 2011) they are used on FreeBSD.`
		395	`*/`
35225	ripley	396
54249	ripley	397	`#ifndef HAVE_CLOG`
54256	ripley	398	`#define clog R_clog`
54275	ripley	399	`/* FIXME: maybe add full IEC60559 support */`
54253	ripley	400	`static double complex clog(double complex x)`
54249	ripley	401	`{`
		402	`double xr = creal(x), xi = cimag(x);`
		403	`return log(hypot(xr, xi)) + atan2(xi, xr)*I;`
		404	`}`
		405	`#endif`
54258	ripley	406
54249	ripley	407	`#ifndef HAVE_CSQRT`
54256	ripley	408	`#define csqrt R_csqrt`
54249	ripley	409	`/* FreeBSD does have this one */`
54253	ripley	410	`static double complex csqrt(double complex x)`
54249	ripley	411	`{`
		412	`return mycpow(x, 0.5+0.0*I);`
		413	`}`
		414	`#endif`
54258	ripley	415
54249	ripley	416	`#ifndef HAVE_CEXP`
54256	ripley	417	`#define cexp R_cexp`
54253	ripley	418	`/* FIXME: check/add full IEC60559 support */`
		419	`static double complex cexp(double complex x)`
54249	ripley	420	`{`
		421	`double expx = exp(creal(x)), y = cimag(x);`
		422	`return expx * cos(y) + (expx * sin(y)) * I;`
		423	`}`
		424	`#endif`
54258	ripley	425
54249	ripley	426	`#ifndef HAVE_CCOS`
54256	ripley	427	`#define ccos R_ccos`
54253	ripley	428	`static double complex ccos(double complex x)`
54249	ripley	429	`{`
54256	ripley	430	`double xr = creal(x), xi = cimag(x);`
54262	ripley	431	`return cos(xr)cosh(xi) - sin(xr)sinh(xi)I; / A&S 4.3.56 */`
54249	ripley	432	`}`
		433	`#endif`
54258	ripley	434
54249	ripley	435	`#ifndef HAVE_CSIN`
54256	ripley	436	`#define csin R_csin`
54253	ripley	437	`static double complex csin(double complex x)`
54249	ripley	438	`{`
54256	ripley	439	`double xr = creal(x), xi = cimag(x);`
54262	ripley	440	`return sin(xr)cosh(xi) + cos(xr)sinh(xi)I; / A&S 4.3.55 */`
54249	ripley	441	`}`
		442	`#endif`
54258	ripley	443
54275	ripley	444	`#ifndef HAVE_CTAN`
		445	`#define ctan R_ctan`
		446	`static double complex ctan(double complex z)`
		447	`{`
		448	`/* A&S 4.3.57 */`
		449	`double x2, y2, den, ri;`
		450	`x2 = 2.0 * creal(z);`
		451	`y2 = 2.0 * cimag(z);`
		452	`den = cos(x2) + cosh(y2);`
		453	`/* any threshold between -log(DBL_EPSILON) and log(DBL_XMAX) will do*/`
		454	`if (ISNAN(y2) \|\| fabs(y2) < 50.0) ri = sinh(y2)/den;`
		455	`else ri = (y2 < 0 ? -1.0 : 1.0);`
		456	`return sin(x2)/den + ri * I;`
		457	`}`
		458	`#endif`
		459
54249	ripley	460	`#ifndef HAVE_CASIN`
54256	ripley	461	`#define casin R_casin`
		462	`static double complex casin(double complex z)`
54249	ripley	463	`{`
54262	ripley	464	`/* A&S 4.4.37 */`
54258	ripley	465	`double alpha, t1, t2, x = creal(z), y = cimag(z), ri;`
54256	ripley	466	`t1 = 0.5 * hypot(x + 1, y);`
		467	`t2 = 0.5 * hypot(x - 1, y);`
		468	`alpha = t1 + t2;`
		469	`ri = log(alpha + sqrt(alpha*alpha - 1));`
68923	ripley	470	`/* This comes from`
54262	ripley	471	`'z_asin() is continuous from below if x >= 1`
68923	ripley	472	`and continuous from above if x <= -1.'`
54262	ripley	473	`*/`
		474	`if(y < 0 \|\| (y == 0 && x > 1)) ri *= -1;`
54256	ripley	475	`return asin(t1 - t2) + ri*I;`
54249	ripley	476	`}`
		477	`#endif`
54258	ripley	478
54249	ripley	479	`#ifndef HAVE_CACOS`
54256	ripley	480	`#define cacos R_cacos`
		481	`static double complex cacos(double complex z)`
54249	ripley	482	`{`
54256	ripley	483	`return M_PI_2 - casin(z);`
54249	ripley	484	`}`
		485	`#endif`
54258	ripley	486
54249	ripley	487	`#ifndef HAVE_CATAN`
54256	ripley	488	`#define catan R_catan`
		489	`static double complex catan(double complex z)`
54249	ripley	490	`{`
54256	ripley	491	`double x = creal(z), y = cimag(z), rr, ri;`
		492	`rr = 0.5 * atan2(2 * x, (1 - x * x - y * y));`
		493	`ri = 0.25 * log((x * x + (y + 1) * (y + 1)) /`
		494	`(x * x + (y - 1) * (y - 1)));`
54275	ripley	495	`return rr + ri*I;`
54249	ripley	496	`}`
		497	`#endif`
54258	ripley	498
54249	ripley	499	`#ifndef HAVE_CCOSH`
54256	ripley	500	`#define ccosh R_ccosh`
		501	`static double complex ccosh(double complex z)`
54249	ripley	502	`{`
54262	ripley	503	`return ccos(z * I); /* A&S 4.5.8 */`
54249	ripley	504	`}`
		505	`#endif`
54258	ripley	506
54249	ripley	507	`#ifndef HAVE_CSINH`
54256	ripley	508	`#define csinh R_csinh`
		509	`static double complex csinh(double complex z)`
54249	ripley	510	`{`
54262	ripley	511	`return -I * csin(z * I); /* A&S 4.5.7 */`
54249	ripley	512	`}`
		513	`#endif`
54258	ripley	514
54275	ripley	515	`static double complex z_tan(double complex z)`
54249	ripley	516	`{`
54275	ripley	517	`double y = cimag(z);`
		518	`double complex r = ctan(z);`
		519	`if(R_FINITE(y) && fabs(y) > 25.0) {`
		520	`/* at this point the real part is nearly zero, and the`
		521	`imaginary part is one: but some OSes get the imag as NaN */`
		522	`#if __GNUC__`
		523	`__imag__ r = y < 0 ? -1.0 : 1.0;`
		524	`#else`
		525	`r = creal(r) + (y < 0 ? -1.0 : 1.0) * I;`
		526	`#endif`
		527	`}`
54262	ripley	528	`return r;`
54249	ripley	529	`}`
54258	ripley	530
72412	kalibera	531	`#ifndef HAVE_CTANH`
		532	`#define ctanh R_ctanh`
		533	`static double complex ctanh(double complex z)`
		534	`{`
		535	`return -I * z_tan(z * I); /* A&S 4.5.9 */`
		536	`}`
		537	`#endif`
		538
		539
54275	ripley	540	`/* Don't rely on the OS at the branch cuts */`
		541
		542	`static double complex z_asin(double complex z)`
54249	ripley	543	`{`
54275	ripley	544	`if(cimag(z) == 0 && fabs(creal(z)) > 1) {`
		545	`double alpha, t1, t2, x = creal(z), ri;`
		546	`t1 = 0.5 * fabs(x + 1);`
		547	`t2 = 0.5 * fabs(x - 1);`
		548	`alpha = t1 + t2;`
		549	`ri = log(alpha + sqrt(alpha*alpha - 1));`
		550	`if(x > 1) ri *= -1;`
		551	`return asin(t1 - t2) + ri*I;`
		552	`}`
		553	`return casin(z);`
54249	ripley	554	`}`
54258	ripley	555
54275	ripley	556	`static double complex z_acos(double complex z)`
54249	ripley	557	`{`
54275	ripley	558	`if(cimag(z) == 0 && fabs(creal(z)) > 1) return M_PI_2 - z_asin(z);`
		559	`return cacos(z);`
54249	ripley	560	`}`
		561
54275	ripley	562	`static double complex z_atan(double complex z)`
		563	`{`
		564	`if(creal(z) == 0 && fabs(cimag(z)) > 1) {`
		565	`double y = cimag(z), rr, ri;`
		566	`rr = (y > 0) ? M_PI_2 : -M_PI_2;`
		567	`ri = 0.25 * log(((y + 1) * (y + 1))/((y - 1) * (y - 1)));`
		568	`return rr + ri*I;`
		569	`}`
		570	`return catan(z);`
		571	`}`
54273	ripley	572
54275	ripley	573	`static double complex z_acosh(double complex z)`
		574	`{`
		575	`return z_acos(z) * I;`
		576	`}`
		577
		578	`static double complex z_asinh(double complex z)`
		579	`{`
		580	`return -I * z_asin(z * I);`
		581	`}`
		582
		583	`static double complex z_atanh(double complex z)`
		584	`{`
		585	`return -I * z_atan(z * I);`
		586	`}`
		587
54273	ripley	588	`static Rboolean cmath1(double complex (*f)(double complex),`
73732	luke	589	`const Rcomplex x, Rcomplex y, R_xlen_t n)`
54273	ripley	590	`{`
59090	ripley	591	`R_xlen_t i;`
54273	ripley	592	`Rboolean naflag = FALSE;`
		593	`for (i = 0 ; i < n ; i++) {`
		594	`if (ISNA(x[i].r) \|\| ISNA(x[i].i)) {`
		595	`y[i].r = NA_REAL; y[i].i = NA_REAL;`
		596	`} else {`
		597	`SET_C99_COMPLEX(y, i, f(toC99(x + i)));`
55516	ripley	598	`if ( (ISNAN(y[i].r) \|\| ISNAN(y[i].i)) &&`
		599	`!(ISNAN(x[i].r) \|\| ISNAN(x[i].i)) ) naflag = TRUE;`
54273	ripley	600	`}`
		601	`}`
		602	`return naflag;`
		603	`}`
		604
38505	ripley	605	`SEXP attribute_hidden complex_math1(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	606	`{`
580	ihaka	607	`SEXP x, y;`
59048	ripley	608	`R_xlen_t n;`
19912	duncan	609	`Rboolean naflag = FALSE;`
54273	ripley	610
6098	pd	611	`PROTECT(x = CAR(args));`
68586	ripley	612	`n = XLENGTH(x);`
6098	pd	613	`PROTECT(y = allocVector(CPLXSXP, n));`
2	r	614
73732	luke	615	`const Rcomplex *px = COMPLEX_RO(x);`
		616	`Rcomplex *py = COMPLEX(y);`
		617
580	ihaka	618	`switch (PRIMVAL(op)) {`
73732	luke	619	`case 10003: naflag = cmath1(clog, px, py, n); break;`
77395	maechler	620	`// 1: floor`
		621	`// 2: ceil[ing]`
73732	luke	622	`case 3: naflag = cmath1(csqrt, px, py, n); break;`
77395	maechler	623	`// 4: sign`
73732	luke	624	`case 10: naflag = cmath1(cexp, px, py, n); break;`
77395	maechler	625	`// 11: expm1`
		626	`// 12: log1p`
73732	luke	627	`case 20: naflag = cmath1(ccos, px, py, n); break;`
		628	`case 21: naflag = cmath1(csin, px, py, n); break;`
		629	`case 22: naflag = cmath1(z_tan, px, py, n); break;`
		630	`case 23: naflag = cmath1(z_acos, px, py, n); break;`
		631	`case 24: naflag = cmath1(z_asin, px, py, n); break;`
		632	`case 25: naflag = cmath1(z_atan, px, py, n); break;`
77395	maechler	633
73732	luke	634	`case 30: naflag = cmath1(ccosh, px, py, n); break;`
		635	`case 31: naflag = cmath1(csinh, px, py, n); break;`
		636	`case 32: naflag = cmath1(ctanh, px, py, n); break;`
		637	`case 33: naflag = cmath1(z_acosh, px, py, n); break;`
		638	`case 34: naflag = cmath1(z_asinh, px, py, n); break;`
		639	`case 35: naflag = cmath1(z_atanh, px, py, n); break;`
54174	ripley	640
580	ihaka	641	`default:`
41725	ripley	642	`/* such as sign, gamma */`
		643	`errorcall(call, _("unimplemented complex function"));`
580	ihaka	644	`}`
6098	pd	645	`if (naflag)`
55518	ripley	646	`warningcall(call, "NaNs produced in function \"%s\"", PRIMNAME(op));`
69108	luke	647	`SHALLOW_DUPLICATE_ATTRIB(y, x);`
6098	pd	648	`UNPROTECT(2);`
580	ihaka	649	`return y;`
2	r	650	`}`
		651
54262	ripley	652	`static void z_rround(Rcomplex r, Rcomplex x, Rcomplex *p)`
		653	`{`
54776	ripley	654	`r->r = fround(x->r, p->r);`
		655	`r->i = fround(x->i, p->r);`
54262	ripley	656	`}`
		657
54174	ripley	658	`static void z_prec(Rcomplex r, Rcomplex x, Rcomplex *p)`
		659	`{`
		660	`z_prec_r(r, x, p->r);`
		661	`}`
		662
		663	`static void z_logbase(Rcomplex r, Rcomplex z, Rcomplex *base)`
		664	`{`
		665	`double complex dz = toC99(z), dbase = toC99(base);`
		666	`SET_C99_COMPLEX(r, 0, clog(dz)/clog(dbase));`
		667	`}`
		668
		669	`static void z_atan2(Rcomplex r, Rcomplex csn, Rcomplex *ccs)`
		670	`{`
		671	`double complex dr, dcsn = toC99(csn), dccs = toC99(ccs);`
		672	`if (dccs == 0) {`
		673	`if(dcsn == 0) {`
54256	ripley	674	`r->r = NA_REAL; r->i = NA_REAL; /* Why not R_NaN? */`
54174	ripley	675	`return;`
54293	ripley	676	`} else {`
		677	`double y = creal(dcsn);`
		678	`if (ISNAN(y)) dr = y;`
		679	`else dr = ((y >= 0) ? M_PI_2 : -M_PI_2);`
		680	`}`
54174	ripley	681	`} else {`
		682	`dr = catan(dcsn / dccs);`
		683	`if(creal(dccs) < 0) dr += M_PI;`
		684	`if(creal(dr) > M_PI) dr -= 2 * M_PI;`
		685	`}`
		686	`SET_C99_COMPLEX(r, 0, dr);`
		687	`}`
		688
54273	ripley	689
		690	`/* Complex Functions of Two Arguments */`
		691
		692	`typedef void (cm2_fun)(Rcomplex , Rcomplex , Rcomplex );`
		693	`SEXP attribute_hidden complex_math2(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	694	`{`
68663	luke	695	`R_xlen_t i, n, na, nb, ia, ib;`
73732	luke	696	`Rcomplex ai, bi, *y;`
		697	`const Rcomplex a, b;`
54273	ripley	698	`SEXP sa, sb, sy;`
		699	`Rboolean naflag = FALSE;`
		700	`cm2_fun f;`
		701
		702	`switch (PRIMVAL(op)) {`
		703	`case 0: /* atan2 */`
		704	`f = z_atan2; break;`
		705	`case 10001: /* round */`
		706	`f = z_rround; break;`
		707	`case 2: /* passed from do_log1arg */`
		708	`case 10:`
		709	`case 10003: /* passed from do_log */`
		710	`f = z_logbase; break;`
		711	`case 10004: /* signif */`
		712	`f = z_prec; break;`
		713	`default:`
70830	luke	714	`error_return(_("unimplemented complex function"));`
54273	ripley	715	`}`
		716
		717	`PROTECT(sa = coerceVector(CAR(args), CPLXSXP));`
		718	`PROTECT(sb = coerceVector(CADR(args), CPLXSXP));`
59048	ripley	719	`na = XLENGTH(sa); nb = XLENGTH(sb);`
67024	luke	720	`if ((na == 0) \|\| (nb == 0)) {`
68923	ripley	721	`UNPROTECT(2);`
		722	`return(allocVector(CPLXSXP, 0));`
67024	luke	723	`}`
580	ihaka	724	`n = (na < nb) ? nb : na;`
		725	`PROTECT(sy = allocVector(CPLXSXP, n));`
73732	luke	726	`a = COMPLEX_RO(sa); b = COMPLEX_RO(sb);`
		727	`y = COMPLEX(sy);`
68663	luke	728	`MOD_ITERATE2(n, na, nb, i, ia, ib, {`
		729	`ai = a[ia]; bi = b[ib];`
6098	pd	730	`if(ISNA(ai.r) && ISNA(ai.i) &&`
		731	`ISNA(bi.r) && ISNA(bi.i)) {`
54273	ripley	732	`y[i].r = NA_REAL; y[i].i = NA_REAL;`
54170	ripley	733	`} else {`
6098	pd	734	`f(&y[i], &ai, &bi);`
55516	ripley	735	`if ( (ISNAN(y[i].r) \|\| ISNAN(y[i].i)) &&`
		736	`!(ISNAN(ai.r) \|\| ISNAN(ai.i) \|\| ISNAN(bi.r) \|\| ISNAN(bi.i)) )`
		737	`naflag = TRUE;`
2	r	738	`}`
68663	luke	739	`});`
6098	pd	740	`if (naflag)`
70830	luke	741	`warning("NaNs produced in function \"%s\"", PRIMNAME(op));`
580	ihaka	742	`if(n == na) {`
69108	luke	743	`SHALLOW_DUPLICATE_ATTRIB(sy, sa);`
54170	ripley	744	`} else if(n == nb) {`
69108	luke	745	`SHALLOW_DUPLICATE_ATTRIB(sy, sb);`
580	ihaka	746	`}`
		747	`UNPROTECT(3);`
		748	`return sy;`
2	r	749	`}`
		750
36990	ripley	751	`SEXP attribute_hidden do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	752	`{`
1839	ihaka	753	`/* complex(length, real, imaginary) */`
580	ihaka	754	`SEXP ans, re, im;`
59048	ripley	755	`R_xlen_t i, na, nr, ni;`
69326	luke	756
		757	`checkArity(op, args);`
580	ihaka	758	`na = asInteger(CAR(args));`
		759	`if(na == NA_INTEGER \|\| na < 0)`
41686	ripley	760	`error(_("invalid length"));`
580	ihaka	761	`PROTECT(re = coerceVector(CADR(args), REALSXP));`
		762	`PROTECT(im = coerceVector(CADDR(args), REALSXP));`
59048	ripley	763	`nr = XLENGTH(re);`
		764	`ni = XLENGTH(im);`
1839	ihaka	765	`/* is always true: if (na >= 0) {*/`
		766	`na = (nr > na) ? nr : na;`
		767	`na = (ni > na) ? ni : na;`
		768	`/* }*/`
580	ihaka	769	`ans = allocVector(CPLXSXP, na);`
73732	luke	770	`Rcomplex *pans = COMPLEX(ans);`
580	ihaka	771	`for(i=0 ; i<na ; i++) {`
73732	luke	772	`pans[i].r = 0;`
		773	`pans[i].i = 0;`
580	ihaka	774	`}`
		775	`UNPROTECT(2);`
		776	`if(na > 0 && nr > 0) {`
73732	luke	777	`const double *p_re = REAL_RO(re);`
580	ihaka	778	`for(i=0 ; i<na ; i++)`
73732	luke	779	`pans[i].r = p_re[i%nr];`
580	ihaka	780	`}`
		781	`if(na > 0 && ni > 0) {`
73732	luke	782	`const double *p_im = REAL_RO(im);`
580	ihaka	783	`for(i=0 ; i<na ; i++)`
73732	luke	784	`pans[i].i = p_im[i%ni];`
580	ihaka	785	`}`
		786	`return ans;`
2	r	787	`}`
		788
61750	ripley	789	`static void R_cpolyroot(double opr, double opi, int *degree,`
		790	`double zeror, double zeroi, Rboolean *fail);`
580	ihaka	791
36990	ripley	792	`SEXP attribute_hidden do_polyroot(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	793	`{`
580	ihaka	794	`SEXP z, zr, zi, r, rr, ri;`
10866	maechler	795	`Rboolean fail;`
		796	`int degree, i, n;`
		797
580	ihaka	798	`checkArity(op, args);`
		799	`z = CAR(args);`
		800	`switch(TYPEOF(z)) {`
		801	`case CPLXSXP:`
		802	`PROTECT(z);`
		803	`break;`
		804	`case REALSXP:`
		805	`case INTSXP:`
		806	`case LGLSXP:`
		807	`PROTECT(z = coerceVector(z, CPLXSXP));`
		808	`break;`
		809	`default:`
31908	ripley	810	`UNIMPLEMENTED_TYPE("polyroot", z);`
580	ihaka	811	`}`
61750	ripley	812	`#ifdef LONG_VECTOR_SUPPORT`
		813	`R_xlen_t nn = XLENGTH(z);`
		814	`if (nn > R_SHORT_LEN_MAX) error("long vectors are not supported");`
		815	`n = (int) nn;`
		816	`#else`
		817	`n = LENGTH(z);`
		818	`#endif`
73732	luke	819	`const Rcomplex *pz = COMPLEX_RO(z);`
18490	ripley	820	`degree = 0;`
		821	`for(i = 0; i < n; i++) {`
73732	luke	822	`if(pz[i].r!= 0.0 \|\| pz[i].i != 0.0) degree = i;`
18490	ripley	823	`}`
		824	`n = degree + 1; /* omit trailing zeroes */`
580	ihaka	825	`if(degree >= 1) {`
		826	`PROTECT(rr = allocVector(REALSXP, n));`
		827	`PROTECT(ri = allocVector(REALSXP, n));`
		828	`PROTECT(zr = allocVector(REALSXP, n));`
		829	`PROTECT(zi = allocVector(REALSXP, n));`
2	r	830
73732	luke	831	`double *p_rr = REAL(rr);`
		832	`double *p_ri = REAL(ri);`
		833	`double *p_zr = REAL(zr);`
		834	`double *p_zi = REAL(zi);`
		835
61750	ripley	836	`for(i = 0 ; i < n ; i++) {`
73732	luke	837	`if(!R_FINITE(pz[i].r) \|\| !R_FINITE(pz[i].i))`
41686	ripley	838	`error(_("invalid polynomial coefficient"));`
73732	luke	839	`p_zr[degree-i] = pz[i].r;`
		840	`p_zi[degree-i] = pz[i].i;`
2	r	841	`}`
73732	luke	842	`R_cpolyroot(p_zr, p_zi, &degree, p_rr, p_ri, &fail);`
41686	ripley	843	`if(fail) error(_("root finding code failed"));`
580	ihaka	844	`UNPROTECT(2);`
		845	`r = allocVector(CPLXSXP, degree);`
73732	luke	846	`Rcomplex *pr = COMPLEX(r);`
61750	ripley	847	`for(i = 0 ; i < degree ; i++) {`
73732	luke	848	`pr[i].r = p_rr[i];`
		849	`pr[i].i = p_ri[i];`
2	r	850	`}`
580	ihaka	851	`UNPROTECT(3);`
		852	`}`
		853	`else {`
		854	`UNPROTECT(1);`
		855	`r = allocVector(CPLXSXP, 0);`
		856	`}`
		857	`return r;`
2	r	858	`}`
61750	ripley	859
		860	`/* Formerly src/appl/cpoly.c:`
		861	`*`
		862	`* Copyright (C) 1997-1998 Ross Ihaka`
		863	`* Copyright (C) 1999-2001 R Core Team`
		864	`*`
		865	`* cpoly finds the zeros of a complex polynomial.`
		866	`*`
		867	`* On Entry`
		868	`*`
		869	`* opr, opi - double precision vectors of real and`
		870	`* imaginary parts of the coefficients in`
		871	`* order of decreasing powers.`
		872	`*`
		873	`* degree - int degree of polynomial.`
		874	`*`
		875	`*`
		876	`* On Return`
		877	`*`
		878	`* zeror, zeroi - output double precision vectors of`
		879	`* real and imaginary parts of the zeros.`
		880	`*`
		881	`* fail - output int parameter, true only if`
		882	`* leading coefficient is zero or if cpoly`
		883	`* has found fewer than degree zeros.`
		884	`*`
		885	`* The program has been written to reduce the chance of overflow`
		886	`* occurring. If it does occur, there is still a possibility that`
		887	`* the zerofinder will work provided the overflowed quantity is`
		888	`* replaced by a large number.`
		889	`*`
		890	`* This is a C translation of the following.`
		891	`*`
		892	`* TOMS Algorithm 419`
		893	`* Jenkins and Traub.`
		894	`* Comm. ACM 15 (1972) 97-99.`
		895	`*`
		896	`* Ross Ihaka`
		897	`* February 1997`
		898	`*/`
		899
		900	`#include <R_ext/Arith.h> /* for declaration of hypot */`
		901	`#include <R_ext/Memory.h> /* for declaration of R_alloc */`
		902
		903	`#include <float.h> /* for FLT_RADIX */`
		904
		905	`#include <Rmath.h> /* for R_pow_di */`
		906
		907	`static void calct(Rboolean *);`
		908	`static Rboolean fxshft(int, double , double );`
		909	`static Rboolean vrshft(int, double , double );`
		910	`static void nexth(Rboolean);`
		911	`static void noshft(int);`
		912
		913	`static void polyev(int, double, double,`
		914	`double , double , double , double , double , double );`
		915	`static double errev(int, double , double , double, double, double, double);`
		916	`static double cpoly_cauchy(int, double , double );`
		917	`static double cpoly_scale(int, double *, double, double, double, double);`
		918	`static void cdivid(double, double, double, double, double , double );`
		919
		920	`/* Global Variables (too many!) */`
		921
		922	`static int nn;`
		923	`static double pr, pi, hr, hi, qpr, qpi, qhr, qhi, shr, shi;`
		924	`static double sr, si;`
		925	`static double tr, ti;`
		926	`static double pvr, pvi;`
		927
		928	`static const double eta = DBL_EPSILON;`
		929	`static const double are = /* eta = */DBL_EPSILON;`
		930	`static const double mre = 2. * M_SQRT2 * /* eta, i.e. */DBL_EPSILON;`
		931	`static const double infin = DBL_MAX;`
		932
		933	`static void R_cpolyroot(double opr, double opi, int *degree,`
		934	`double zeror, double zeroi, Rboolean *fail)`
		935	`{`
		936	`static const double smalno = DBL_MIN;`
		937	`static const double base = (double)FLT_RADIX;`
		938	`static int d_n, i, i1, i2;`
		939	`static double zi, zr, xx, yy;`
		940	`static double bnd, xxx;`
		941	`Rboolean conv;`
		942	`int d1;`
		943	`double *tmp;`
		944	`static const double cosr =/* cos 94 */ -0.06975647374412529990;`
		945	`static const double sinr =/* sin 94 */ 0.99756405025982424767;`
		946	`xx = M_SQRT1_2;/* 1/sqrt(2) = 0.707.... */`
		947
		948	`yy = -xx;`
		949	`*fail = FALSE;`
		950
		951	`nn = *degree;`
		952	`d1 = nn - 1;`
		953
		954	`/* algorithm fails if the leading coefficient is zero. */`
		955
		956	`if (opr[0] == 0. && opi[0] == 0.) {`
		957	`*fail = TRUE;`
		958	`return;`
		959	`}`
		960
		961	`/* remove the zeros at the origin if any. */`
		962
		963	`while (opr[nn] == 0. && opi[nn] == 0.) {`
		964	`d_n = d1-nn+1;`
		965	`zeror[d_n] = 0.;`
		966	`zeroi[d_n] = 0.;`
		967	`nn--;`
		968	`}`
		969	`nn++;`
		970	`/-- Now, global var. nn := #{coefficients} = (relevant degree)+1 /`
		971
		972	`if (nn == 1) return;`
		973
		974	`/* Use a single allocation as these as small */`
63181	ripley	975	`const void *vmax = vmaxget();`
61750	ripley	976	`tmp = (double ) R_alloc((size_t) (10nn), sizeof(double));`
		977	`pr = tmp; pi = tmp + nn; hr = tmp + 2nn; hi = tmp + 3nn;`
		978	`qpr = tmp + 4nn; qpi = tmp + 5nn; qhr = tmp + 6nn; qhi = tmp + 7nn;`
		979	`shr = tmp + 8nn; shi = tmp + 9nn;`
68923	ripley	980
61750	ripley	981	`/* make a copy of the coefficients and shr[] = \| p[] \| */`
		982	`for (i = 0; i < nn; i++) {`
		983	`pr[i] = opr[i];`
		984	`pi[i] = opi[i];`
		985	`shr[i] = hypot(pr[i], pi[i]);`
		986	`}`
		987
		988	`/* scale the polynomial with factor 'bnd'. */`
		989	`bnd = cpoly_scale(nn, shr, eta, infin, smalno, base);`
		990	`if (bnd != 1.) {`
		991	`for (i=0; i < nn; i++) {`
		992	`pr[i] *= bnd;`
		993	`pi[i] *= bnd;`
		994	`}`
		995	`}`
		996
		997	`/* start the algorithm for one zero */`
		998
		999	`while (nn > 2) {`
		1000
		1001	`/* calculate bnd, a lower bound on the modulus of the zeros. */`
		1002
		1003	`for (i=0 ; i < nn ; i++)`
		1004	`shr[i] = hypot(pr[i], pi[i]);`
		1005	`bnd = cpoly_cauchy(nn, shr, shi);`
		1006
		1007	`/* outer loop to control 2 major passes */`
		1008	`/* with different sequences of shifts */`
		1009
		1010	`for (i1 = 1; i1 <= 2; i1++) {`
		1011
		1012	`/* first stage calculation, no shift */`
		1013
		1014	`noshft(5);`
		1015
		1016	`/* inner loop to select a shift */`
		1017	`for (i2 = 1; i2 <= 9; i2++) {`
		1018
		1019	`/* shift is chosen with modulus bnd */`
		1020	`/* and amplitude rotated by 94 degrees */`
		1021	`/* from the previous shift */`
		1022
		1023	`xxx= cosr * xx - sinr * yy;`
		1024	`yy = sinr * xx + cosr * yy;`
		1025	`xx = xxx;`
		1026	`sr = bnd * xx;`
		1027	`si = bnd * yy;`
		1028
		1029	`/* second stage calculation, fixed shift */`
		1030
		1031	`conv = fxshft(i2 * 10, &zr, &zi);`
		1032	`if (conv)`
		1033	`goto L10;`
		1034	`}`
		1035	`}`
		1036
		1037	`/* the zerofinder has failed on two major passes */`
		1038	`/* return empty handed */`
		1039
		1040	`*fail = TRUE;`
63171	ripley	1041	`vmaxset(vmax);`
61750	ripley	1042	`return;`
		1043
		1044	`/* the second stage jumps directly to the third stage iteration.`
		1045	`* if successful, the zero is stored and the polynomial deflated.`
		1046	`*/`
		1047	`L10:`
		1048	`d_n = d1+2 - nn;`
		1049	`zeror[d_n] = zr;`
		1050	`zeroi[d_n] = zi;`
		1051	`--nn;`
		1052	`for (i=0; i < nn ; i++) {`
		1053	`pr[i] = qpr[i];`
		1054	`pi[i] = qpi[i];`
		1055	`}`
		1056	`}/while/`
		1057
		1058	`/* calculate the final zero and return */`
		1059	`cdivid(-pr[1], -pi[1], pr[0], pi[0], &zeror[d1], &zeroi[d1]);`
		1060
63171	ripley	1061	`vmaxset(vmax);`
61750	ripley	1062	`return;`
		1063	`}`
		1064
		1065
		1066	`/* Computes the derivative polynomial as the initial`
		1067	`* polynomial and computes l1 no-shift h polynomials. */`
		1068
		1069	`static void noshft(int l1)`
		1070	`{`
		1071	`int i, j, jj, n = nn - 1, nm1 = n - 1;`
		1072
		1073	`double t1, t2, xni;`
		1074
		1075	`for (i=0; i < n; i++) {`
		1076	`xni = (double)(nn - i - 1);`
		1077	`hr[i] = xni * pr[i] / n;`
		1078	`hi[i] = xni * pi[i] / n;`
		1079	`}`
		1080
		1081	`for (jj = 1; jj <= l1; jj++) {`
		1082
		1083	`if (hypot(hr[n-1], hi[n-1]) <=`
		1084	`eta * 10.0 * hypot(pr[n-1], pi[n-1])) {`
		1085	`/* If the constant term is essentially zero, */`
		1086	`/* shift h coefficients. */`
		1087
		1088	`for (i = 1; i <= nm1; i++) {`
		1089	`j = nn - i;`
		1090	`hr[j-1] = hr[j-2];`
		1091	`hi[j-1] = hi[j-2];`
		1092	`}`
		1093	`hr[0] = 0.;`
		1094	`hi[0] = 0.;`
		1095	`}`
		1096	`else {`
		1097	`cdivid(-pr[nn-1], -pi[nn-1], hr[n-1], hi[n-1], &tr, &ti);`
		1098	`for (i = 1; i <= nm1; i++) {`
		1099	`j = nn - i;`
		1100	`t1 = hr[j-2];`
		1101	`t2 = hi[j-2];`
		1102	`hr[j-1] = tr * t1 - ti * t2 + pr[j-1];`
		1103	`hi[j-1] = tr * t2 + ti * t1 + pi[j-1];`
		1104	`}`
		1105	`hr[0] = pr[0];`
		1106	`hi[0] = pi[0];`
		1107	`}`
		1108	`}`
		1109	`}`
		1110
		1111
		1112	`/* Computes l2 fixed-shift h polynomials and tests for convergence.`
		1113	`* initiates a variable-shift iteration and returns with the`
		1114	`* approximate zero if successful.`
		1115	`*/`
		1116	`static Rboolean fxshft(int l2, double zr, double zi)`
		1117	`{`
		1118	`/* l2 - limit of fixed shift steps`
		1119	`* zr,zi - approximate zero if convergence (result TRUE)`
		1120	`*`
		1121	`* Return value indicates convergence of stage 3 iteration`
		1122	`*`
		1123	`* Uses global (sr,si), nn, pr[], pi[], .. (all args of polyev() !)`
		1124	`*/`
		1125
		1126	`Rboolean pasd, bool, test;`
		1127	`static double svsi, svsr;`
		1128	`static int i, j, n;`
		1129	`static double oti, otr;`
		1130
		1131	`n = nn - 1;`
		1132
		1133	`/* evaluate p at s. */`
		1134
		1135	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
		1136
		1137	`test = TRUE;`
		1138	`pasd = FALSE;`
		1139
		1140	`/* calculate first t = -p(s)/h(s). */`
		1141
		1142	`calct(&bool);`
		1143
		1144	`/* main loop for one second stage step. */`
		1145
		1146	`for (j=1; j<=l2; j++) {`
		1147
		1148	`otr = tr;`
		1149	`oti = ti;`
		1150
		1151	`/* compute next h polynomial and new t. */`
		1152
		1153	`nexth(bool);`
		1154	`calct(&bool);`
		1155	`*zr = sr + tr;`
		1156	`*zi = si + ti;`
		1157
		1158	`/* test for convergence unless stage 3 has */`
		1159	`/* failed once or this is the last h polynomial. */`
		1160
		1161	`if (!bool && test && j != l2) {`
		1162	`if (hypot(tr - otr, ti - oti) >= hypot(zr, zi) * 0.5) {`
		1163	`pasd = FALSE;`
		1164	`}`
		1165	`else if (! pasd) {`
		1166	`pasd = TRUE;`
		1167	`}`
		1168	`else {`
		1169
		1170	`/* the weak convergence test has been */`
		1171	`/* passed twice, start the third stage */`
		1172	`/* iteration, after saving the current */`
		1173	`/* h polynomial and shift. */`
		1174
		1175	`for (i = 0; i < n; i++) {`
		1176	`shr[i] = hr[i];`
		1177	`shi[i] = hi[i];`
		1178	`}`
		1179	`svsr = sr;`
		1180	`svsi = si;`
		1181	`if (vrshft(10, zr, zi)) {`
		1182	`return TRUE;`
		1183	`}`
		1184
		1185	`/* the iteration failed to converge. */`
		1186	`/* turn off testing and restore */`
		1187	`/* h, s, pv and t. */`
		1188
		1189	`test = FALSE;`
		1190	`for (i=1 ; i<=n ; i++) {`
		1191	`hr[i-1] = shr[i-1];`
		1192	`hi[i-1] = shi[i-1];`
		1193	`}`
		1194	`sr = svsr;`
		1195	`si = svsi;`
		1196	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
		1197	`calct(&bool);`
		1198	`}`
		1199	`}`
		1200	`}`
		1201
		1202	`/* attempt an iteration with final h polynomial */`
		1203	`/* from second stage. */`
		1204
		1205	`return(vrshft(10, zr, zi));`
		1206	`}`
		1207
		1208
		1209	`/* carries out the third stage iteration.`
		1210	`*/`
		1211	`static Rboolean vrshft(int l3, double zr, double zi)`
		1212	`{`
		1213	`/* l3 - limit of steps in stage 3.`
		1214	`* zr,zi - on entry contains the initial iterate;`
		1215	`* if the iteration converges it contains`
		1216	`* the final iterate on exit.`
		1217	`* Returns TRUE if iteration converges`
		1218	`*`
		1219	`* Assign and uses GLOBAL sr, si`
		1220	`*/`
		1221	`Rboolean bool, b;`
		1222	`static int i, j;`
		1223	`static double r1, r2, mp, ms, tp, relstp;`
		1224	`static double omp;`
		1225
		1226	`b = FALSE;`
		1227	`sr = *zr;`
		1228	`si = *zi;`
		1229
		1230	`/* main loop for stage three */`
		1231
		1232	`for (i = 1; i <= l3; i++) {`
		1233
		1234	`/* evaluate p at s and test for convergence. */`
		1235	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
		1236
		1237	`mp = hypot(pvr, pvi);`
		1238	`ms = hypot(sr, si);`
		1239	`if (mp <= 20. * errev(nn, qpr, qpi, ms, mp, /are=/eta, mre)) {`
		1240	`goto L_conv;`
		1241	`}`
		1242
		1243	`/* polynomial value is smaller in value than */`
		1244	`/* a bound on the error in evaluating p, */`
		1245	`/* terminate the iteration. */`
		1246
		1247	`if (i != 1) {`
		1248
		1249	`if (!b && mp >= omp && relstp < .05) {`
		1250
		1251	`/* iteration has stalled. probably a */`
		1252	`/* cluster of zeros. do 5 fixed shift */`
		1253	`/* steps into the cluster to force */`
		1254	`/* one zero to dominate. */`
		1255
		1256	`tp = relstp;`
		1257	`b = TRUE;`
		1258	`if (relstp < eta)`
		1259	`tp = eta;`
		1260	`r1 = sqrt(tp);`
		1261	`r2 = sr * (r1 + 1.) - si * r1;`
		1262	`si = sr * r1 + si * (r1 + 1.);`
		1263	`sr = r2;`
		1264	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
		1265	`for (j = 1; j <= 5; ++j) {`
		1266	`calct(&bool);`
		1267	`nexth(bool);`
		1268	`}`
		1269	`omp = infin;`
		1270	`goto L10;`
		1271	`}`
		1272	`else {`
		1273
		1274	`/* exit if polynomial value */`
		1275	`/* increases significantly. */`
		1276
		1277	`if (mp * .1 > omp)`
		1278	`return FALSE;`
		1279	`}`
		1280	`}`
		1281	`omp = mp;`
		1282
		1283	`/* calculate next iterate. */`
		1284
		1285	`L10:`
		1286	`calct(&bool);`
		1287	`nexth(bool);`
		1288	`calct(&bool);`
		1289	`if (!bool) {`
		1290	`relstp = hypot(tr, ti) / hypot(sr, si);`
		1291	`sr += tr;`
		1292	`si += ti;`
		1293	`}`
		1294	`}`
		1295	`return FALSE;`
		1296
		1297	`L_conv:`
		1298	`*zr = sr;`
		1299	`*zi = si;`
		1300	`return TRUE;`
		1301	`}`
		1302
		1303	`static void calct(Rboolean *bool)`
		1304	`{`
		1305	`/* computes t = -p(s)/h(s).`
		1306	`* bool - logical, set true if h(s) is essentially zero. */`
		1307
		1308	`int n = nn - 1;`
		1309	`double hvi, hvr;`
		1310
		1311	`/* evaluate h(s). */`
		1312	`polyev(n, sr, si, hr, hi,`
		1313	`qhr, qhi, &hvr, &hvi);`
		1314
		1315	`bool = hypot(hvr, hvi) <= are 10. * hypot(hr[n-1], hi[n-1]);`
		1316	`if (!*bool) {`
		1317	`cdivid(-pvr, -pvi, hvr, hvi, &tr, &ti);`
		1318	`}`
		1319	`else {`
		1320	`tr = 0.;`
		1321	`ti = 0.;`
		1322	`}`
		1323	`}`
		1324
		1325	`static void nexth(Rboolean bool)`
		1326	`{`
		1327	`/* calculates the next shifted h polynomial.`
		1328	`* bool : if TRUE h(s) is essentially zero`
		1329	`*/`
		1330	`int j, n = nn - 1;`
		1331	`double t1, t2;`
		1332
		1333	`if (!bool) {`
		1334	`for (j=1; j < n; j++) {`
		1335	`t1 = qhr[j - 1];`
		1336	`t2 = qhi[j - 1];`
		1337	`hr[j] = tr * t1 - ti * t2 + qpr[j];`
		1338	`hi[j] = tr * t2 + ti * t1 + qpi[j];`
		1339	`}`
		1340	`hr[0] = qpr[0];`
		1341	`hi[0] = qpi[0];`
		1342	`}`
		1343	`else {`
		1344	`/* if h(s) is zero replace h with qh. */`
		1345
		1346	`for (j=1; j < n; j++) {`
		1347	`hr[j] = qhr[j-1];`
		1348	`hi[j] = qhi[j-1];`
		1349	`}`
		1350	`hr[0] = 0.;`
		1351	`hi[0] = 0.;`
		1352	`}`
		1353	`}`
		1354
		1355	`/--------------------- Independent Complex Polynomial Utilities ----------/`
		1356
		1357	`static`
		1358	`void polyev(int n,`
		1359	`double s_r, double s_i,`
		1360	`double p_r, double p_i,`
		1361	`double q_r, double q_i,`
		1362	`double v_r, double v_i)`
		1363	`{`
		1364	`/* evaluates a polynomial p at s by the horner recurrence`
		1365	`* placing the partial sums in q and the computed value in v_.`
		1366	`*/`
		1367	`int i;`
		1368	`double t;`
		1369
		1370	`q_r[0] = p_r[0];`
		1371	`q_i[0] = p_i[0];`
		1372	`*v_r = q_r[0];`
		1373	`*v_i = q_i[0];`
		1374	`for (i = 1; i < n; i++) {`
		1375	`t = v_r s_r - v_i s_i + p_r[i];`
		1376	`q_i[i] = v_i = v_r * s_i + v_i s_r + p_i[i];`
		1377	`q_r[i] = *v_r = t;`
		1378	`}`
		1379	`}`
		1380
		1381	`static`
		1382	`double errev(int n, double qr, double qi,`
		1383	`double ms, double mp, double a_re, double m_re)`
		1384	`{`
		1385	`/* bounds the error in evaluating the polynomial by the horner`
		1386	`* recurrence.`
		1387	`*`
		1388	`* qr,qi - the partial sum vectors`
		1389	`* ms - modulus of the point`
		1390	`* mp - modulus of polynomial value`
		1391	`* a_re,m_re - error bounds on complex addition and multiplication`
		1392	`*/`
		1393	`double e;`
		1394	`int i;`
		1395
		1396	`e = hypot(qr[0], qi[0]) * m_re / (a_re + m_re);`
		1397	`for (i=0; i < n; i++)`
		1398	`e = e*ms + hypot(qr[i], qi[i]);`
		1399
		1400	`return e * (a_re + m_re) - mp * m_re;`
		1401	`}`
		1402
		1403
		1404	`static`
		1405	`double cpoly_cauchy(int n, double pot, double q)`
		1406	`{`
		1407	`/* Computes a lower bound on the moduli of the zeros of a polynomial`
		1408	`* pot[1:nn] is the modulus of the coefficients.`
		1409	`*/`
		1410	`double f, x, delf, dx, xm;`
		1411	`int i, n1 = n - 1;`
		1412
		1413	`pot[n1] = -pot[n1];`
		1414
		1415	`/* compute upper estimate of bound. */`
		1416
		1417	`x = exp((log(-pot[n1]) - log(pot[0])) / (double) n1);`
		1418
		1419	`/* if newton step at the origin is better, use it. */`
		1420
		1421	`if (pot[n1-1] != 0.) {`
		1422	`xm = -pot[n1] / pot[n1-1];`
		1423	`if (xm < x)`
		1424	`x = xm;`
		1425	`}`
		1426
		1427	`/* chop the interval (0,x) unitl f le 0. */`
		1428
		1429	`for(;;) {`
		1430	`xm = x * 0.1;`
		1431	`f = pot[0];`
		1432	`for (i = 1; i < n; i++)`
		1433	`f = f * xm + pot[i];`
		1434	`if (f <= 0.0) {`
		1435	`break;`
		1436	`}`
		1437	`x = xm;`
		1438	`}`
		1439
		1440	`dx = x;`
		1441
		1442	`/* do Newton iteration until x converges to two decimal places. */`
		1443
		1444	`while (fabs(dx / x) > 0.005) {`
		1445	`q[0] = pot[0];`
		1446	`for(i = 1; i < n; i++)`
		1447	`q[i] = q[i-1] * x + pot[i];`
		1448	`f = q[n1];`
		1449	`delf = q[0];`
		1450	`for(i = 1; i < n1; i++)`
		1451	`delf = delf * x + q[i];`
		1452	`dx = f / delf;`
		1453	`x -= dx;`
		1454	`}`
		1455	`return x;`
		1456	`}`
		1457
		1458	`static`
		1459	`double cpoly_scale(int n, double *pot,`
		1460	`double eps, double BIG, double small, double base)`
		1461	`{`
		1462	`/* Returns a scale factor to multiply the coefficients of the polynomial.`
		1463	`* The scaling is done to avoid overflow and to avoid`
		1464	`* undetected underflow interfering with the convergence criterion.`
		1465	`* The factor is a power of the base.`
		1466
		1467	`* pot[1:n] : modulus of coefficients of p`
		1468	`* eps,BIG,`
		1469	`* small,base - constants describing the floating point arithmetic.`
		1470	`*/`
		1471
		1472	`int i, ell;`
		1473	`double x, high, sc, lo, min_, max_;`
		1474
		1475	`/* find largest and smallest moduli of coefficients. */`
		1476	`high = sqrt(BIG);`
		1477	`lo = small / eps;`
		1478	`max_ = 0.;`
		1479	`min_ = BIG;`
		1480	`for (i = 0; i < n; i++) {`
		1481	`x = pot[i];`
		1482	`if (x > max_) max_ = x;`
		1483	`if (x != 0. && x < min_)`
		1484	`min_ = x;`
		1485	`}`
		1486
		1487	`/* scale only if there are very large or very small components. */`
		1488
		1489	`if (min_ < lo \|\| max_ > high) {`
		1490	`x = lo / min_;`
		1491	`if (x <= 1.)`
		1492	`sc = 1. / (sqrt(max_) * sqrt(min_));`
		1493	`else {`
		1494	`sc = x;`
		1495	`if (BIG / sc > max_)`
		1496	`sc = 1.0;`
		1497	`}`
		1498	`ell = (int) (log(sc) / log(base) + 0.5);`
		1499	`return R_pow_di(base, ell);`
		1500	`}`
		1501	`else return 1.0;`
		1502	`}`
		1503
		1504
		1505	`static`
		1506	`void cdivid(double ar, double ai, double br, double bi,`
		1507	`double cr, double ci)`
		1508	`{`
		1509	`/* complex division c = a/b, i.e., (cr +ici) = (ar +iai) / (br +i*bi),`
		1510	`avoiding overflow. */`
		1511
		1512	`double d, r;`
		1513
		1514	`if (br == 0. && bi == 0.) {`
		1515	`/* division by zero, c = infinity. */`
		1516	`cr = ci = R_PosInf;`
		1517	`}`
		1518	`else if (fabs(br) >= fabs(bi)) {`
		1519	`r = bi / br;`
		1520	`d = br + r * bi;`
		1521	`cr = (ar + ai r) / d;`
		1522	`ci = (ai - ar r) / d;`
		1523	`}`
		1524	`else {`
		1525	`r = br / bi;`
		1526	`d = bi + r * br;`
		1527	`cr = (ar r + ai) / d;`
		1528	`ci = (ai r - ar) / d;`
		1529	`}`
		1530	`}`
		1531
		1532	`/* static double cpoly_cmod(double r, double i)`
		1533	`* --> replaced by hypot() everywhere`
		1534	`*/`

The R Project SVN R

(root)/trunk/src/main/complex.c – Rev 77458