1 /* Formerly src/appl/cpoly.c:
2 *
3 * Copyright (C) 1997-1998 Ross Ihaka
4 * Copyright (C) 1999-2001 R Core Team
5 *
6 * cpoly finds the zeros of a complex polynomial.
7 *
8 * On Entry
9 *
10 * opr, opi - double precision vectors of real and
11 * imaginary parts of the coefficients in
12 * order of decreasing powers.
13 *
14 * degree - int degree of polynomial.
15 *
16 *
17 * On Return
18 *
19 * zeror, zeroi - output double precision vectors of
20 * real and imaginary parts of the zeros.
21 *
22 * fail - output int parameter, true only if
23 * leading coefficient is zero or if cpoly
24 * has found fewer than degree zeros.
25 *
26 * The program has been written to reduce the chance of overflow
27 * occurring. If it does occur, there is still a possibility that
28 * the zerofinder will work provided the overflowed quantity is
29 * replaced by a large number.
30 *
31 * This is a C translation of the following.
32 *
33 * TOMS Algorithm 419
34 * Jenkins and Traub.
35 * Comm. ACM 15 (1972) 97-99.
36 *
37 * Ross Ihaka
38 * February 1997
39 */
40
41 #include <R_ext/Arith.h> /* for declaration of hypot */
42 #include <R_ext/Memory.h> /* for declaration of R_alloc */
43
44 #include <float.h> /* for FLT_RADIX */
45
46 #include <Rmath.h> /* for R_pow_di */
47
48 static void calct(Rboolean *);
49 static Rboolean fxshft(int, double *, double *);
50 static Rboolean vrshft(int, double *, double *);
51 static void nexth(Rboolean);
52 static void noshft(int);
53
54 static void polyev(int, double, double,
55 double *, double *, double *, double *, double *, double *);
56 static double errev(int, double *, double *, double, double, double, double);
57 static double cpoly_cauchy(int, double *, double *);
58 static double cpoly_scale(int, double *, double, double, double, double);
59 static void cdivid(double, double, double, double, double *, double *);
60
61 /* Global Variables (too many!) */
62
63 static int nn;
64 static double *pr, *pi, *hr, *hi, *qpr, *qpi, *qhr, *qhi, *shr, *shi;
65 static double sr, si;
66 static double tr, ti;
67 static double pvr, pvi;
68
69 static const double eta = DBL_EPSILON;
70 static const double are = /* eta = */DBL_EPSILON;
71 static const double mre = 2. * M_SQRT2 * /* eta, i.e. */DBL_EPSILON;
72 static const double infin = DBL_MAX;
73
cpolyroot(double * opr,double * opi,int * degree,double * zeror,double * zeroi,Rboolean * fail)74 void cpolyroot(double *opr, double *opi, int *degree,
75 double *zeror, double *zeroi, Rboolean *fail)
76 {
77 static const double smalno = DBL_MIN;
78 static const double base = (double)FLT_RADIX;
79 static int d_n, i, i1, i2;
80 static double zi, zr, xx, yy;
81 static double bnd, xxx;
82 Rboolean conv;
83 int d1;
84 double *tmp;
85 static const double cosr =/* cos 94 */ -0.06975647374412529990;
86 static const double sinr =/* sin 94 */ 0.99756405025982424767;
87 xx = M_SQRT1_2;/* 1/sqrt(2) = 0.707.... */
88
89 yy = -xx;
90 *fail = FALSE;
91
92 nn = *degree;
93 d1 = nn - 1;
94
95 /* algorithm fails if the leading coefficient is zero. */
96
97 if (opr[0] == 0. && opi[0] == 0.) {
98 *fail = TRUE;
99 return;
100 }
101
102 /* remove the zeros at the origin if any. */
103
104 while (opr[nn] == 0. && opi[nn] == 0.) {
105 d_n = d1-nn+1;
106 zeror[d_n] = 0.;
107 zeroi[d_n] = 0.;
108 nn--;
109 }
110 nn++;
111 /*-- Now, global var. nn := #{coefficients} = (relevant degree)+1 */
112
113 if (nn == 1) return;
114
115 /* Use a single allocation as these as small */
116 tmp = (double *) R_alloc((size_t) (10*nn), sizeof(double));
117 pr = tmp; pi = tmp + nn; hr = tmp + 2*nn; hi = tmp + 3*nn;
118 qpr = tmp + 4*nn; qpi = tmp + 5*nn; qhr = tmp + 6*nn; qhi = tmp + 7*nn;
119 shr = tmp + 8*nn; shi = tmp + 9*nn;
120
121 /* make a copy of the coefficients and shr[] = | p[] | */
122 for (i = 0; i < nn; i++) {
123 pr[i] = opr[i];
124 pi[i] = opi[i];
125 shr[i] = hypot(pr[i], pi[i]);
126 }
127
128 /* scale the polynomial with factor 'bnd'. */
129 bnd = cpoly_scale(nn, shr, eta, infin, smalno, base);
130 if (bnd != 1.) {
131 for (i=0; i < nn; i++) {
132 pr[i] *= bnd;
133 pi[i] *= bnd;
134 }
135 }
136
137 /* start the algorithm for one zero */
138
139 while (nn > 2) {
140
141 /* calculate bnd, a lower bound on the modulus of the zeros. */
142
143 for (i=0 ; i < nn ; i++)
144 shr[i] = hypot(pr[i], pi[i]);
145 bnd = cpoly_cauchy(nn, shr, shi);
146
147 /* outer loop to control 2 major passes */
148 /* with different sequences of shifts */
149
150 for (i1 = 1; i1 <= 2; i1++) {
151
152 /* first stage calculation, no shift */
153
154 noshft(5);
155
156 /* inner loop to select a shift */
157 for (i2 = 1; i2 <= 9; i2++) {
158
159 /* shift is chosen with modulus bnd */
160 /* and amplitude rotated by 94 degrees */
161 /* from the previous shift */
162
163 xxx= cosr * xx - sinr * yy;
164 yy = sinr * xx + cosr * yy;
165 xx = xxx;
166 sr = bnd * xx;
167 si = bnd * yy;
168
169 /* second stage calculation, fixed shift */
170
171 conv = fxshft(i2 * 10, &zr, &zi);
172 if (conv)
173 goto L10;
174 }
175 }
176
177 /* the zerofinder has failed on two major passes */
178 /* return empty handed */
179
180 *fail = TRUE;
181 return;
182
183 /* the second stage jumps directly to the third stage iteration.
184 * if successful, the zero is stored and the polynomial deflated.
185 */
186 L10:
187 d_n = d1+2 - nn;
188 zeror[d_n] = zr;
189 zeroi[d_n] = zi;
190 --nn;
191 for (i=0; i < nn ; i++) {
192 pr[i] = qpr[i];
193 pi[i] = qpi[i];
194 }
195 }/*while*/
196
197 /* calculate the final zero and return */
198 cdivid(-pr[1], -pi[1], pr[0], pi[0], &zeror[d1], &zeroi[d1]);
199
200 return;
201 }
202
203
204 /* Computes the derivative polynomial as the initial
205 * polynomial and computes l1 no-shift h polynomials. */
206
noshft(int l1)207 static void noshft(int l1)
208 {
209 int i, j, jj, n = nn - 1, nm1 = n - 1;
210
211 double t1, t2, xni;
212
213 for (i=0; i < n; i++) {
214 xni = (double)(nn - i - 1);
215 hr[i] = xni * pr[i] / n;
216 hi[i] = xni * pi[i] / n;
217 }
218
219 for (jj = 1; jj <= l1; jj++) {
220
221 if (hypot(hr[n-1], hi[n-1]) <=
222 eta * 10.0 * hypot(pr[n-1], pi[n-1])) {
223 /* If the constant term is essentially zero, */
224 /* shift h coefficients. */
225
226 for (i = 1; i <= nm1; i++) {
227 j = nn - i;
228 hr[j-1] = hr[j-2];
229 hi[j-1] = hi[j-2];
230 }
231 hr[0] = 0.;
232 hi[0] = 0.;
233 }
234 else {
235 cdivid(-pr[nn-1], -pi[nn-1], hr[n-1], hi[n-1], &tr, &ti);
236 for (i = 1; i <= nm1; i++) {
237 j = nn - i;
238 t1 = hr[j-2];
239 t2 = hi[j-2];
240 hr[j-1] = tr * t1 - ti * t2 + pr[j-1];
241 hi[j-1] = tr * t2 + ti * t1 + pi[j-1];
242 }
243 hr[0] = pr[0];
244 hi[0] = pi[0];
245 }
246 }
247 }
248
249
250 /* Computes l2 fixed-shift h polynomials and tests for convergence.
251 * initiates a variable-shift iteration and returns with the
252 * approximate zero if successful.
253 */
fxshft(int l2,double * zr,double * zi)254 static Rboolean fxshft(int l2, double *zr, double *zi)
255 {
256 /* l2 - limit of fixed shift steps
257 * zr,zi - approximate zero if convergence (result TRUE)
258 *
259 * Return value indicates convergence of stage 3 iteration
260 *
261 * Uses global (sr,si), nn, pr[], pi[], .. (all args of polyev() !)
262 */
263
264 Rboolean pasd, bool, test;
265 static double svsi, svsr;
266 static int i, j, n;
267 static double oti, otr;
268
269 n = nn - 1;
270
271 /* evaluate p at s. */
272
273 polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);
274
275 test = TRUE;
276 pasd = FALSE;
277
278 /* calculate first t = -p(s)/h(s). */
279
280 calct(&bool);
281
282 /* main loop for one second stage step. */
283
284 for (j=1; j<=l2; j++) {
285
286 otr = tr;
287 oti = ti;
288
289 /* compute next h polynomial and new t. */
290
291 nexth(bool);
292 calct(&bool);
293 *zr = sr + tr;
294 *zi = si + ti;
295
296 /* test for convergence unless stage 3 has */
297 /* failed once or this is the last h polynomial. */
298
299 if (!bool && test && j != l2) {
300 if (hypot(tr - otr, ti - oti) >= hypot(*zr, *zi) * 0.5) {
301 pasd = FALSE;
302 }
303 else if (! pasd) {
304 pasd = TRUE;
305 }
306 else {
307
308 /* the weak convergence test has been */
309 /* passed twice, start the third stage */
310 /* iteration, after saving the current */
311 /* h polynomial and shift. */
312
313 for (i = 0; i < n; i++) {
314 shr[i] = hr[i];
315 shi[i] = hi[i];
316 }
317 svsr = sr;
318 svsi = si;
319 if (vrshft(10, zr, zi)) {
320 return TRUE;
321 }
322
323 /* the iteration failed to converge. */
324 /* turn off testing and restore */
325 /* h, s, pv and t. */
326
327 test = FALSE;
328 for (i=1 ; i<=n ; i++) {
329 hr[i-1] = shr[i-1];
330 hi[i-1] = shi[i-1];
331 }
332 sr = svsr;
333 si = svsi;
334 polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);
335 calct(&bool);
336 }
337 }
338 }
339
340 /* attempt an iteration with final h polynomial */
341 /* from second stage. */
342
343 return(vrshft(10, zr, zi));
344 }
345
346
347 /* carries out the third stage iteration.
348 */
vrshft(int l3,double * zr,double * zi)349 static Rboolean vrshft(int l3, double *zr, double *zi)
350 {
351 /* l3 - limit of steps in stage 3.
352 * zr,zi - on entry contains the initial iterate;
353 * if the iteration converges it contains
354 * the final iterate on exit.
355 * Returns TRUE if iteration converges
356 *
357 * Assign and uses GLOBAL sr, si
358 */
359 Rboolean bool, b;
360 static int i, j;
361 static double r1, r2, mp, ms, tp, relstp;
362 static double omp;
363
364 b = FALSE;
365 sr = *zr;
366 si = *zi;
367
368 /* main loop for stage three */
369
370 for (i = 1; i <= l3; i++) {
371
372 /* evaluate p at s and test for convergence. */
373 polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);
374
375 mp = hypot(pvr, pvi);
376 ms = hypot(sr, si);
377 if (mp <= 20. * errev(nn, qpr, qpi, ms, mp, /*are=*/eta, mre)) {
378 goto L_conv;
379 }
380
381 /* polynomial value is smaller in value than */
382 /* a bound on the error in evaluating p, */
383 /* terminate the iteration. */
384
385 if (i != 1) {
386
387 if (!b && mp >= omp && relstp < .05) {
388
389 /* iteration has stalled. probably a */
390 /* cluster of zeros. do 5 fixed shift */
391 /* steps into the cluster to force */
392 /* one zero to dominate. */
393
394 tp = relstp;
395 b = TRUE;
396 if (relstp < eta)
397 tp = eta;
398 r1 = sqrt(tp);
399 r2 = sr * (r1 + 1.) - si * r1;
400 si = sr * r1 + si * (r1 + 1.);
401 sr = r2;
402 polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);
403 for (j = 1; j <= 5; ++j) {
404 calct(&bool);
405 nexth(bool);
406 }
407 omp = infin;
408 goto L10;
409 }
410 else {
411
412 /* exit if polynomial value */
413 /* increases significantly. */
414
415 if (mp * .1 > omp)
416 return FALSE;
417 }
418 }
419 omp = mp;
420
421 /* calculate next iterate. */
422
423 L10:
424 calct(&bool);
425 nexth(bool);
426 calct(&bool);
427 if (!bool) {
428 relstp = hypot(tr, ti) / hypot(sr, si);
429 sr += tr;
430 si += ti;
431 }
432 }
433 return FALSE;
434
435 L_conv:
436 *zr = sr;
437 *zi = si;
438 return TRUE;
439 }
440
calct(Rboolean * bool)441 static void calct(Rboolean *bool)
442 {
443 /* computes t = -p(s)/h(s).
444 * bool - logical, set true if h(s) is essentially zero. */
445
446 int n = nn - 1;
447 double hvi, hvr;
448
449 /* evaluate h(s). */
450 polyev(n, sr, si, hr, hi,
451 qhr, qhi, &hvr, &hvi);
452
453 *bool = hypot(hvr, hvi) <= are * 10. * hypot(hr[n-1], hi[n-1]);
454 if (!*bool) {
455 cdivid(-pvr, -pvi, hvr, hvi, &tr, &ti);
456 }
457 else {
458 tr = 0.;
459 ti = 0.;
460 }
461 }
462
nexth(Rboolean bool)463 static void nexth(Rboolean bool)
464 {
465 /* calculates the next shifted h polynomial.
466 * bool : if TRUE h(s) is essentially zero
467 */
468 int j, n = nn - 1;
469 double t1, t2;
470
471 if (!bool) {
472 for (j=1; j < n; j++) {
473 t1 = qhr[j - 1];
474 t2 = qhi[j - 1];
475 hr[j] = tr * t1 - ti * t2 + qpr[j];
476 hi[j] = tr * t2 + ti * t1 + qpi[j];
477 }
478 hr[0] = qpr[0];
479 hi[0] = qpi[0];
480 }
481 else {
482 /* if h(s) is zero replace h with qh. */
483
484 for (j=1; j < n; j++) {
485 hr[j] = qhr[j-1];
486 hi[j] = qhi[j-1];
487 }
488 hr[0] = 0.;
489 hi[0] = 0.;
490 }
491 }
492
493 /*--------------------- Independent Complex Polynomial Utilities ----------*/
494
495 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)496 void polyev(int n,
497 double s_r, double s_i,
498 double *p_r, double *p_i,
499 double *q_r, double *q_i,
500 double *v_r, double *v_i)
501 {
502 /* evaluates a polynomial p at s by the horner recurrence
503 * placing the partial sums in q and the computed value in v_.
504 */
505 int i;
506 double t;
507
508 q_r[0] = p_r[0];
509 q_i[0] = p_i[0];
510 *v_r = q_r[0];
511 *v_i = q_i[0];
512 for (i = 1; i < n; i++) {
513 t = *v_r * s_r - *v_i * s_i + p_r[i];
514 q_i[i] = *v_i = *v_r * s_i + *v_i * s_r + p_i[i];
515 q_r[i] = *v_r = t;
516 }
517 }
518
519 static
errev(int n,double * qr,double * qi,double ms,double mp,double a_re,double m_re)520 double errev(int n, double *qr, double *qi,
521 double ms, double mp, double a_re, double m_re)
522 {
523 /* bounds the error in evaluating the polynomial by the horner
524 * recurrence.
525 *
526 * qr,qi - the partial sum vectors
527 * ms - modulus of the point
528 * mp - modulus of polynomial value
529 * a_re,m_re - error bounds on complex addition and multiplication
530 */
531 double e;
532 int i;
533
534 e = hypot(qr[0], qi[0]) * m_re / (a_re + m_re);
535 for (i=0; i < n; i++)
536 e = e*ms + hypot(qr[i], qi[i]);
537
538 return e * (a_re + m_re) - mp * m_re;
539 }
540
541
542 static
cpoly_cauchy(int n,double * pot,double * q)543 double cpoly_cauchy(int n, double *pot, double *q)
544 {
545 /* Computes a lower bound on the moduli of the zeros of a polynomial
546 * pot[1:nn] is the modulus of the coefficients.
547 */
548 double f, x, delf, dx, xm;
549 int i, n1 = n - 1;
550
551 pot[n1] = -pot[n1];
552
553 /* compute upper estimate of bound. */
554
555 x = exp((log(-pot[n1]) - log(pot[0])) / (double) n1);
556
557 /* if newton step at the origin is better, use it. */
558
559 if (pot[n1-1] != 0.) {
560 xm = -pot[n1] / pot[n1-1];
561 if (xm < x)
562 x = xm;
563 }
564
565 /* chop the interval (0,x) unitl f le 0. */
566
567 for(;;) {
568 xm = x * 0.1;
569 f = pot[0];
570 for (i = 1; i < n; i++)
571 f = f * xm + pot[i];
572 if (f <= 0.0) {
573 break;
574 }
575 x = xm;
576 }
577
578 dx = x;
579
580 /* do Newton iteration until x converges to two decimal places. */
581
582 while (fabs(dx / x) > 0.005) {
583 q[0] = pot[0];
584 for(i = 1; i < n; i++)
585 q[i] = q[i-1] * x + pot[i];
586 f = q[n1];
587 delf = q[0];
588 for(i = 1; i < n1; i++)
589 delf = delf * x + q[i];
590 dx = f / delf;
591 x -= dx;
592 }
593 return x;
594 }
595
596 static
cpoly_scale(int n,double * pot,double eps,double BIG,double small,double base)597 double cpoly_scale(int n, double *pot,
598 double eps, double BIG, double small, double base)
599 {
600 /* Returns a scale factor to multiply the coefficients of the polynomial.
601 * The scaling is done to avoid overflow and to avoid
602 * undetected underflow interfering with the convergence criterion.
603 * The factor is a power of the base.
604
605 * pot[1:n] : modulus of coefficients of p
606 * eps,BIG,
607 * small,base - constants describing the floating point arithmetic.
608 */
609
610 int i, ell;
611 double x, high, sc, lo, min_, max_;
612
613 /* find largest and smallest moduli of coefficients. */
614 high = sqrt(BIG);
615 lo = small / eps;
616 max_ = 0.;
617 min_ = BIG;
618 for (i = 0; i < n; i++) {
619 x = pot[i];
620 if (x > max_) max_ = x;
621 if (x != 0. && x < min_)
622 min_ = x;
623 }
624
625 /* scale only if there are very large or very small components. */
626
627 if (min_ < lo || max_ > high) {
628 x = lo / min_;
629 if (x <= 1.)
630 sc = 1. / (sqrt(max_) * sqrt(min_));
631 else {
632 sc = x;
633 if (BIG / sc > max_)
634 sc = 1.0;
635 }
636 ell = (int) (log(sc) / log(base) + 0.5);
637 return R_pow_di(base, ell);
638 }
639 else return 1.0;
640 }
641
642
643 static
cdivid(double ar,double ai,double br,double bi,double * cr,double * ci)644 void cdivid(double ar, double ai, double br, double bi,
645 double *cr, double *ci)
646 {
647 /* complex division c = a/b, i.e., (cr +i*ci) = (ar +i*ai) / (br +i*bi),
648 avoiding overflow. */
649
650 double d, r;
651
652 if (br == 0. && bi == 0.) {
653 /* division by zero, c = infinity. */
654 *cr = *ci = R_PosInf;
655 }
656 else if (fabs(br) >= fabs(bi)) {
657 r = bi / br;
658 d = br + r * bi;
659 *cr = (ar + ai * r) / d;
660 *ci = (ai - ar * r) / d;
661 }
662 else {
663 r = br / bi;
664 d = bi + r * br;
665 *cr = (ar * r + ai) / d;
666 *ci = (ai * r - ar) / d;
667 }
668 }
669
670 /* static double cpoly_cmod(double *r, double *i)
671 * --> replaced by hypot() everywhere
672 */
673