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`
87891	ripley	3	`* Copyright (C) 2000-2025 The R Core Team`
77395	maechler	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_* */`
88866	maechler	66	`#include <complex.h> /* incl "complex math.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
83446	ripley	75	`attribute_hidden SEXP 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
83446	ripley	172	`attribute_hidden SEXP 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
83446	ripley	245	`attribute_hidden SEXP 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`
83446	ripley	360	`attribute_hidden void 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 */`
88866	maechler	465	`double x = creal(z), y = cimag(z),`
		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
88866	maechler	588	`#ifdef HAVE_CABS`
		589	`# define R_CABS(Z) cabs(Z)`
		590	`#else`
		591	`# define R_CABS(Z) hypot(creal(Z), cimag(Z))`
		592	`#endif`
		593
		594
87891	ripley	595	`static bool cmath1(double complex (*f)(double complex),`
73732	luke	596	`const Rcomplex x, Rcomplex y, R_xlen_t n)`
54273	ripley	597	`{`
59090	ripley	598	`R_xlen_t i;`
87891	ripley	599	`bool naflag = false;`
54273	ripley	600	`for (i = 0 ; i < n ; i++) {`
		601	`if (ISNA(x[i].r) \|\| ISNA(x[i].i)) {`
		602	`y[i].r = NA_REAL; y[i].i = NA_REAL;`
		603	`} else {`
		604	`SET_C99_COMPLEX(y, i, f(toC99(x + i)));`
55516	ripley	605	`if ( (ISNAN(y[i].r) \|\| ISNAN(y[i].i)) &&`
87891	ripley	606	`!(ISNAN(x[i].r) \|\| ISNAN(x[i].i)) ) naflag = true;`
54273	ripley	607	`}`
		608	`}`
		609	`return naflag;`
		610	`}`
		611
83446	ripley	612	`attribute_hidden SEXP complex_math1(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	613	`{`
580	ihaka	614	`SEXP x, y;`
59048	ripley	615	`R_xlen_t n;`
87891	ripley	616	`bool naflag = false;`
54273	ripley	617
6098	pd	618	`PROTECT(x = CAR(args));`
68586	ripley	619	`n = XLENGTH(x);`
6098	pd	620	`PROTECT(y = allocVector(CPLXSXP, n));`
2	r	621
73732	luke	622	`const Rcomplex *px = COMPLEX_RO(x);`
		623	`Rcomplex *py = COMPLEX(y);`
		624
580	ihaka	625	`switch (PRIMVAL(op)) {`
73732	luke	626	`case 10003: naflag = cmath1(clog, px, py, n); break;`
77395	maechler	627	`// 1: floor`
		628	`// 2: ceil[ing]`
73732	luke	629	`case 3: naflag = cmath1(csqrt, px, py, n); break;`
77395	maechler	630	`// 4: sign`
73732	luke	631	`case 10: naflag = cmath1(cexp, px, py, n); break;`
77395	maechler	632	`// 11: expm1`
		633	`// 12: log1p`
73732	luke	634	`case 20: naflag = cmath1(ccos, px, py, n); break;`
		635	`case 21: naflag = cmath1(csin, px, py, n); break;`
		636	`case 22: naflag = cmath1(z_tan, px, py, n); break;`
		637	`case 23: naflag = cmath1(z_acos, px, py, n); break;`
		638	`case 24: naflag = cmath1(z_asin, px, py, n); break;`
		639	`case 25: naflag = cmath1(z_atan, px, py, n); break;`
77395	maechler	640
73732	luke	641	`case 30: naflag = cmath1(ccosh, px, py, n); break;`
		642	`case 31: naflag = cmath1(csinh, px, py, n); break;`
		643	`case 32: naflag = cmath1(ctanh, px, py, n); break;`
		644	`case 33: naflag = cmath1(z_acosh, px, py, n); break;`
		645	`case 34: naflag = cmath1(z_asinh, px, py, n); break;`
		646	`case 35: naflag = cmath1(z_atanh, px, py, n); break;`
54174	ripley	647
580	ihaka	648	`default:`
41725	ripley	649	`/* such as sign, gamma */`
		650	`errorcall(call, _("unimplemented complex function"));`
580	ihaka	651	`}`
6098	pd	652	`if (naflag)`
55518	ripley	653	`warningcall(call, "NaNs produced in function \"%s\"", PRIMNAME(op));`
69108	luke	654	`SHALLOW_DUPLICATE_ATTRIB(y, x);`
6098	pd	655	`UNPROTECT(2);`
580	ihaka	656	`return y;`
2	r	657	`}`
		658
54262	ripley	659	`static void z_rround(Rcomplex r, Rcomplex x, Rcomplex *p)`
		660	`{`
54776	ripley	661	`r->r = fround(x->r, p->r);`
		662	`r->i = fround(x->i, p->r);`
54262	ripley	663	`}`
		664
54174	ripley	665	`static void z_prec(Rcomplex r, Rcomplex x, Rcomplex *p)`
		666	`{`
		667	`z_prec_r(r, x, p->r);`
		668	`}`
		669
		670	`static void z_logbase(Rcomplex r, Rcomplex z, Rcomplex *base)`
		671	`{`
		672	`double complex dz = toC99(z), dbase = toC99(base);`
		673	`SET_C99_COMPLEX(r, 0, clog(dz)/clog(dbase));`
		674	`}`
		675
		676	`static void z_atan2(Rcomplex r, Rcomplex csn, Rcomplex *ccs)`
		677	`{`
		678	`double complex dr, dcsn = toC99(csn), dccs = toC99(ccs);`
		679	`if (dccs == 0) {`
		680	`if(dcsn == 0) {`
54256	ripley	681	`r->r = NA_REAL; r->i = NA_REAL; /* Why not R_NaN? */`
54174	ripley	682	`return;`
54293	ripley	683	`} else {`
		684	`double y = creal(dcsn);`
		685	`if (ISNAN(y)) dr = y;`
		686	`else dr = ((y >= 0) ? M_PI_2 : -M_PI_2);`
		687	`}`
54174	ripley	688	`} else {`
		689	`dr = catan(dcsn / dccs);`
		690	`if(creal(dccs) < 0) dr += M_PI;`
		691	`if(creal(dr) > M_PI) dr -= 2 * M_PI;`
		692	`}`
		693	`SET_C99_COMPLEX(r, 0, dr);`
		694	`}`
		695
54273	ripley	696
		697	`/* Complex Functions of Two Arguments */`
		698
		699	`typedef void (cm2_fun)(Rcomplex , Rcomplex , Rcomplex );`
83446	ripley	700	`attribute_hidden SEXP complex_math2(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	701	`{`
54273	ripley	702	`cm2_fun f;`
		703	`switch (PRIMVAL(op)) {`
		704	`case 0: /* atan2 */`
		705	`f = z_atan2; break;`
		706	`case 10001: /* round */`
		707	`f = z_rround; break;`
88828	maechler	708	`case 10002: /* passed from do_log1arg */`
88866	maechler	709	`case 10010: /* " " " */`
54273	ripley	710	`case 10003: /* passed from do_log */`
		711	`f = z_logbase; break;`
		712	`case 10004: /* signif */`
		713	`f = z_prec; break;`
		714	`default:`
70830	luke	715	`error_return(_("unimplemented complex function"));`
54273	ripley	716	`}`
		717
88866	maechler	718	`SEXP`
		719	`sa = PROTECT(coerceVector(CAR(args), CPLXSXP)),`
		720	`sb = PROTECT(coerceVector(CADR(args), CPLXSXP));`
		721	`R_xlen_t na = XLENGTH(sa), nb = XLENGTH(sb), n;`
67024	luke	722	`if ((na == 0) \|\| (nb == 0)) {`
68923	ripley	723	`UNPROTECT(2);`
		724	`return(allocVector(CPLXSXP, 0));`
67024	luke	725	`}`
580	ihaka	726	`n = (na < nb) ? nb : na;`
88866	maechler	727	`SEXP sy = PROTECT(allocVector(CPLXSXP, n));`
		728	`const Rcomplex`
		729	`*a = COMPLEX_RO(sa),`
		730	`*b = COMPLEX_RO(sb);`
		731	`Rcomplex ai, bi, *y = COMPLEX(sy);`
		732	`R_xlen_t i, ia, ib;`
		733	`bool naflag = false;`
68663	luke	734	`MOD_ITERATE2(n, na, nb, i, ia, ib, {`
		735	`ai = a[ia]; bi = b[ib];`
6098	pd	736	`if(ISNA(ai.r) && ISNA(ai.i) &&`
		737	`ISNA(bi.r) && ISNA(bi.i)) {`
54273	ripley	738	`y[i].r = NA_REAL; y[i].i = NA_REAL;`
54170	ripley	739	`} else {`
6098	pd	740	`f(&y[i], &ai, &bi);`
55516	ripley	741	`if ( (ISNAN(y[i].r) \|\| ISNAN(y[i].i)) &&`
		742	`!(ISNAN(ai.r) \|\| ISNAN(ai.i) \|\| ISNAN(bi.r) \|\| ISNAN(bi.i)) )`
87891	ripley	743	`naflag = true;`
2	r	744	`}`
68663	luke	745	`});`
6098	pd	746	`if (naflag)`
70830	luke	747	`warning("NaNs produced in function \"%s\"", PRIMNAME(op));`
580	ihaka	748	`if(n == na) {`
69108	luke	749	`SHALLOW_DUPLICATE_ATTRIB(sy, sa);`
54170	ripley	750	`} else if(n == nb) {`
69108	luke	751	`SHALLOW_DUPLICATE_ATTRIB(sy, sb);`
580	ihaka	752	`}`
		753	`UNPROTECT(3);`
		754	`return sy;`
2	r	755	`}`
		756
83446	ripley	757	`attribute_hidden SEXP do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	758	`{`
1839	ihaka	759	`/* complex(length, real, imaginary) */`
580	ihaka	760	`SEXP ans, re, im;`
59048	ripley	761	`R_xlen_t i, na, nr, ni;`
69326	luke	762
		763	`checkArity(op, args);`
580	ihaka	764	`na = asInteger(CAR(args));`
		765	`if(na == NA_INTEGER \|\| na < 0)`
41686	ripley	766	`error(_("invalid length"));`
580	ihaka	767	`PROTECT(re = coerceVector(CADR(args), REALSXP));`
		768	`PROTECT(im = coerceVector(CADDR(args), REALSXP));`
59048	ripley	769	`nr = XLENGTH(re);`
		770	`ni = XLENGTH(im);`
1839	ihaka	771	`/* is always true: if (na >= 0) {*/`
		772	`na = (nr > na) ? nr : na;`
		773	`na = (ni > na) ? ni : na;`
		774	`/* }*/`
580	ihaka	775	`ans = allocVector(CPLXSXP, na);`
73732	luke	776	`Rcomplex *pans = COMPLEX(ans);`
580	ihaka	777	`for(i=0 ; i<na ; i++) {`
73732	luke	778	`pans[i].r = 0;`
		779	`pans[i].i = 0;`
580	ihaka	780	`}`
		781	`UNPROTECT(2);`
		782	`if(na > 0 && nr > 0) {`
73732	luke	783	`const double *p_re = REAL_RO(re);`
580	ihaka	784	`for(i=0 ; i<na ; i++)`
73732	luke	785	`pans[i].r = p_re[i%nr];`
580	ihaka	786	`}`
		787	`if(na > 0 && ni > 0) {`
73732	luke	788	`const double *p_im = REAL_RO(im);`
580	ihaka	789	`for(i=0 ; i<na ; i++)`
73732	luke	790	`pans[i].i = p_im[i%ni];`
580	ihaka	791	`}`
		792	`return ans;`
2	r	793	`}`
		794
61750	ripley	795	`static void R_cpolyroot(double opr, double opi, int *degree,`
87891	ripley	796	`double zeror, double zeroi, bool *fail);`
580	ihaka	797
83446	ripley	798	`attribute_hidden SEXP do_polyroot(SEXP call, SEXP op, SEXP args, SEXP rho)`
2	r	799	`{`
580	ihaka	800	`checkArity(op, args);`
88866	maechler	801	`SEXP z = CAR(args);`
580	ihaka	802	`switch(TYPEOF(z)) {`
		803	`case CPLXSXP:`
		804	`PROTECT(z);`
		805	`break;`
		806	`case REALSXP:`
		807	`case INTSXP:`
		808	`case LGLSXP:`
		809	`PROTECT(z = coerceVector(z, CPLXSXP));`
		810	`break;`
		811	`default:`
31908	ripley	812	`UNIMPLEMENTED_TYPE("polyroot", z);`
580	ihaka	813	`}`
88866	maechler	814
		815	`int i, n;`
61750	ripley	816	`#ifdef LONG_VECTOR_SUPPORT`
		817	`R_xlen_t nn = XLENGTH(z);`
		818	`if (nn > R_SHORT_LEN_MAX) error("long vectors are not supported");`
		819	`n = (int) nn;`
		820	`#else`
		821	`n = LENGTH(z);`
		822	`#endif`
73732	luke	823	`const Rcomplex *pz = COMPLEX_RO(z);`
88866	maechler	824	`int degree = 0; // := max{i; z[i] != 0}`
18490	ripley	825	`for(i = 0; i < n; i++) {`
73732	luke	826	`if(pz[i].r!= 0.0 \|\| pz[i].i != 0.0) degree = i;`
18490	ripley	827	`}`
		828	`n = degree + 1; /* omit trailing zeroes */`
88866	maechler	829	`SEXP r;`
580	ihaka	830	`if(degree >= 1) {`
88866	maechler	831	`SEXP zr, zi, rr, ri;`
580	ihaka	832	`PROTECT(rr = allocVector(REALSXP, n));`
		833	`PROTECT(ri = allocVector(REALSXP, n));`
		834	`PROTECT(zr = allocVector(REALSXP, n));`
		835	`PROTECT(zi = allocVector(REALSXP, n));`
2	r	836
73732	luke	837	`double *p_rr = REAL(rr);`
		838	`double *p_ri = REAL(ri);`
		839	`double *p_zr = REAL(zr);`
		840	`double *p_zi = REAL(zi);`
		841
61750	ripley	842	`for(i = 0 ; i < n ; i++) {`
73732	luke	843	`if(!R_FINITE(pz[i].r) \|\| !R_FINITE(pz[i].i))`
41686	ripley	844	`error(_("invalid polynomial coefficient"));`
73732	luke	845	`p_zr[degree-i] = pz[i].r;`
		846	`p_zi[degree-i] = pz[i].i;`
2	r	847	`}`
88866	maechler	848	`bool fail;`
73732	luke	849	`R_cpolyroot(p_zr, p_zi, &degree, p_rr, p_ri, &fail);`
41686	ripley	850	`if(fail) error(_("root finding code failed"));`
580	ihaka	851	`UNPROTECT(2);`
		852	`r = allocVector(CPLXSXP, degree);`
73732	luke	853	`Rcomplex *pr = COMPLEX(r);`
61750	ripley	854	`for(i = 0 ; i < degree ; i++) {`
73732	luke	855	`pr[i].r = p_rr[i];`
		856	`pr[i].i = p_ri[i];`
2	r	857	`}`
580	ihaka	858	`UNPROTECT(3);`
		859	`}`
		860	`else {`
		861	`UNPROTECT(1);`
		862	`r = allocVector(CPLXSXP, 0);`
		863	`}`
		864	`return r;`
2	r	865	`}`
61750	ripley	866
		867	`/* Formerly src/appl/cpoly.c:`
		868	`*`
		869	`* Copyright (C) 1997-1998 Ross Ihaka`
		870	`* Copyright (C) 1999-2001 R Core Team`
		871	`*`
		872	`* cpoly finds the zeros of a complex polynomial.`
		873	`*`
		874	`* On Entry`
		875	`*`
		876	`* opr, opi - double precision vectors of real and`
		877	`* imaginary parts of the coefficients in`
		878	`* order of decreasing powers.`
		879	`*`
		880	`* degree - int degree of polynomial.`
		881	`*`
		882	`*`
		883	`* On Return`
		884	`*`
		885	`* zeror, zeroi - output double precision vectors of`
		886	`* real and imaginary parts of the zeros.`
		887	`*`
		888	`* fail - output int parameter, true only if`
		889	`* leading coefficient is zero or if cpoly`
		890	`* has found fewer than degree zeros.`
		891	`*`
		892	`* The program has been written to reduce the chance of overflow`
		893	`* occurring. If it does occur, there is still a possibility that`
		894	`* the zerofinder will work provided the overflowed quantity is`
		895	`* replaced by a large number.`
		896	`*`
		897	`* This is a C translation of the following.`
		898	`*`
		899	`* TOMS Algorithm 419`
		900	`* Jenkins and Traub.`
		901	`* Comm. ACM 15 (1972) 97-99.`
		902	`*`
		903	`* Ross Ihaka`
		904	`* February 1997`
		905	`*/`
		906
		907	`#include <R_ext/Arith.h> /* for declaration of hypot */`
		908	`#include <R_ext/Memory.h> /* for declaration of R_alloc */`
		909
		910	`#include <float.h> /* for FLT_RADIX */`
		911
		912	`#include <Rmath.h> /* for R_pow_di */`
		913
87891	ripley	914	`static void calct(bool *);`
		915	`static bool fxshft(int, double , double );`
		916	`static bool vrshft(int, double , double );`
		917	`static void nexth(bool);`
61750	ripley	918	`static void noshft(int);`
		919
		920	`static void polyev(int, double, double,`
		921	`double , double , double , double , double , double );`
		922	`static double errev(int, double , double , double, double, double, double);`
		923	`static double cpoly_cauchy(int, double , double );`
		924	`static double cpoly_scale(int, double *, double, double, double, double);`
		925	`static void cdivid(double, double, double, double, double , double );`
		926
		927	`/* Global Variables (too many!) */`
		928
		929	`static int nn;`
		930	`static double pr, pi, hr, hi, qpr, qpi, qhr, qhi, shr, shi;`
		931	`static double sr, si;`
		932	`static double tr, ti;`
		933	`static double pvr, pvi;`
		934
		935	`static const double eta = DBL_EPSILON;`
		936	`static const double are = /* eta = */DBL_EPSILON;`
		937	`static const double mre = 2. * M_SQRT2 * /* eta, i.e. */DBL_EPSILON;`
		938	`static const double infin = DBL_MAX;`
		939
		940	`static void R_cpolyroot(double opr, double opi, int *degree,`
87891	ripley	941	`double zeror, double zeroi, bool *fail)`
61750	ripley	942	`{`
		943	`static const double smalno = DBL_MIN;`
		944	`static const double base = (double)FLT_RADIX;`
		945	`static int d_n, i, i1, i2;`
		946	`static double zi, zr, xx, yy;`
		947	`static double bnd, xxx;`
87891	ripley	948	`bool conv;`
61750	ripley	949	`int d1;`
		950	`double *tmp;`
		951	`static const double cosr =/* cos 94 */ -0.06975647374412529990;`
		952	`static const double sinr =/* sin 94 */ 0.99756405025982424767;`
		953	`xx = M_SQRT1_2;/* 1/sqrt(2) = 0.707.... */`
		954
		955	`yy = -xx;`
87891	ripley	956	`*fail = false;`
61750	ripley	957
		958	`nn = *degree;`
		959	`d1 = nn - 1;`
		960
		961	`/* algorithm fails if the leading coefficient is zero. */`
		962
		963	`if (opr[0] == 0. && opi[0] == 0.) {`
87891	ripley	964	`*fail = true;`
61750	ripley	965	`return;`
		966	`}`
		967
		968	`/* remove the zeros at the origin if any. */`
		969
		970	`while (opr[nn] == 0. && opi[nn] == 0.) {`
		971	`d_n = d1-nn+1;`
		972	`zeror[d_n] = 0.;`
		973	`zeroi[d_n] = 0.;`
		974	`nn--;`
		975	`}`
		976	`nn++;`
		977	`/-- Now, global var. nn := #{coefficients} = (relevant degree)+1 /`
		978
		979	`if (nn == 1) return;`
		980
		981	`/* Use a single allocation as these as small */`
63181	ripley	982	`const void *vmax = vmaxget();`
61750	ripley	983	`tmp = (double ) R_alloc((size_t) (10nn), sizeof(double));`
		984	`pr = tmp; pi = tmp + nn; hr = tmp + 2nn; hi = tmp + 3nn;`
		985	`qpr = tmp + 4nn; qpi = tmp + 5nn; qhr = tmp + 6nn; qhi = tmp + 7nn;`
		986	`shr = tmp + 8nn; shi = tmp + 9nn;`
68923	ripley	987
61750	ripley	988	`/* make a copy of the coefficients and shr[] = \| p[] \| */`
		989	`for (i = 0; i < nn; i++) {`
		990	`pr[i] = opr[i];`
		991	`pi[i] = opi[i];`
		992	`shr[i] = hypot(pr[i], pi[i]);`
		993	`}`
		994
		995	`/* scale the polynomial with factor 'bnd'. */`
		996	`bnd = cpoly_scale(nn, shr, eta, infin, smalno, base);`
		997	`if (bnd != 1.) {`
		998	`for (i=0; i < nn; i++) {`
		999	`pr[i] *= bnd;`
		1000	`pi[i] *= bnd;`
		1001	`}`
		1002	`}`
		1003
		1004	`/* start the algorithm for one zero */`
		1005
		1006	`while (nn > 2) {`
		1007
		1008	`/* calculate bnd, a lower bound on the modulus of the zeros. */`
		1009
		1010	`for (i=0 ; i < nn ; i++)`
		1011	`shr[i] = hypot(pr[i], pi[i]);`
		1012	`bnd = cpoly_cauchy(nn, shr, shi);`
		1013
		1014	`/* outer loop to control 2 major passes */`
		1015	`/* with different sequences of shifts */`
		1016
		1017	`for (i1 = 1; i1 <= 2; i1++) {`
		1018
		1019	`/* first stage calculation, no shift */`
		1020
		1021	`noshft(5);`
		1022
		1023	`/* inner loop to select a shift */`
		1024	`for (i2 = 1; i2 <= 9; i2++) {`
		1025
		1026	`/* shift is chosen with modulus bnd */`
		1027	`/* and amplitude rotated by 94 degrees */`
		1028	`/* from the previous shift */`
		1029
		1030	`xxx= cosr * xx - sinr * yy;`
		1031	`yy = sinr * xx + cosr * yy;`
		1032	`xx = xxx;`
		1033	`sr = bnd * xx;`
		1034	`si = bnd * yy;`
		1035
		1036	`/* second stage calculation, fixed shift */`
		1037
		1038	`conv = fxshft(i2 * 10, &zr, &zi);`
		1039	`if (conv)`
		1040	`goto L10;`
		1041	`}`
		1042	`}`
		1043
		1044	`/* the zerofinder has failed on two major passes */`
		1045	`/* return empty handed */`
		1046
87891	ripley	1047	`*fail = true;`
63171	ripley	1048	`vmaxset(vmax);`
61750	ripley	1049	`return;`
		1050
		1051	`/* the second stage jumps directly to the third stage iteration.`
		1052	`* if successful, the zero is stored and the polynomial deflated.`
		1053	`*/`
		1054	`L10:`
		1055	`d_n = d1+2 - nn;`
		1056	`zeror[d_n] = zr;`
		1057	`zeroi[d_n] = zi;`
		1058	`--nn;`
		1059	`for (i=0; i < nn ; i++) {`
		1060	`pr[i] = qpr[i];`
		1061	`pi[i] = qpi[i];`
		1062	`}`
		1063	`}/while/`
		1064
		1065	`/* calculate the final zero and return */`
		1066	`cdivid(-pr[1], -pi[1], pr[0], pi[0], &zeror[d1], &zeroi[d1]);`
		1067
63171	ripley	1068	`vmaxset(vmax);`
61750	ripley	1069	`return;`
		1070	`}`
		1071
		1072
		1073	`/* Computes the derivative polynomial as the initial`
		1074	`* polynomial and computes l1 no-shift h polynomials. */`
		1075
		1076	`static void noshft(int l1)`
		1077	`{`
		1078	`int i, j, jj, n = nn - 1, nm1 = n - 1;`
		1079
		1080	`double t1, t2, xni;`
		1081
		1082	`for (i=0; i < n; i++) {`
		1083	`xni = (double)(nn - i - 1);`
		1084	`hr[i] = xni * pr[i] / n;`
		1085	`hi[i] = xni * pi[i] / n;`
		1086	`}`
		1087
		1088	`for (jj = 1; jj <= l1; jj++) {`
		1089
		1090	`if (hypot(hr[n-1], hi[n-1]) <=`
		1091	`eta * 10.0 * hypot(pr[n-1], pi[n-1])) {`
		1092	`/* If the constant term is essentially zero, */`
		1093	`/* shift h coefficients. */`
		1094
		1095	`for (i = 1; i <= nm1; i++) {`
		1096	`j = nn - i;`
		1097	`hr[j-1] = hr[j-2];`
		1098	`hi[j-1] = hi[j-2];`
		1099	`}`
		1100	`hr[0] = 0.;`
		1101	`hi[0] = 0.;`
		1102	`}`
		1103	`else {`
		1104	`cdivid(-pr[nn-1], -pi[nn-1], hr[n-1], hi[n-1], &tr, &ti);`
		1105	`for (i = 1; i <= nm1; i++) {`
		1106	`j = nn - i;`
		1107	`t1 = hr[j-2];`
		1108	`t2 = hi[j-2];`
		1109	`hr[j-1] = tr * t1 - ti * t2 + pr[j-1];`
		1110	`hi[j-1] = tr * t2 + ti * t1 + pi[j-1];`
		1111	`}`
		1112	`hr[0] = pr[0];`
		1113	`hi[0] = pi[0];`
		1114	`}`
		1115	`}`
		1116	`}`
		1117
		1118
		1119	`/* Computes l2 fixed-shift h polynomials and tests for convergence.`
		1120	`* initiates a variable-shift iteration and returns with the`
		1121	`* approximate zero if successful.`
		1122	`*/`
87891	ripley	1123	`static bool fxshft(int l2, double zr, double zi)`
61750	ripley	1124	`{`
		1125	`/* l2 - limit of fixed shift steps`
		1126	`* zr,zi - approximate zero if convergence (result TRUE)`
		1127	`*`
		1128	`* Return value indicates convergence of stage 3 iteration`
		1129	`*`
		1130	`* Uses global (sr,si), nn, pr[], pi[], .. (all args of polyev() !)`
		1131	`*/`
		1132
87891	ripley	1133	`bool pasd, h_s_0, test;`
61750	ripley	1134	`static double svsi, svsr;`
		1135	`static int i, j, n;`
		1136	`static double oti, otr;`
		1137
		1138	`n = nn - 1;`
		1139
		1140	`/* evaluate p at s. */`
		1141
		1142	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
		1143
87891	ripley	1144	`test = true;`
		1145	`pasd = false;`
61750	ripley	1146
		1147	`/* calculate first t = -p(s)/h(s). */`
		1148
82945	maechler	1149	`calct(&h_s_0);`
61750	ripley	1150
		1151	`/* main loop for one second stage step. */`
		1152
		1153	`for (j=1; j<=l2; j++) {`
		1154
		1155	`otr = tr;`
		1156	`oti = ti;`
		1157
		1158	`/* compute next h polynomial and new t. */`
		1159
82945	maechler	1160	`nexth(h_s_0);`
		1161	`calct(&h_s_0);`
61750	ripley	1162	`*zr = sr + tr;`
		1163	`*zi = si + ti;`
		1164
		1165	`/* test for convergence unless stage 3 has */`
		1166	`/* failed once or this is the last h polynomial. */`
		1167
82945	maechler	1168	`if (!h_s_0 && test && j != l2) {`
61750	ripley	1169	`if (hypot(tr - otr, ti - oti) >= hypot(zr, zi) * 0.5) {`
87891	ripley	1170	`pasd = false;`
61750	ripley	1171	`}`
		1172	`else if (! pasd) {`
87891	ripley	1173	`pasd = true;`
61750	ripley	1174	`}`
		1175	`else {`
		1176
		1177	`/* the weak convergence test has been */`
		1178	`/* passed twice, start the third stage */`
		1179	`/* iteration, after saving the current */`
		1180	`/* h polynomial and shift. */`
		1181
		1182	`for (i = 0; i < n; i++) {`
		1183	`shr[i] = hr[i];`
		1184	`shi[i] = hi[i];`
		1185	`}`
		1186	`svsr = sr;`
		1187	`svsi = si;`
		1188	`if (vrshft(10, zr, zi)) {`
87891	ripley	1189	`return true;`
61750	ripley	1190	`}`
		1191
		1192	`/* the iteration failed to converge. */`
		1193	`/* turn off testing and restore */`
		1194	`/* h, s, pv and t. */`
		1195
87891	ripley	1196	`test = false;`
61750	ripley	1197	`for (i=1 ; i<=n ; i++) {`
		1198	`hr[i-1] = shr[i-1];`
		1199	`hi[i-1] = shi[i-1];`
		1200	`}`
		1201	`sr = svsr;`
		1202	`si = svsi;`
		1203	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
82945	maechler	1204	`calct(&h_s_0);`
61750	ripley	1205	`}`
		1206	`}`
		1207	`}`
		1208
		1209	`/* attempt an iteration with final h polynomial */`
		1210	`/* from second stage. */`
		1211
		1212	`return(vrshft(10, zr, zi));`
		1213	`}`
		1214
		1215
		1216	`/* carries out the third stage iteration.`
		1217	`*/`
87891	ripley	1218	`static bool vrshft(int l3, double zr, double zi)`
61750	ripley	1219	`{`
		1220	`/* l3 - limit of steps in stage 3.`
		1221	`* zr,zi - on entry contains the initial iterate;`
		1222	`* if the iteration converges it contains`
		1223	`* the final iterate on exit.`
		1224	`* Returns TRUE if iteration converges`
		1225	`*`
		1226	`* Assign and uses GLOBAL sr, si`
		1227	`*/`
87891	ripley	1228	`bool h_s_0, b;`
61750	ripley	1229	`static int i, j;`
		1230	`static double r1, r2, mp, ms, tp, relstp;`
		1231	`static double omp;`
		1232
87891	ripley	1233	`b = false;`
61750	ripley	1234	`sr = *zr;`
		1235	`si = *zi;`
		1236
		1237	`/* main loop for stage three */`
		1238
		1239	`for (i = 1; i <= l3; i++) {`
		1240
		1241	`/* evaluate p at s and test for convergence. */`
		1242	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
		1243
		1244	`mp = hypot(pvr, pvi);`
		1245	`ms = hypot(sr, si);`
		1246	`if (mp <= 20. * errev(nn, qpr, qpi, ms, mp, /are=/eta, mre)) {`
		1247	`goto L_conv;`
		1248	`}`
		1249
		1250	`/* polynomial value is smaller in value than */`
		1251	`/* a bound on the error in evaluating p, */`
		1252	`/* terminate the iteration. */`
		1253
		1254	`if (i != 1) {`
		1255
		1256	`if (!b && mp >= omp && relstp < .05) {`
		1257
		1258	`/* iteration has stalled. probably a */`
		1259	`/* cluster of zeros. do 5 fixed shift */`
		1260	`/* steps into the cluster to force */`
		1261	`/* one zero to dominate. */`
		1262
		1263	`tp = relstp;`
87891	ripley	1264	`b = true;`
61750	ripley	1265	`if (relstp < eta)`
		1266	`tp = eta;`
		1267	`r1 = sqrt(tp);`
		1268	`r2 = sr * (r1 + 1.) - si * r1;`
		1269	`si = sr * r1 + si * (r1 + 1.);`
		1270	`sr = r2;`
		1271	`polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);`
		1272	`for (j = 1; j <= 5; ++j) {`
82945	maechler	1273	`calct(&h_s_0);`
		1274	`nexth(h_s_0);`
61750	ripley	1275	`}`
		1276	`omp = infin;`
		1277	`goto L10;`
		1278	`}`
		1279	`else {`
		1280
		1281	`/* exit if polynomial value */`
		1282	`/* increases significantly. */`
		1283
		1284	`if (mp * .1 > omp)`
87891	ripley	1285	`return false;`
61750	ripley	1286	`}`
		1287	`}`
		1288	`omp = mp;`
		1289
		1290	`/* calculate next iterate. */`
		1291
		1292	`L10:`
82945	maechler	1293	`calct(&h_s_0);`
		1294	`nexth(h_s_0);`
		1295	`calct(&h_s_0);`
		1296	`if (!h_s_0) {`
61750	ripley	1297	`relstp = hypot(tr, ti) / hypot(sr, si);`
		1298	`sr += tr;`
		1299	`si += ti;`
		1300	`}`
		1301	`}`
87891	ripley	1302	`return false;`
61750	ripley	1303
		1304	`L_conv:`
		1305	`*zr = sr;`
		1306	`*zi = si;`
87891	ripley	1307	`return true;`
61750	ripley	1308	`}`
		1309
87891	ripley	1310	`static void calct(bool *h_s_0)`
61750	ripley	1311	`{`
		1312	`/* computes t = -p(s)/h(s).`
82945	maechler	1313	`* h_s_0 - logical, set true if h(s) is essentially zero. */`
61750	ripley	1314
		1315	`int n = nn - 1;`
		1316	`double hvi, hvr;`
		1317
		1318	`/* evaluate h(s). */`
		1319	`polyev(n, sr, si, hr, hi,`
		1320	`qhr, qhi, &hvr, &hvi);`
		1321
82945	maechler	1322	`h_s_0 = hypot(hvr, hvi) <= are 10. * hypot(hr[n-1], hi[n-1]);`
		1323	`if (!*h_s_0) {`
61750	ripley	1324	`cdivid(-pvr, -pvi, hvr, hvi, &tr, &ti);`
		1325	`}`
		1326	`else {`
		1327	`tr = 0.;`
		1328	`ti = 0.;`
		1329	`}`
		1330	`}`
		1331
87891	ripley	1332	`static void nexth(bool h_s_0)`
61750	ripley	1333	`{`
		1334	`/* calculates the next shifted h polynomial.`
82945	maechler	1335	`* h_s_0 : if TRUE h(s) is essentially zero`
61750	ripley	1336	`*/`
		1337	`int j, n = nn - 1;`
		1338	`double t1, t2;`
		1339
82945	maechler	1340	`if (!h_s_0) {`
61750	ripley	1341	`for (j=1; j < n; j++) {`
		1342	`t1 = qhr[j - 1];`
		1343	`t2 = qhi[j - 1];`
		1344	`hr[j] = tr * t1 - ti * t2 + qpr[j];`
		1345	`hi[j] = tr * t2 + ti * t1 + qpi[j];`
		1346	`}`
		1347	`hr[0] = qpr[0];`
		1348	`hi[0] = qpi[0];`
		1349	`}`
		1350	`else {`
		1351	`/* if h(s) is zero replace h with qh. */`
		1352
		1353	`for (j=1; j < n; j++) {`
		1354	`hr[j] = qhr[j-1];`
		1355	`hi[j] = qhi[j-1];`
		1356	`}`
		1357	`hr[0] = 0.;`
		1358	`hi[0] = 0.;`
		1359	`}`
		1360	`}`
		1361
		1362	`/--------------------- Independent Complex Polynomial Utilities ----------/`
		1363
		1364	`static`
		1365	`void polyev(int n,`
		1366	`double s_r, double s_i,`
		1367	`double p_r, double p_i,`
		1368	`double q_r, double q_i,`
		1369	`double v_r, double v_i)`
		1370	`{`
		1371	`/* evaluates a polynomial p at s by the horner recurrence`
		1372	`* placing the partial sums in q and the computed value in v_.`
		1373	`*/`
		1374	`int i;`
		1375	`double t;`
		1376
		1377	`q_r[0] = p_r[0];`
		1378	`q_i[0] = p_i[0];`
		1379	`*v_r = q_r[0];`
		1380	`*v_i = q_i[0];`
		1381	`for (i = 1; i < n; i++) {`
		1382	`t = v_r s_r - v_i s_i + p_r[i];`
		1383	`q_i[i] = v_i = v_r * s_i + v_i s_r + p_i[i];`
		1384	`q_r[i] = *v_r = t;`
		1385	`}`
		1386	`}`
		1387
		1388	`static`
		1389	`double errev(int n, double qr, double qi,`
		1390	`double ms, double mp, double a_re, double m_re)`
		1391	`{`
		1392	`/* bounds the error in evaluating the polynomial by the horner`
		1393	`* recurrence.`
		1394	`*`
		1395	`* qr,qi - the partial sum vectors`
		1396	`* ms - modulus of the point`
		1397	`* mp - modulus of polynomial value`
		1398	`* a_re,m_re - error bounds on complex addition and multiplication`
		1399	`*/`
		1400	`double e;`
		1401	`int i;`
		1402
		1403	`e = hypot(qr[0], qi[0]) * m_re / (a_re + m_re);`
		1404	`for (i=0; i < n; i++)`
		1405	`e = e*ms + hypot(qr[i], qi[i]);`
		1406
		1407	`return e * (a_re + m_re) - mp * m_re;`
		1408	`}`
		1409
		1410
		1411	`static`
		1412	`double cpoly_cauchy(int n, double pot, double q)`
		1413	`{`
		1414	`/* Computes a lower bound on the moduli of the zeros of a polynomial`
		1415	`* pot[1:nn] is the modulus of the coefficients.`
		1416	`*/`
		1417	`double f, x, delf, dx, xm;`
		1418	`int i, n1 = n - 1;`
		1419
		1420	`pot[n1] = -pot[n1];`
		1421
		1422	`/* compute upper estimate of bound. */`
		1423
		1424	`x = exp((log(-pot[n1]) - log(pot[0])) / (double) n1);`
		1425
		1426	`/* if newton step at the origin is better, use it. */`
		1427
		1428	`if (pot[n1-1] != 0.) {`
		1429	`xm = -pot[n1] / pot[n1-1];`
		1430	`if (xm < x)`
		1431	`x = xm;`
		1432	`}`
		1433
		1434	`/* chop the interval (0,x) unitl f le 0. */`
		1435
		1436	`for(;;) {`
		1437	`xm = x * 0.1;`
		1438	`f = pot[0];`
		1439	`for (i = 1; i < n; i++)`
		1440	`f = f * xm + pot[i];`
		1441	`if (f <= 0.0) {`
		1442	`break;`
		1443	`}`
		1444	`x = xm;`
		1445	`}`
		1446
		1447	`dx = x;`
		1448
		1449	`/* do Newton iteration until x converges to two decimal places. */`
		1450
		1451	`while (fabs(dx / x) > 0.005) {`
		1452	`q[0] = pot[0];`
		1453	`for(i = 1; i < n; i++)`
		1454	`q[i] = q[i-1] * x + pot[i];`
		1455	`f = q[n1];`
		1456	`delf = q[0];`
		1457	`for(i = 1; i < n1; i++)`
		1458	`delf = delf * x + q[i];`
		1459	`dx = f / delf;`
		1460	`x -= dx;`
		1461	`}`
		1462	`return x;`
		1463	`}`
		1464
		1465	`static`
		1466	`double cpoly_scale(int n, double *pot,`
		1467	`double eps, double BIG, double small, double base)`
		1468	`{`
		1469	`/* Returns a scale factor to multiply the coefficients of the polynomial.`
		1470	`* The scaling is done to avoid overflow and to avoid`
		1471	`* undetected underflow interfering with the convergence criterion.`
		1472	`* The factor is a power of the base.`
		1473
		1474	`* pot[1:n] : modulus of coefficients of p`
		1475	`* eps,BIG,`
		1476	`* small,base - constants describing the floating point arithmetic.`
		1477	`*/`
		1478
		1479	`int i, ell;`
		1480	`double x, high, sc, lo, min_, max_;`
		1481
		1482	`/* find largest and smallest moduli of coefficients. */`
		1483	`high = sqrt(BIG);`
		1484	`lo = small / eps;`
		1485	`max_ = 0.;`
		1486	`min_ = BIG;`
		1487	`for (i = 0; i < n; i++) {`
		1488	`x = pot[i];`
		1489	`if (x > max_) max_ = x;`
		1490	`if (x != 0. && x < min_)`
		1491	`min_ = x;`
		1492	`}`
		1493
		1494	`/* scale only if there are very large or very small components. */`
		1495
		1496	`if (min_ < lo \|\| max_ > high) {`
		1497	`x = lo / min_;`
		1498	`if (x <= 1.)`
		1499	`sc = 1. / (sqrt(max_) * sqrt(min_));`
		1500	`else {`
		1501	`sc = x;`
		1502	`if (BIG / sc > max_)`
		1503	`sc = 1.0;`
		1504	`}`
		1505	`ell = (int) (log(sc) / log(base) + 0.5);`
		1506	`return R_pow_di(base, ell);`
		1507	`}`
		1508	`else return 1.0;`
		1509	`}`
		1510
		1511
		1512	`static`
		1513	`void cdivid(double ar, double ai, double br, double bi,`
		1514	`double cr, double ci)`
		1515	`{`
		1516	`/* complex division c = a/b, i.e., (cr +ici) = (ar +iai) / (br +i*bi),`
		1517	`avoiding overflow. */`
		1518
		1519	`double d, r;`
		1520
		1521	`if (br == 0. && bi == 0.) {`
		1522	`/* division by zero, c = infinity. */`
		1523	`cr = ci = R_PosInf;`
		1524	`}`
		1525	`else if (fabs(br) >= fabs(bi)) {`
		1526	`r = bi / br;`
		1527	`d = br + r * bi;`
		1528	`cr = (ar + ai r) / d;`
		1529	`ci = (ai - ar r) / d;`
		1530	`}`
		1531	`else {`
		1532	`r = br / bi;`
		1533	`d = bi + r * br;`
		1534	`cr = (ar r + ai) / d;`
		1535	`ci = (ai r - ar) / d;`
		1536	`}`
		1537	`}`
		1538
		1539	`/* static double cpoly_cmod(double r, double i)`
		1540	`* --> replaced by hypot() everywhere`
		1541	`*/`

The R Project SVN R

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