The R Project SVN
R
(root)/branches/djm-parseRd/src/appl/cpoly.c – Rev 46659
Rev 45774 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
/** R : A Computer Language for Statistical Data Analysis* Copyright (C) 1997-1998 Ross Ihaka* Copyright (C) 1999-2001 R Development Core Team** This program is free software; you can redistribute it and/or modify* it under the terms of the GNU General Public License as published by* the Free Software Foundation; either version 2 of the License, or* (at your option) any later version.** This program is distributed in the hope that it will be useful,* but WITHOUT ANY WARRANTY; without even the implied warranty of* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the* GNU General Public License for more details.** You should have received a copy of the GNU General Public License* along with this program; if not, a copy is available at* http://www.r-project.org/Licenses/** cpoly finds the zeros of a complex polynomial.** On Entry** opr, opi - double precision vectors of real and* imaginary parts of the coefficients in* order of decreasing powers.** degree - int degree of polynomial.*** On Return** zeror, zeroi - output double precision vectors of* real and imaginary parts of the zeros.** fail - output int parameter, true only if* leading coefficient is zero or if cpoly* has found fewer than degree zeros.** The program has been written to reduce the chance of overflow* occurring. If it does occur, there is still a possibility that* the zerofinder will work provided the overflowed quantity is* replaced by a large number.** This is a C translation of the following.** TOMS Algorithm 419* Jenkins and Traub.* Comm. ACM 15 (1972) 97-99.** Ross Ihaka* February 1997*/#ifdef HAVE_CONFIG_H#include <config.h>#endif#include <R_ext/Arith.h> /* for declaration of hypot */#include <R_ext/Memory.h> /* for declaration of R_alloc */#include <float.h> /* for FLT_RADIX */#include <Rmath.h> /* for R_pow_di */#include <R_ext/Applic.h>#ifdef HAVE_VISIBILITY_ATTRIBUTE# define attribute_hidden __attribute__ ((visibility ("hidden")))#else# define attribute_hidden#endif#ifndef HAVE_HYPOT# define hypot pythag#endifstatic void calct(Rboolean *);static Rboolean fxshft(int, double *, double *);static Rboolean vrshft(int, double *, double *);static void nexth(Rboolean);static void noshft(int);/* Consider exporting these (via Applic.h): */static void polyev(int, double, double,double *, double *, double *, double *, double *, double *);static double errev(int, double *, double *, double, double, double, double);static double cpoly_cauchy(int, double *, double *);static double cpoly_scale(int, double *, double, double, double, double);static void cdivid(double, double, double, double, double *, double *);/* Global Variables (too many!) */static int nn;#if 0#define NMAX 50static double pr[NMAX];static double pi[NMAX];static double hr[NMAX];static double hi[NMAX];static double qpr[NMAX];static double qpi[NMAX];static double qhr[NMAX];static double qhi[NMAX];static double shr[NMAX];static double shi[NMAX];#elsestatic double *pr, *pi, *hr, *hi, *qpr, *qpi, *qhr, *qhi, *shr, *shi;#endifstatic double sr, si;static double tr, ti;static double pvr, pvi;static const double eta = DBL_EPSILON;static const double are = /* eta = */DBL_EPSILON;static const double mre = 2. * M_SQRT2 * /* eta, i.e. */DBL_EPSILON;static const double infin = DBL_MAX;attribute_hiddenvoid R_cpolyroot(double *opr, double *opi, int *degree,double *zeror, double *zeroi, Rboolean *fail){static const double smalno = DBL_MIN;static const double base = (double)FLT_RADIX;static int d_n, i, i1, i2;static double zi, zr, xx, yy;static double bnd, xxx;Rboolean conv;int d1;double *tmp;static const double cosr =/* cos 94 */ -0.06975647374412529990;static const double sinr =/* sin 94 */ 0.99756405025982424767;xx = M_SQRT1_2;/* 1/sqrt(2) = 0.707.... */yy = -xx;*fail = FALSE;nn = *degree;d1 = nn - 1;/* algorithm fails if the leading coefficient is zero. */if (opr[0] == 0. && opi[0] == 0.) {*fail = TRUE;return;}/* remove the zeros at the origin if any. */while (opr[nn] == 0. && opi[nn] == 0.) {d_n = d1-nn+1;zeror[d_n] = 0.;zeroi[d_n] = 0.;nn--;}nn++;/*-- Now, global var. nn := #{coefficients} = (relevant degree)+1 */if (nn == 1) return;/* Use a single allocation as these as small */tmp = (double *) R_alloc(10*nn, sizeof(double));pr = tmp; pi = tmp + nn; hr = tmp + 2*nn; hi = tmp + 3*nn;qpr = tmp + 4*nn; qpi = tmp + 5*nn; qhr = tmp + 6*nn; qhi = tmp + 7*nn;shr = tmp + 8*nn; shi = tmp + 9*nn;/* make a copy of the coefficients and shr[] = | p[] | */for (i = 0; i < nn; i++) {pr[i] = opr[i];pi[i] = opi[i];shr[i] = hypot(pr[i], pi[i]);}/* scale the polynomial with factor 'bnd'. */bnd = cpoly_scale(nn, shr, eta, infin, smalno, base);if (bnd != 1.) {for (i=0; i < nn; i++) {pr[i] *= bnd;pi[i] *= bnd;}}/* start the algorithm for one zero */while (nn > 2) {/* calculate bnd, a lower bound on the modulus of the zeros. */for (i=0 ; i < nn ; i++)shr[i] = hypot(pr[i], pi[i]);bnd = cpoly_cauchy(nn, shr, shi);/* outer loop to control 2 major passes *//* with different sequences of shifts */for (i1 = 1; i1 <= 2; i1++) {/* first stage calculation, no shift */noshft(5);/* inner loop to select a shift */for (i2 = 1; i2 <= 9; i2++) {/* shift is chosen with modulus bnd *//* and amplitude rotated by 94 degrees *//* from the previous shift */xxx= cosr * xx - sinr * yy;yy = sinr * xx + cosr * yy;xx = xxx;sr = bnd * xx;si = bnd * yy;/* second stage calculation, fixed shift */conv = fxshft(i2 * 10, &zr, &zi);if (conv)goto L10;}}/* the zerofinder has failed on two major passes *//* return empty handed */*fail = TRUE;return;/* the second stage jumps directly to the third stage iteration.* if successful, the zero is stored and the polynomial deflated.*/L10:d_n = d1+2 - nn;zeror[d_n] = zr;zeroi[d_n] = zi;--nn;for (i=0; i < nn ; i++) {pr[i] = qpr[i];pi[i] = qpi[i];}}/*while*//* calculate the final zero and return */cdivid(-pr[1], -pi[1], pr[0], pi[0], &zeror[d1], &zeroi[d1]);return;}/* Computes the derivative polynomial as the initial* polynomial and computes l1 no-shift h polynomials. */static void noshft(int l1){int i, j, jj, n = nn - 1, nm1 = n - 1;double t1, t2, xni;for (i=0; i < n; i++) {xni = (double)(nn - i - 1);hr[i] = xni * pr[i] / n;hi[i] = xni * pi[i] / n;}for (jj = 1; jj <= l1; jj++) {if (hypot(hr[n-1], hi[n-1]) <=eta * 10.0 * hypot(pr[n-1], pi[n-1])) {/* If the constant term is essentially zero, *//* shift h coefficients. */for (i = 1; i <= nm1; i++) {j = nn - i;hr[j-1] = hr[j-2];hi[j-1] = hi[j-2];}hr[0] = 0.;hi[0] = 0.;}else {cdivid(-pr[nn-1], -pi[nn-1], hr[n-1], hi[n-1], &tr, &ti);for (i = 1; i <= nm1; i++) {j = nn - i;t1 = hr[j-2];t2 = hi[j-2];hr[j-1] = tr * t1 - ti * t2 + pr[j-1];hi[j-1] = tr * t2 + ti * t1 + pi[j-1];}hr[0] = pr[0];hi[0] = pi[0];}}}/* Computes l2 fixed-shift h polynomials and tests for convergence.* initiates a variable-shift iteration and returns with the* approximate zero if successful.*/static Rboolean fxshft(int l2, double *zr, double *zi){/* l2 - limit of fixed shift steps* zr,zi - approximate zero if convergence (result TRUE)** Return value indicates convergence of stage 3 iteration** Uses global (sr,si), nn, pr[], pi[], .. (all args of polyev() !)*/Rboolean pasd, bool, test;static double svsi, svsr;static int i, j, n;static double oti, otr;n = nn - 1;/* evaluate p at s. */polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);test = TRUE;pasd = FALSE;/* calculate first t = -p(s)/h(s). */calct(&bool);/* main loop for one second stage step. */for (j=1; j<=l2; j++) {otr = tr;oti = ti;/* compute next h polynomial and new t. */nexth(bool);calct(&bool);*zr = sr + tr;*zi = si + ti;/* test for convergence unless stage 3 has *//* failed once or this is the last h polynomial. */if (!bool && test && j != l2) {if (hypot(tr - otr, ti - oti) >= hypot(*zr, *zi) * 0.5) {pasd = FALSE;}else if (! pasd) {pasd = TRUE;}else {/* the weak convergence test has been *//* passed twice, start the third stage *//* iteration, after saving the current *//* h polynomial and shift. */for (i = 0; i < n; i++) {shr[i] = hr[i];shi[i] = hi[i];}svsr = sr;svsi = si;if (vrshft(10, zr, zi)) {return TRUE;}/* the iteration failed to converge. *//* turn off testing and restore *//* h, s, pv and t. */test = FALSE;for (i=1 ; i<=n ; i++) {hr[i-1] = shr[i-1];hi[i-1] = shi[i-1];}sr = svsr;si = svsi;polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);calct(&bool);}}}/* attempt an iteration with final h polynomial *//* from second stage. */return(vrshft(10, zr, zi));}/* carries out the third stage iteration.*/static Rboolean vrshft(int l3, double *zr, double *zi){/* l3 - limit of steps in stage 3.* zr,zi - on entry contains the initial iterate;* if the iteration converges it contains* the final iterate on exit.* Returns TRUE if iteration converges** Assign and uses GLOBAL sr, si*/Rboolean bool, b;static int i, j;static double r1, r2, mp, ms, tp, relstp;static double omp;b = FALSE;sr = *zr;si = *zi;/* main loop for stage three */for (i = 1; i <= l3; i++) {/* evaluate p at s and test for convergence. */polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);mp = hypot(pvr, pvi);ms = hypot(sr, si);if (mp <= 20. * errev(nn, qpr, qpi, ms, mp, /*are=*/eta, mre)) {goto L_conv;}/* polynomial value is smaller in value than *//* a bound on the error in evaluating p, *//* terminate the iteration. */if (i != 1) {if (!b && mp >= omp && relstp < .05) {/* iteration has stalled. probably a *//* cluster of zeros. do 5 fixed shift *//* steps into the cluster to force *//* one zero to dominate. */tp = relstp;b = TRUE;if (relstp < eta)tp = eta;r1 = sqrt(tp);r2 = sr * (r1 + 1.) - si * r1;si = sr * r1 + si * (r1 + 1.);sr = r2;polyev(nn, sr, si, pr, pi, qpr, qpi, &pvr, &pvi);for (j = 1; j <= 5; ++j) {calct(&bool);nexth(bool);}omp = infin;goto L10;}else {/* exit if polynomial value *//* increases significantly. */if (mp * .1 > omp)return FALSE;}}omp = mp;/* calculate next iterate. */L10:calct(&bool);nexth(bool);calct(&bool);if (!bool) {relstp = hypot(tr, ti) / hypot(sr, si);sr += tr;si += ti;}}return FALSE;L_conv:*zr = sr;*zi = si;return TRUE;}static void calct(Rboolean *bool){/* computes t = -p(s)/h(s).* bool - logical, set true if h(s) is essentially zero. */int n = nn - 1;double hvi, hvr;/* evaluate h(s). */polyev(n, sr, si, hr, hi,qhr, qhi, &hvr, &hvi);*bool = hypot(hvr, hvi) <= are * 10. * hypot(hr[n-1], hi[n-1]);if (!*bool) {cdivid(-pvr, -pvi, hvr, hvi, &tr, &ti);}else {tr = 0.;ti = 0.;}}static void nexth(Rboolean bool){/* calculates the next shifted h polynomial.* bool : if TRUE h(s) is essentially zero*/int j, n = nn - 1;double t1, t2;if (!bool) {for (j=1; j < n; j++) {t1 = qhr[j - 1];t2 = qhi[j - 1];hr[j] = tr * t1 - ti * t2 + qpr[j];hi[j] = tr * t2 + ti * t1 + qpi[j];}hr[0] = qpr[0];hi[0] = qpi[0];}else {/* if h(s) is zero replace h with qh. */for (j=1; j < n; j++) {hr[j] = qhr[j-1];hi[j] = qhi[j-1];}hr[0] = 0.;hi[0] = 0.;}}/*--------------------- Independent Complex Polynomial Utilities ----------*/staticvoid 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){/* evaluates a polynomial p at s by the horner recurrence* placing the partial sums in q and the computed value in v_.*/int i;double t;q_r[0] = p_r[0];q_i[0] = p_i[0];*v_r = q_r[0];*v_i = q_i[0];for (i = 1; i < n; i++) {t = *v_r * s_r - *v_i * s_i + p_r[i];q_i[i] = *v_i = *v_r * s_i + *v_i * s_r + p_i[i];q_r[i] = *v_r = t;}}staticdouble errev(int n, double *qr, double *qi,double ms, double mp, double a_re, double m_re){/* bounds the error in evaluating the polynomial by the horner* recurrence.** qr,qi - the partial sum vectors* ms - modulus of the point* mp - modulus of polynomial value* a_re,m_re - error bounds on complex addition and multiplication*/double e;int i;e = hypot(qr[0], qi[0]) * m_re / (a_re + m_re);for (i=0; i < n; i++)e = e*ms + hypot(qr[i], qi[i]);return e * (a_re + m_re) - mp * m_re;}staticdouble cpoly_cauchy(int n, double *pot, double *q){/* Computes a lower bound on the moduli of the zeros of a polynomial* pot[1:nn] is the modulus of the coefficients.*/double f, x, delf, dx, xm;int i, n1 = n - 1;pot[n1] = -pot[n1];/* compute upper estimate of bound. */x = exp((log(-pot[n1]) - log(pot[0])) / (double) n1);/* if newton step at the origin is better, use it. */if (pot[n1-1] != 0.) {xm = -pot[n1] / pot[n1-1];if (xm < x)x = xm;}/* chop the interval (0,x) unitl f le 0. */for(;;) {xm = x * 0.1;f = pot[0];for (i = 1; i < n; i++)f = f * xm + pot[i];if (f <= 0.0) {break;}x = xm;}dx = x;/* do Newton iteration until x converges to two decimal places. */while (fabs(dx / x) > 0.005) {q[0] = pot[0];for(i = 1; i < n; i++)q[i] = q[i-1] * x + pot[i];f = q[n1];delf = q[0];for(i = 1; i < n1; i++)delf = delf * x + q[i];dx = f / delf;x -= dx;}return x;}staticdouble cpoly_scale(int n, double *pot,double eps, double BIG, double small, double base){/* Returns a scale factor to multiply the coefficients of the polynomial.* The scaling is done to avoid overflow and to avoid* undetected underflow interfering with the convergence criterion.* The factor is a power of the base.* pot[1:n] : modulus of coefficients of p* eps,BIG,* small,base - constants describing the floating point arithmetic.*/int i, ell;double x, high, sc, lo, min_, max_;/* find largest and smallest moduli of coefficients. */high = sqrt(BIG);lo = small / eps;max_ = 0.;min_ = BIG;for (i = 0; i < n; i++) {x = pot[i];if (x > max_) max_ = x;if (x != 0. && x < min_)min_ = x;}/* scale only if there are very large or very small components. */if (min_ < lo || max_ > high) {x = lo / min_;if (x <= 1.)sc = 1. / (sqrt(max_) * sqrt(min_));else {sc = x;if (BIG / sc > max_)sc = 1.0;}ell = (int) (log(sc) / log(base) + 0.5);return R_pow_di(base, ell);}else return 1.0;}staticvoid cdivid(double ar, double ai, double br, double bi,double *cr, double *ci){/* complex division c = a/b, i.e., (cr +i*ci) = (ar +i*ai) / (br +i*bi),avoiding overflow. */double d, r;if (br == 0. && bi == 0.) {/* division by zero, c = infinity. */*cr = *ci = R_PosInf;}else if (fabs(br) >= fabs(bi)) {r = bi / br;d = br + r * bi;*cr = (ar + ai * r) / d;*ci = (ai - ar * r) / d;}else {r = br / bi;d = bi + r * br;*cr = (ar * r + ai) / d;*ci = (ai * r - ar) / d;}}/* static double cpoly_cmod(double *r, double *i)* --> replaced by hypot() everywhere*/