1 /*
2 * R : A Computer Language for Statistical Data Analysis
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
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
18 * along with this program; if not, a copy is available at
19 * https://www.R-project.org/Licenses/
20 */
21
22 #ifdef HAVE_CONFIG_H
23 #include <config.h>
24 #endif
25
26 /* Note: gcc -pedantic may warn in several places about C99 features
27 as extensions.
28 This was a very-long-standing GCC bug, http://gcc.gnu.org/PR7263
29 The system <complex.h> header can work around it: some do.
30 It should have been resolved (after a decade) in 2012.
31 */
32
33 #if defined(HAVE_CTANH) && !defined(HAVE_WORKING_CTANH)
34 #undef HAVE_CTANH
35 #endif
36
37 #if 0
38 /* For testing substitute fns */
39 #undef HAVE_CARG
40 #undef HAVE_CABS
41 #undef HAVE_CPOW
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
56 #ifdef __SUNPRO_C
57 /* segfaults in Solaris Studio 12.3 */
58 #undef HAVE_CPOW
59 #endif
60
61 #include <Defn.h> /* -> ../include/R_ext/Complex.h */
62 #include <Internal.h>
63 #include <Rmath.h>
64
65 #include "arithmetic.h" /* complex_* */
66 #include <complex.h>
67 #include "Rcomplex.h" /* I, SET_C99_COMPLEX, toC99 */
68 #include <R_ext/Itermacros.h>
69
70
71 /* interval at which to check interrupts, a guess */
72 #define NINTERRUPT 10000000
73
74
complex_unary(ARITHOP_TYPE code,SEXP s1,SEXP call)75 SEXP attribute_hidden complex_unary(ARITHOP_TYPE code, SEXP s1, SEXP call)
76 {
77 R_xlen_t i, n;
78 SEXP ans;
79
80 switch(code) {
81 case PLUSOP:
82 return s1;
83 case MINUSOP:
84 ans = NO_REFERENCES(s1) ? s1 : duplicate(s1);
85 Rcomplex *pans = COMPLEX(ans);
86 const Rcomplex *ps1 = COMPLEX_RO(s1);
87 n = XLENGTH(s1);
88 for (i = 0; i < n; i++) {
89 Rcomplex x = ps1[i];
90 pans[i].r = -x.r;
91 pans[i].i = -x.i;
92 }
93 return ans;
94 default:
95 errorcall(call, _("invalid complex unary operator"));
96 }
97 return R_NilValue; /* -Wall */
98 }
99
R_cpow_n(double complex X,int k)100 static R_INLINE double complex R_cpow_n(double complex X, int k)
101 {
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);
105 else {/* k > 0 */
106 double complex z = (double complex) 1.;;
107 while (k > 0) {
108 if (k & 1) z = z * X;
109 if (k == 1) break;
110 k >>= 1; /* efficient division by 2; now have k >= 1 */
111 X = X * X;
112 }
113 return z;
114 }
115 }
116
117 #if defined(Win32)
118 # undef HAVE_CPOW
119 #endif
120 /* reason for this:
121 1) X^n (e.g. for n = +/- 2, 3) is unnecessarily inaccurate in glibc;
122 cut-off 65536 : guided from empirical speed measurements
123
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.
126
127 3) PPC macOS crashed on powers of 0+0i (at least under Rosetta).
128 Really 0i^-1 should by Inf+NaNi, but getting that portably seems too hard.
129 (C1x's CMPLX will eventually be possible.)
130 */
131
mycpow(double complex X,double complex Y)132 static double complex mycpow (double complex X, double complex Y)
133 {
134 double complex Z;
135 double yr = creal(Y), yi = cimag(Y);
136 int k;
137 if (X == 0.0) {
138 if (yi == 0.0) Z = R_pow(0.0, yr); else Z = R_NaN + R_NaN*I;
139 } else if (yi == 0.0 && yr == (k = (int) yr) && abs(k) <= 65536)
140 Z = R_cpow_n(X, k);
141 else
142 #if defined(HAVE_CPOW) && !defined(__FreeBSD__)
143 Z = cpow(X, Y);
144 #else
145 {
146 /* Used for FreeBSD and MingGW, hence mainly with gcc */
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)) */
155 r = log(r);
156 theta += r * yi;
157 rho = exp(r * yr - i * yi);
158 }
159 #ifdef __GNUC__
160 __real__ Z = rho * cos(theta);
161 __imag__ Z = rho * sin(theta);
162 #else
163 Z = rho * cos(theta) + (rho * sin(theta)) * I;
164 #endif
165 }
166 #endif
167 return Z;
168 }
169
170
171
complex_binary(ARITHOP_TYPE code,SEXP s1,SEXP s2)172 SEXP attribute_hidden complex_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2)
173 {
174 R_xlen_t i, i1, i2, n, n1, n2;
175 SEXP ans;
176
177 /* Note: "s1" and "s2" are protected in the calling code. */
178 n1 = XLENGTH(s1);
179 n2 = XLENGTH(s2);
180 /* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */
181 if (n1 == 0 || n2 == 0) return(allocVector(CPLXSXP, 0));
182
183 n = (n1 > n2) ? n1 : n2;
184 ans = R_allocOrReuseVector(s1, s2, CPLXSXP, n);
185 PROTECT(ans);
186
187 Rcomplex *pans = COMPLEX(ans);
188 const Rcomplex *ps1 = COMPLEX_RO(s1);
189 const Rcomplex *ps2 = COMPLEX_RO(s2);
190
191 switch (code) {
192 case PLUSOP:
193 MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {
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;
198 });
199 break;
200 case MINUSOP:
201 MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {
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;
206 });
207 break;
208 case TIMESOP:
209 MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {
210 SET_C99_COMPLEX(pans, i,
211 toC99(&ps1[i1]) * toC99(&ps2[i2]));
212 });
213 break;
214 case DIVOP:
215 MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {
216 SET_C99_COMPLEX(pans, i,
217 toC99(&ps1[i1]) / toC99(&ps2[i2]));
218 });
219 break;
220 case POWOP:
221 MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {
222 SET_C99_COMPLEX(pans, i,
223 mycpow(toC99(&ps1[i1]), toC99(&ps2[i2])));
224 });
225 break;
226 default:
227 error(_("unimplemented complex operation"));
228 }
229 UNPROTECT(1);
230
231 /* quick return if there are no attributes */
232 if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue)
233 return ans;
234
235 /* Copy attributes from longer argument. */
236
237 if (ans != s2 && n == n2 && ATTRIB(s2) != R_NilValue)
238 copyMostAttrib(s2, ans);
239 if (ans != s1 && n == n1 && ATTRIB(s1) != R_NilValue)
240 copyMostAttrib(s1, ans); /* Done 2nd so s1's attrs overwrite s2's */
241
242 return ans;
243 }
244
do_cmathfuns(SEXP call,SEXP op,SEXP args,SEXP env)245 SEXP attribute_hidden do_cmathfuns(SEXP call, SEXP op, SEXP args, SEXP env)
246 {
247 SEXP x, y = R_NilValue; /* -Wall*/
248 R_xlen_t i, n;
249
250 checkArity(op, args);
251 check1arg(args, call, "z");
252 if (DispatchGroup("Complex", call, op, args, env, &x))
253 return x;
254 x = CAR(args);
255 if (isComplex(x)) {
256 n = XLENGTH(x);
257 const Rcomplex *px = COMPLEX_RO(x);
258 switch(PRIMVAL(op)) {
259 case 1: /* Re */
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 }
266 break;
267 case 2: /* Im */
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 }
274 break;
275 case 3: /* Mod */
276 case 6: /* abs */
277 {
278 y = allocVector(REALSXP, n);
279 double *py = REAL(y);
280 for(i = 0 ; i < n ; i++)
281 #if HAVE_CABS
282 py[i] = cabs(toC99(&px[i]));
283 #else
284 py[i] = hypot(px[i].r, px[i].i);
285 #endif
286 }
287 break;
288 case 4: /* Arg */
289 {
290 y = allocVector(REALSXP, n);
291 double *py = REAL(y);
292 for(i = 0 ; i < n ; i++)
293 #if HAVE_CARG
294 py[i] = carg(toC99(&px[i]));
295 #else
296 py[i] = atan2(px[i].i, px[i].r);
297 #endif
298 }
299 break;
300 case 5: /* Conj */
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 }
308 }
309 break;
310 }
311 }
312 else if(isNumeric(x)) { /* so no complex numbers involved */
313 n = XLENGTH(x);
314 if(isReal(x)) PROTECT(x);
315 else PROTECT(x = coerceVector(x, REALSXP));
316 y = NO_REFERENCES(x) ? x : allocVector(REALSXP, n);
317 double *py = REAL(y);
318 const double *px = REAL_RO(x);
319
320 switch(PRIMVAL(op)) {
321 case 1: /* Re */
322 case 5: /* Conj */
323 for(i = 0 ; i < n ; i++)
324 py[i] = px[i];
325 break;
326 case 2: /* Im */
327 for(i = 0 ; i < n ; i++)
328 py[i] = 0.0;
329 break;
330 case 4: /* Arg */
331 for(i = 0 ; i < n ; i++)
332 if(ISNAN(px[i]))
333 py[i] = px[i];
334 else if (px[i] >= 0)
335 py[i] = 0;
336 else
337 py[i] = M_PI;
338 break;
339 case 3: /* Mod */
340 case 6: /* abs */
341 for(i = 0 ; i < n ; i++)
342 py[i] = fabs(px[i]);
343 break;
344 }
345 UNPROTECT(1);
346 }
347 else errorcall(call, _("non-numeric argument to function"));
348
349 if (x != y && ATTRIB(x) != R_NilValue) {
350 PROTECT(x);
351 PROTECT(y);
352 SHALLOW_DUPLICATE_ATTRIB(y, x);
353 UNPROTECT(2);
354 }
355 return y;
356 }
357
358 /* Implementing signif(<complex>) *and* used in format.c and printutils.c */
359 #define MAX_DIGITS 22
z_prec_r(Rcomplex * r,const Rcomplex * x,double digits)360 void attribute_hidden z_prec_r(Rcomplex *r, const Rcomplex *x, double digits)
361 {
362 // Implement r <- signif(x, digits)
363
364 r->r = x->r; r->i = x->i;
365 double m = 0.0,
366 m1 = fabs(x->r),
367 m2 = fabs(x->i);
368 if(R_FINITE(m1)) m = m1;
369 if(R_FINITE(m2) && m2 > m) m = m2;
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 }
374 int dig = (int)floor(digits+0.5);
375 if (dig > MAX_DIGITS) return; else if (dig < 1) dig = 1;
376 int mag = (int)floor(log10(m));
377 dig = dig - mag - 1;
378 if (dig > 306) {
379 double pow10 = 1.0e4;
380 digits = (double)(dig - 4);
381 r->r = fround(pow10 * x->r, digits)/pow10;
382 r->i = fround(pow10 * x->i, digits)/pow10;
383 } else {
384 digits = (double)(dig);
385 r->r = fround(x->r, digits);
386 r->i = fround(x->i, digits);
387 }
388 }
389
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!
393
394 Currently (Feb 2011) they are used on FreeBSD.
395 */
396
397 #if !defined(HAVE_CLOG) || defined(__FreeBSD__)
398 #define clog R_clog
399 /* FIXME: maybe add full IEC60559 support */
clog(double complex x)400 static double complex clog(double complex x)
401 {
402 double xr = creal(x), xi = cimag(x);
403 return log(hypot(xr, xi)) + atan2(xi, xr)*I;
404 }
405 #endif
406
407 #ifndef HAVE_CSQRT
408 #define csqrt R_csqrt
409 /* FreeBSD does have this one */
csqrt(double complex x)410 static double complex csqrt(double complex x)
411 {
412 return mycpow(x, 0.5+0.0*I);
413 }
414 #endif
415
416 #ifndef HAVE_CEXP
417 #define cexp R_cexp
418 /* FIXME: check/add full IEC60559 support */
cexp(double complex x)419 static double complex cexp(double complex x)
420 {
421 double expx = exp(creal(x)), y = cimag(x);
422 return expx * cos(y) + (expx * sin(y)) * I;
423 }
424 #endif
425
426 #ifndef HAVE_CCOS
427 #define ccos R_ccos
ccos(double complex x)428 static double complex ccos(double complex x)
429 {
430 double xr = creal(x), xi = cimag(x);
431 return cos(xr)*cosh(xi) - sin(xr)*sinh(xi)*I; /* A&S 4.3.56 */
432 }
433 #endif
434
435 #ifndef HAVE_CSIN
436 #define csin R_csin
csin(double complex x)437 static double complex csin(double complex x)
438 {
439 double xr = creal(x), xi = cimag(x);
440 return sin(xr)*cosh(xi) + cos(xr)*sinh(xi)*I; /* A&S 4.3.55 */
441 }
442 #endif
443
444 #ifndef HAVE_CTAN
445 #define ctan R_ctan
ctan(double complex z)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
460 #ifndef HAVE_CASIN
461 #define casin R_casin
casin(double complex z)462 static double complex casin(double complex z)
463 {
464 /* A&S 4.4.37 */
465 double alpha, t1, t2, x = creal(z), y = cimag(z), ri;
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));
470 /* This comes from
471 'z_asin() is continuous from below if x >= 1
472 and continuous from above if x <= -1.'
473 */
474 if(y < 0 || (y == 0 && x > 1)) ri *= -1;
475 return asin(t1 - t2) + ri*I;
476 }
477 #endif
478
479 #ifndef HAVE_CACOS
480 #define cacos R_cacos
cacos(double complex z)481 static double complex cacos(double complex z)
482 {
483 return M_PI_2 - casin(z);
484 }
485 #endif
486
487 #ifndef HAVE_CATAN
488 #define catan R_catan
catan(double complex z)489 static double complex catan(double complex z)
490 {
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)));
495 return rr + ri*I;
496 }
497 #endif
498
499 #ifndef HAVE_CCOSH
500 #define ccosh R_ccosh
ccosh(double complex z)501 static double complex ccosh(double complex z)
502 {
503 return ccos(z * I); /* A&S 4.5.8 */
504 }
505 #endif
506
507 #ifndef HAVE_CSINH
508 #define csinh R_csinh
csinh(double complex z)509 static double complex csinh(double complex z)
510 {
511 return -I * csin(z * I); /* A&S 4.5.7 */
512 }
513 #endif
514
z_tan(double complex z)515 static double complex z_tan(double complex z)
516 {
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 }
528 return r;
529 }
530
531 #ifndef HAVE_CTANH
532 #define ctanh R_ctanh
ctanh(double complex z)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
540 /* Don't rely on the OS at the branch cuts */
541
z_asin(double complex z)542 static double complex z_asin(double complex z)
543 {
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);
554 }
555
z_acos(double complex z)556 static double complex z_acos(double complex z)
557 {
558 if(cimag(z) == 0 && fabs(creal(z)) > 1) return M_PI_2 - z_asin(z);
559 return cacos(z);
560 }
561
z_atan(double complex z)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 }
572
z_acosh(double complex z)573 static double complex z_acosh(double complex z)
574 {
575 return z_acos(z) * I;
576 }
577
z_asinh(double complex z)578 static double complex z_asinh(double complex z)
579 {
580 return -I * z_asin(z * I);
581 }
582
z_atanh(double complex z)583 static double complex z_atanh(double complex z)
584 {
585 return -I * z_atan(z * I);
586 }
587
cmath1(double complex (* f)(double complex),const Rcomplex * x,Rcomplex * y,R_xlen_t n)588 static Rboolean cmath1(double complex (*f)(double complex),
589 const Rcomplex *x, Rcomplex *y, R_xlen_t n)
590 {
591 R_xlen_t i;
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)));
598 if ( (ISNAN(y[i].r) || ISNAN(y[i].i)) &&
599 !(ISNAN(x[i].r) || ISNAN(x[i].i)) ) naflag = TRUE;
600 }
601 }
602 return naflag;
603 }
604
complex_math1(SEXP call,SEXP op,SEXP args,SEXP env)605 SEXP attribute_hidden complex_math1(SEXP call, SEXP op, SEXP args, SEXP env)
606 {
607 SEXP x, y;
608 R_xlen_t n;
609 Rboolean naflag = FALSE;
610
611 PROTECT(x = CAR(args));
612 n = XLENGTH(x);
613 PROTECT(y = allocVector(CPLXSXP, n));
614
615 const Rcomplex *px = COMPLEX_RO(x);
616 Rcomplex *py = COMPLEX(y);
617
618 switch (PRIMVAL(op)) {
619 case 10003: naflag = cmath1(clog, px, py, n); break;
620 // 1: floor
621 // 2: ceil[ing]
622 case 3: naflag = cmath1(csqrt, px, py, n); break;
623 // 4: sign
624 case 10: naflag = cmath1(cexp, px, py, n); break;
625 // 11: expm1
626 // 12: log1p
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;
633
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;
640
641 default:
642 /* such as sign, gamma */
643 errorcall(call, _("unimplemented complex function"));
644 }
645 if (naflag)
646 warningcall(call, "NaNs produced in function \"%s\"", PRIMNAME(op));
647 SHALLOW_DUPLICATE_ATTRIB(y, x);
648 UNPROTECT(2);
649 return y;
650 }
651
z_rround(Rcomplex * r,Rcomplex * x,Rcomplex * p)652 static void z_rround(Rcomplex *r, Rcomplex *x, Rcomplex *p)
653 {
654 r->r = fround(x->r, p->r);
655 r->i = fround(x->i, p->r);
656 }
657
z_prec(Rcomplex * r,Rcomplex * x,Rcomplex * p)658 static void z_prec(Rcomplex *r, Rcomplex *x, Rcomplex *p)
659 {
660 z_prec_r(r, x, p->r);
661 }
662
z_logbase(Rcomplex * r,Rcomplex * z,Rcomplex * base)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
z_atan2(Rcomplex * r,Rcomplex * csn,Rcomplex * ccs)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) {
674 r->r = NA_REAL; r->i = NA_REAL; /* Why not R_NaN? */
675 return;
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 }
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
689
690 /* Complex Functions of Two Arguments */
691
692 typedef void (*cm2_fun)(Rcomplex *, Rcomplex *, Rcomplex *);
complex_math2(SEXP call,SEXP op,SEXP args,SEXP env)693 SEXP attribute_hidden complex_math2(SEXP call, SEXP op, SEXP args, SEXP env)
694 {
695 R_xlen_t i, n, na, nb, ia, ib;
696 Rcomplex ai, bi, *y;
697 const Rcomplex *a, *b;
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:
714 error_return(_("unimplemented complex function"));
715 }
716
717 PROTECT(sa = coerceVector(CAR(args), CPLXSXP));
718 PROTECT(sb = coerceVector(CADR(args), CPLXSXP));
719 na = XLENGTH(sa); nb = XLENGTH(sb);
720 if ((na == 0) || (nb == 0)) {
721 UNPROTECT(2);
722 return(allocVector(CPLXSXP, 0));
723 }
724 n = (na < nb) ? nb : na;
725 PROTECT(sy = allocVector(CPLXSXP, n));
726 a = COMPLEX_RO(sa); b = COMPLEX_RO(sb);
727 y = COMPLEX(sy);
728 MOD_ITERATE2(n, na, nb, i, ia, ib, {
729 ai = a[ia]; bi = b[ib];
730 if(ISNA(ai.r) && ISNA(ai.i) &&
731 ISNA(bi.r) && ISNA(bi.i)) {
732 y[i].r = NA_REAL; y[i].i = NA_REAL;
733 } else {
734 f(&y[i], &ai, &bi);
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;
738 }
739 });
740 if (naflag)
741 warning("NaNs produced in function \"%s\"", PRIMNAME(op));
742 if(n == na) {
743 SHALLOW_DUPLICATE_ATTRIB(sy, sa);
744 } else if(n == nb) {
745 SHALLOW_DUPLICATE_ATTRIB(sy, sb);
746 }
747 UNPROTECT(3);
748 return sy;
749 }
750
do_complex(SEXP call,SEXP op,SEXP args,SEXP rho)751 SEXP attribute_hidden do_complex(SEXP call, SEXP op, SEXP args, SEXP rho)
752 {
753 /* complex(length, real, imaginary) */
754 SEXP ans, re, im;
755 R_xlen_t i, na, nr, ni;
756
757 checkArity(op, args);
758 na = asInteger(CAR(args));
759 if(na == NA_INTEGER || na < 0)
760 error(_("invalid length"));
761 PROTECT(re = coerceVector(CADR(args), REALSXP));
762 PROTECT(im = coerceVector(CADDR(args), REALSXP));
763 nr = XLENGTH(re);
764 ni = XLENGTH(im);
765 /* is always true: if (na >= 0) {*/
766 na = (nr > na) ? nr : na;
767 na = (ni > na) ? ni : na;
768 /* }*/
769 ans = allocVector(CPLXSXP, na);
770 Rcomplex *pans = COMPLEX(ans);
771 for(i=0 ; i<na ; i++) {
772 pans[i].r = 0;
773 pans[i].i = 0;
774 }
775 UNPROTECT(2);
776 if(na > 0 && nr > 0) {
777 const double *p_re = REAL_RO(re);
778 for(i=0 ; i<na ; i++)
779 pans[i].r = p_re[i%nr];
780 }
781 if(na > 0 && ni > 0) {
782 const double *p_im = REAL_RO(im);
783 for(i=0 ; i<na ; i++)
784 pans[i].i = p_im[i%ni];
785 }
786 return ans;
787 }
788
789 static void R_cpolyroot(double *opr, double *opi, int *degree,
790 double *zeror, double *zeroi, Rboolean *fail);
791
do_polyroot(SEXP call,SEXP op,SEXP args,SEXP rho)792 SEXP attribute_hidden do_polyroot(SEXP call, SEXP op, SEXP args, SEXP rho)
793 {
794 SEXP z, zr, zi, r, rr, ri;
795 Rboolean fail;
796 int degree, i, n;
797
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:
810 UNIMPLEMENTED_TYPE("polyroot", z);
811 }
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
819 const Rcomplex *pz = COMPLEX_RO(z);
820 degree = 0;
821 for(i = 0; i < n; i++) {
822 if(pz[i].r!= 0.0 || pz[i].i != 0.0) degree = i;
823 }
824 n = degree + 1; /* omit trailing zeroes */
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));
830
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
836 for(i = 0 ; i < n ; i++) {
837 if(!R_FINITE(pz[i].r) || !R_FINITE(pz[i].i))
838 error(_("invalid polynomial coefficient"));
839 p_zr[degree-i] = pz[i].r;
840 p_zi[degree-i] = pz[i].i;
841 }
842 R_cpolyroot(p_zr, p_zi, °ree, p_rr, p_ri, &fail);
843 if(fail) error(_("root finding code failed"));
844 UNPROTECT(2);
845 r = allocVector(CPLXSXP, degree);
846 Rcomplex *pr = COMPLEX(r);
847 for(i = 0 ; i < degree ; i++) {
848 pr[i].r = p_rr[i];
849 pr[i].i = p_ri[i];
850 }
851 UNPROTECT(3);
852 }
853 else {
854 UNPROTECT(1);
855 r = allocVector(CPLXSXP, 0);
856 }
857 return r;
858 }
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
R_cpolyroot(double * opr,double * opi,int * degree,double * zeror,double * zeroi,Rboolean * fail)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 */
975 const void *vmax = vmaxget();
976 tmp = (double *) R_alloc((size_t) (10*nn), sizeof(double));
977 pr = tmp; pi = tmp + nn; hr = tmp + 2*nn; hi = tmp + 3*nn;
978 qpr = tmp + 4*nn; qpi = tmp + 5*nn; qhr = tmp + 6*nn; qhi = tmp + 7*nn;
979 shr = tmp + 8*nn; shi = tmp + 9*nn;
980
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;
1041 vmaxset(vmax);
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
1061 vmaxset(vmax);
1062 return;
1063 }
1064
1065
1066 /* Computes the derivative polynomial as the initial
1067 * polynomial and computes l1 no-shift h polynomials. */
1068
noshft(int l1)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 */
fxshft(int l2,double * zr,double * zi)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 */
vrshft(int l3,double * zr,double * zi)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
calct(Rboolean * bool)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
nexth(Rboolean bool)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
polyev(int n,double s_r,double s_i,double * p_r,double * p_i,double * q_r,double * q_i,double * v_r,double * v_i)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
errev(int n,double * qr,double * qi,double ms,double mp,double a_re,double m_re)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
cpoly_cauchy(int n,double * pot,double * q)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
cpoly_scale(int n,double * pot,double eps,double BIG,double small,double base)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
cdivid(double ar,double ai,double br,double bi,double * cr,double * ci)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 +i*ci) = (ar +i*ai) / (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 */
1535