BRL-CAD
roots.c
Go to the documentation of this file.
1 /* R O O T S . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1985-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 
21 
22 #include "common.h"
23 
24 #include <math.h>
25 #include "bio.h"
26 
27 
28 #include "vmath.h"
29 #include "bn.h"
30 #include "raytrace.h"
31 
32 
33 static const bn_poly_t bn_Zero_poly = { BN_POLY_MAGIC, 0, {0.0} };
34 
35 
36 /**
37  * Evaluates p(Z), p'(Z), and p''(Z) for any Z (real or complex).
38  * Given an equation of the form:
39  *
40  * p(Z) = a0*Z^n + a1*Z^(n-1) +... an != 0,
41  *
42  * the function value and derivatives needed for the iterations can be
43  * computed by synthetic division using the formulas
44  *
45  * p(Z) = bn, p'(Z) = c(n-1), p''(Z) = d(n-2),
46  *
47  * where
48  *
49  * b0 = a0, bi = b(i-1)*Z + ai, i = 1, 2, ...n
50  * c0 = b0, ci = c(i-1)*Z + bi, i = 1, 2, ...n-1
51  * d0 = c0, di = d(i-1)*Z + ci, i = 1, 2, ...n-2
52  */
53 void
54 rt_poly_eval_w_2derivatives(register bn_complex_t *cZ, register bn_poly_t *eqn, register bn_complex_t *b, register bn_complex_t *c, register bn_complex_t *d)
55 /* input */
56 /* input */
57 /* outputs */
58 {
59  register size_t n;
60  register int m;
61 
62  bn_cx_cons(b, eqn->cf[0], 0.0);
63  *c = *b;
64  *d = *c;
65 
66  for (n=1; (m = eqn->dgr - n) >= 0; ++n) {
67  bn_cx_mul(b, cZ);
68  b->re += eqn->cf[n];
69  if (m > 0) {
70  bn_cx_mul(c, cZ);
71  bn_cx_add(c, b);
72  }
73  if (m > 1) {
74  bn_cx_mul(d, cZ);
75  bn_cx_add(d, c);
76  }
77  }
78 }
79 
80 
81 /**
82  * Calculates one root of a polynomial (p(Z)) using Laguerre's
83  * method. This is an iterative technique which has very good global
84  * convergence properties. The formulas for this method are
85  *
86  * n * p(Z)
87  * newZ = Z - ----------------------- ,
88  * p'(Z) +- sqrt(H(Z))
89  *
90  * where
91  * H(Z) = (n-1) [ (n-1)(p'(Z))^2 - n p(Z)p''(Z) ],
92  *
93  * where n is the degree of the polynomial. The sign in the
94  * denominator is chosen so that |newZ - Z| is as small as possible.
95  */
96 int
97 rt_poly_findroot(register bn_poly_t *eqn, /* polynomial */
98  register bn_complex_t *nxZ, /* initial guess for root */
99  const char *str)
100 {
101  bn_complex_t p0 = {0.0, 0.0}, p1 = {0.0, 0.0}, p2 = {0.0, 0.0}; /* evaluated polynomial+derivatives */
102  bn_complex_t p1_H; /* p1 - H, temporary */
103  bn_complex_t cZ, cH; /* 'Z' and H(Z) in comment */
104  bn_complex_t T; /* temporary for making H */
105  fastf_t diff=0.0; /* test values for convergence */
106  fastf_t b=0.0; /* floating temps */
107  int n;
108  int i;
109 
110  for (i=0; i < 100; i++) {
111  cZ = *nxZ;
112  rt_poly_eval_w_2derivatives(&cZ, eqn, &p0, &p1, &p2);
113 
114  /* Compute H for Laguerre's method. */
115  n = eqn->dgr-1;
116  bn_cx_mul2(&cH, &p1, &p1);
117  bn_cx_scal(&cH, (double)(n*n));
118  bn_cx_mul2(&T, &p2, &p0);
119  bn_cx_scal(&T, (double)(eqn->dgr*n));
120  bn_cx_sub(&cH, &T);
121 
122  /* Calculate the next iteration for Laguerre's method. Test
123  * to see whether addition or subtraction gives the larger
124  * denominator for the next 'Z', and use the appropriate value
125  * in the formula.
126  */
127  bn_cx_sqrt(&cH, &cH);
128  p1_H = p1;
129  bn_cx_sub(&p1_H, &cH);
130  bn_cx_add(&p1, &cH); /* p1 <== p1+H */
131  bn_cx_scal(&p0, (double)(eqn->dgr));
132  if (bn_cx_amplsq(&p1_H) > bn_cx_amplsq(&p1)) {
133  bn_cx_div(&p0, &p1_H);
134  bn_cx_sub(nxZ, &p0);
135  } else {
136  bn_cx_div(&p0, &p1);
137  bn_cx_sub(nxZ, &p0);
138  }
139 
140  /* Use proportional convergence test to allow very small roots
141  * and avoid wasting time on large roots. The original
142  * version used bn_cx_ampl(), which requires a square root.
143  * Using bn_cx_amplsq() saves lots of cycles, but changes loop
144  * termination conditions somewhat.
145  *
146  * diff is |p0|**2. nxZ = Z - p0.
147  */
148  b = bn_cx_amplsq(nxZ);
149  diff = bn_cx_amplsq(&p0);
150 
151  if (b < diff)
152  continue;
153 
154  if (ZERO(diff))
155  return i; /* OK -- can't do better */
156 
157  /* FIXME: figure out why SMALL_FASTF is too sensitive, why
158  * anything smaller than 1.0e-5 is too sensitive and causes
159  * eto off-by-many differences.
160  */
161  if (diff > (b - diff) * 1.0e-5)
162  continue;
163 
164  return i; /* OK */
165  }
166 
167  /* If the thing hasn't converged yet, it probably won't. */
168  bu_log("rt_poly_findroot: solving for %s didn't converge in %d iterations, b=%g, diff=%g\n",
169  str, i, b, diff);
170  bu_log("nxZ=%gR+%gI, p0=%gR+%gI\n", nxZ->re, nxZ->im, p0.re, p0.im);
171  return -1; /* ERROR */
172 }
173 
174 
175 /**
176  * Evaluates p(Z) for any Z (real or complex). In this case, test all
177  * "nroots" entries of roots[] to ensure that they are roots (zeros)
178  * of this polynomial.
179  *
180  * Returns -
181  * 0 all roots are true zeros
182  * 1 at least one "root[]" entry is not a true zero
183  *
184  * Given an equation of the form
185  *
186  * p(Z) = a0*Z^n + a1*Z^(n-1) +... an != 0,
187  *
188  * the function value can be computed using the formula
189  *
190  * p(Z) = bn, where
191  *
192  * b0 = a0, bi = b(i-1)*Z + ai, i = 1, 2, ...n
193  */
194 int
195 rt_poly_checkroots(register bn_poly_t *eqn, bn_complex_t *roots, register int nroots)
196 {
197  register fastf_t er, ei; /* "epoly" */
198  register fastf_t zr, zi; /* Z value to evaluate at */
199  register size_t n;
200  int m;
201 
202  for (m=0; m < nroots; ++m) {
203  /* Select value of Z to evaluate at */
204  zr = bn_cx_real(&roots[m]);
205  zi = bn_cx_imag(&roots[m]);
206 
207  /* Initialize */
208  er = eqn->cf[0];
209  /* ei = 0.0; */
210 
211  /* n=1 step. Broken out because ei = 0.0 */
212  ei = er*zi; /* must come before er= */
213  er = er*zr + eqn->cf[1];
214 
215  /* Remaining steps */
216  for (n=2; n <= eqn->dgr; ++n) {
217  register fastf_t tr, ti; /* temps */
218  tr = er*zr - ei*zi + eqn->cf[n];
219  ti = er*zi + ei*zr;
220  er = tr;
221  ei = ti;
222  }
223  if (fabs(er) > 1.0e-5 || fabs(ei) > 1.0e-5)
224  return 1; /* FAIL */
225  }
226  /* Evaluating polynomial for all Z values gives zero */
227  return 0; /* OK */
228 }
229 
230 
231 /**
232  * Deflates a polynomial by a given root.
233  */
234 void
235 rt_poly_deflate(register bn_poly_t *oldP, register bn_complex_t *root)
236 {
237  bn_poly_t divisor = bn_Zero_poly;
238  bn_poly_t rem = bn_Zero_poly;
239 
240  /* Make a polynomial out of the given root: Linear for a real
241  * root, Quadratic for a complex root (since they come in con-
242  * jugate pairs).
243  */
244  if (ZERO(root->im)) {
245  /* root is real */
246  divisor.dgr = 1;
247  divisor.cf[0] = 1;
248  divisor.cf[1] = - root->re;
249  } else {
250  /* root is complex */
251  divisor.dgr = 2;
252  divisor.cf[0] = 1;
253  divisor.cf[1] = -2 * root->re;
254  divisor.cf[2] = bn_cx_amplsq(root);
255  }
256 
257  /* Use synthetic division to find the quotient (new polynomial)
258  * and the remainder (should be zero if the root was really a
259  * root).
260  */
261  bn_poly_synthetic_division(oldP, &rem, oldP, &divisor);
262 }
263 
264 
265 int
266 rt_poly_roots(register bn_poly_t *eqn, /* equation to be solved */
267  register bn_complex_t roots[], /* space to put roots found */
268  const char *name) /* name of the primitive being checked */
269 {
270  register size_t n; /* number of roots found */
271  fastf_t factor; /* scaling factor for copy */
272 
273  /* Remove leading coefficients which are too close to zero,
274  * to prevent the polynomial factoring from blowing up, below.
275  */
276  while (ZERO(eqn->cf[0])) {
277  for (n=0; n <= eqn->dgr; n++) {
278  eqn->cf[n] = eqn->cf[n+1];
279  }
280  if (--eqn->dgr <= 0)
281  return 0;
282  }
283 
284  /* Factor the polynomial so the first coefficient is one
285  * for ease of handling.
286  */
287  factor = 1.0 / eqn->cf[0];
288  (void) bn_poly_scale(eqn, factor);
289  n = 0; /* Number of roots found */
290 
291  /* A trailing coefficient of zero indicates that zero
292  * is a root of the equation.
293  */
294  while (ZERO(eqn->cf[eqn->dgr])) {
295  roots[n].re = roots[n].im = 0.0;
296  --eqn->dgr;
297  ++n;
298  }
299 
300  while (eqn->dgr > 2) {
301  if (eqn->dgr == 4) {
302  if (bn_poly_quartic_roots(&roots[n], eqn)) {
303  if (rt_poly_checkroots(eqn, &roots[n], 4) == 0) {
304  return n+4;
305  }
306  }
307  } else if (eqn->dgr == 3) {
308  if (bn_poly_cubic_roots(&roots[n], eqn)) {
309  if (rt_poly_checkroots(eqn, &roots[n], 3) == 0) {
310  return n+3;
311  }
312  }
313  }
314 
315  /*
316  * Set initial guess for root to almost zero.
317  * This method requires a small nudge off the real axis.
318  */
319  bn_cx_cons(&roots[n], 0.0, SMALL);
320  if ((rt_poly_findroot(eqn, &roots[n], name)) < 0)
321  return n; /* return those we found, anyways */
322 
323  if (fabs(roots[n].im) > 1.0e-5* fabs(roots[n].re)) {
324  /* If root is complex, its complex conjugate is
325  * also a root since complex roots come in con-
326  * jugate pairs when all coefficients are real.
327  */
328  ++n;
329  roots[n] = roots[n-1];
330  bn_cx_conj(&roots[n]);
331  } else {
332  /* Change 'practically real' to real */
333  roots[n].im = 0.0;
334  }
335 
336  rt_poly_deflate(eqn, &roots[n]);
337  ++n;
338  }
339 
340  /* For polynomials of lower degree, iterative techniques
341  * are an inefficient way to find the roots.
342  */
343  if (eqn->dgr == 1) {
344  roots[n].re = -(eqn->cf[1]);
345  roots[n].im = 0.0;
346  ++n;
347  } else if (eqn->dgr == 2) {
348  bn_poly_quadratic_roots(&roots[n], eqn);
349  n += 2;
350  }
351  return n;
352 }
353 
354 
355 /*
356  * Local Variables:
357  * mode: C
358  * tab-width: 8
359  * indent-tabs-mode: t
360  * c-file-style: "stroustrup"
361  * End:
362  * ex: shiftwidth=4 tabstop=8
363  */
fastf_t cf[BN_MAX_POLY_DEGREE+1]
Definition: poly.h:50
#define bn_cx_cons(cp, r, i)
Definition: complex.h:55
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define bn_cx_sub(ap, bp)
Definition: complex.h:58
#define SMALL
Definition: defines.h:351
int bn_poly_quartic_roots(struct bn_complex roots[], const struct bn_poly *eqn)
Uses the quartic formula to find the roots (in `complex' form) of any quartic equation with real coef...
Definition: clone.c:90
int bn_poly_quadratic_roots(struct bn_complex roots[], const struct bn_poly *quadrat)
Uses the quadratic formula to find the roots (in `complex' form) of any quadratic equation with real ...
Header file for the BRL-CAD common definitions.
Definition: poly.h:47
#define BN_POLY_MAGIC
Definition: magic.h:70
double im
imaginary part
Definition: complex.h:41
void bn_poly_synthetic_division(struct bn_poly *quo, struct bn_poly *rem, const struct bn_poly *dvdend, const struct bn_poly *dvsor)
Divides any polynomial into any other polynomial using synthetic division. Both polynomials must have...
int rt_poly_roots(register bn_poly_t *eqn, register bn_complex_t roots[], const char *name)
Definition: roots.c:266
double re
real part
Definition: complex.h:40
void rt_poly_eval_w_2derivatives(register bn_complex_t *cZ, register bn_poly_t *eqn, register bn_complex_t *b, register bn_complex_t *c, register bn_complex_t *d)
Definition: roots.c:54
size_t dgr
Definition: poly.h:49
int rt_poly_checkroots(register bn_poly_t *eqn, bn_complex_t *roots, register int nroots)
Definition: roots.c:195
struct bn_poly * bn_poly_scale(struct bn_poly *eqn, double factor)
scale a polynomial
#define bn_cx_real(cp)
Definition: complex.h:48
#define bn_cx_add(ap, bp)
Definition: complex.h:51
#define bn_cx_mul2(ap, bp, cp)
Definition: complex.h:66
#define bn_cx_amplsq(cp)
Definition: complex.h:53
int bn_poly_cubic_roots(struct bn_complex roots[], const struct bn_poly *eqn)
Uses the cubic formula to find the roots (in `complex' form) of any cubic equation with real coeffici...
#define ZERO(val)
Definition: units.c:38
void rt_poly_deflate(register bn_poly_t *oldP, register bn_complex_t *root)
Definition: roots.c:235
void bn_cx_div(bn_complex_t *ap, const bn_complex_t *bp)
Divide one complex by another.
struct bn_poly bn_Zero_poly
Definition: bn_poly_add.c:38
int rt_poly_findroot(register bn_poly_t *eqn, register bn_complex_t *nxZ, const char *str)
Definition: roots.c:97
#define bn_cx_scal(cp, s)
Definition: complex.h:57
#define bn_cx_mul(ap, bp)
Definition: complex.h:60
Complex numbers.
Definition: complex.h:39
bn_poly_t rem[1]
double fastf_t
Definition: defines.h:300
void bn_cx_sqrt(bn_complex_t *op, const bn_complex_t *ip)
Compute square root of complex number.
Definition: complex.c:66
#define bn_cx_imag(cp)
Definition: complex.h:49
#define bn_cx_conj(cp)
Definition: complex.h:54