BRL-CAD
poly.c
Go to the documentation of this file.
1 /* P O L Y . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2004-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 /** @addtogroup poly */
22 /** @{ */
23 /** @file libbn/poly.c
24  *
25  * Library for dealing with polynomials.
26  *
27  */
28 
29 #include "common.h"
30 
31 #include <stdlib.h> /* for abs */
32 #include <stdio.h>
33 #include <math.h>
34 
35 #include "bu/log.h"
36 #include "bu/parallel.h"
37 #include "vmath.h"
38 #include "bn/poly.h"
39 
40 
41 #define CUBEROOT(a) (((a) > 0.0) ? pow(a, THIRD) : -pow(-(a), THIRD))
42 
43 
44 static const fastf_t THIRD = 1.0 / 3.0;
45 static const fastf_t TWENTYSEVENTH = 1.0 / 27.0;
46 static const struct bn_poly bn_Zero_poly = { BN_POLY_MAGIC, 0, {0.0} };
47 
48 
49 struct bn_poly *
50 bn_poly_mul(register struct bn_poly *product, register const struct bn_poly *m1, register const struct bn_poly *m2)
51 {
52  if (m1->dgr == 1 && m2->dgr == 1) {
53  product->dgr = 2;
54  product->cf[0] = m1->cf[0] * m2->cf[0];
55  product->cf[1] = m1->cf[0] * m2->cf[1] +
56  m1->cf[1] * m2->cf[0];
57  product->cf[2] = m1->cf[1] * m2->cf[1];
58  return product;
59  }
60  if (m1->dgr == 2 && m2->dgr == 2) {
61  product->dgr = 4;
62  product->cf[0] = m1->cf[0] * m2->cf[0];
63  product->cf[1] = m1->cf[0] * m2->cf[1] +
64  m1->cf[1] * m2->cf[0];
65  product->cf[2] = m1->cf[0] * m2->cf[2] +
66  m1->cf[1] * m2->cf[1] +
67  m1->cf[2] * m2->cf[0];
68  product->cf[3] = m1->cf[1] * m2->cf[2] +
69  m1->cf[2] * m2->cf[1];
70  product->cf[4] = m1->cf[2] * m2->cf[2];
71  return product;
72  }
73 
74  /* Not one of the common (or easy) cases. */
75  {
76  register size_t ct1, ct2;
77 
78  *product = bn_Zero_poly;
79 
80  /* If the degree of the product will be larger than the
81  * maximum size allowed in "polyno.h", then return a null
82  * pointer to indicate failure.
83  */
84  if ((product->dgr = m1->dgr + m2->dgr) > BN_MAX_POLY_DEGREE)
85  return BN_POLY_NULL;
86 
87  for (ct1=0; ct1 <= m1->dgr; ++ct1) {
88  for (ct2=0; ct2 <= m2->dgr; ++ct2) {
89  product->cf[ct1+ct2] +=
90  m1->cf[ct1] * m2->cf[ct2];
91  }
92  }
93  }
94  return product;
95 }
96 
97 
98 struct bn_poly *
99 bn_poly_scale(register struct bn_poly *eqn, double factor)
100 {
101  register size_t cnt;
102 
103  for (cnt=0; cnt <= eqn->dgr; ++cnt) {
104  eqn->cf[cnt] *= factor;
105  }
106  return eqn;
107 }
108 
109 
110 struct bn_poly *
111 bn_poly_add(register struct bn_poly *sum, register const struct bn_poly *poly1, register const struct bn_poly *poly2)
112 {
113  struct bn_poly tmp;
114  register size_t i, offset;
115 
116  offset = labs((long)poly1->dgr - (long)poly2->dgr);
117 
118  tmp = bn_Zero_poly;
119 
120  if (poly1->dgr >= poly2->dgr) {
121  *sum = *poly1;
122  for (i=0; i <= poly2->dgr; ++i) {
123  tmp.cf[i+offset] = poly2->cf[i];
124  }
125  } else {
126  *sum = *poly2;
127  for (i=0; i <= poly1->dgr; ++i) {
128  tmp.cf[i+offset] = poly1->cf[i];
129  }
130  }
131 
132  for (i=0; i <= sum->dgr; ++i) {
133  sum->cf[i] += tmp.cf[i];
134  }
135  return sum;
136 }
137 
138 
139 struct bn_poly *
140 bn_poly_sub(register struct bn_poly *diff, register const struct bn_poly *poly1, register const struct bn_poly *poly2)
141 {
142  struct bn_poly tmp;
143  register size_t i, offset;
144 
145  offset = labs((long)poly1->dgr - (long)poly2->dgr);
146 
147  *diff = bn_Zero_poly;
148  tmp = bn_Zero_poly;
149 
150  if (poly1->dgr >= poly2->dgr) {
151  *diff = *poly1;
152  for (i=0; i <= poly2->dgr; ++i) {
153  tmp.cf[i+offset] = poly2->cf[i];
154  }
155  } else {
156  diff->dgr = poly2->dgr;
157  for (i=0; i <= poly1->dgr; ++i) {
158  diff->cf[i+offset] = poly1->cf[i];
159  }
160  tmp = *poly2;
161  }
162 
163  for (i=0; i <= diff->dgr; ++i) {
164  diff->cf[i] -= tmp.cf[i];
165  }
166  return diff;
167 }
168 
169 
170 void
171 bn_poly_synthetic_division(register struct bn_poly *quo, register struct bn_poly *rem, register const struct bn_poly *dvdend, register const struct bn_poly *dvsor)
172 {
173  register size_t divisor;
174  register size_t n;
175 
176  *quo = *dvdend;
177  *rem = bn_Zero_poly;
178 
179  if (dvsor->dgr > dvdend->dgr) {
180  quo->dgr = -1;
181  } else {
182  quo->dgr = dvdend->dgr - dvsor->dgr;
183  }
184  if ((rem->dgr = dvsor->dgr - 1) > dvdend->dgr)
185  rem->dgr = dvdend->dgr;
186 
187  for (n=0; n <= quo->dgr; ++n) {
188  quo->cf[n] /= dvsor->cf[0];
189  for (divisor=1; divisor <= dvsor->dgr; ++divisor) {
190  quo->cf[n+divisor] -= quo->cf[n] * dvsor->cf[divisor];
191  }
192  }
193  for (n=1; n<=(rem->dgr+1); ++n) {
194  rem->cf[n-1] = quo->cf[quo->dgr+n];
195  quo->cf[quo->dgr+n] = 0;
196  }
197 }
198 
199 
200 int
201 bn_poly_quadratic_roots(register struct bn_complex *roots, register const struct bn_poly *quadrat)
202 {
203  fastf_t discrim, denom, rad;
204  const fastf_t small = SMALL_FASTF;
205 
206  if (NEAR_ZERO(quadrat->cf[0], small)) {
207  /* root = -cf[2] / cf[1] */
208  if (NEAR_ZERO(quadrat->cf[1], small)) {
209  /* No solution. Now what? */
210  /* bu_log("bn_poly_quadratic_roots(): ERROR, no solution\n"); */
211  return 0;
212  }
213  /* Fake it as a repeated root. */
214  roots[0].re = roots[1].re = -quadrat->cf[2]/quadrat->cf[1];
215  roots[0].im = roots[1].im = 0.0;
216  return 1; /* OK - repeated root */
217  }
218 
219  discrim = quadrat->cf[1]*quadrat->cf[1] - 4.0* quadrat->cf[0]*quadrat->cf[2];
220  denom = 0.5 / quadrat->cf[0];
221 
222  if (discrim > 0.0) {
223  rad = sqrt(discrim);
224 
225  if (NEAR_ZERO(quadrat->cf[1], small)) {
226  double r = fabs(rad * denom);
227  roots[0].re = r;
228  roots[1].re = -r;
229  } else {
230  double t, r1, r2;
231 
232  if (quadrat->cf[1] > 0.0) {
233  t = -0.5 * (quadrat->cf[1] + rad);
234  } else {
235  t = -0.5 * (quadrat->cf[1] - rad);
236  }
237  r1 = t / quadrat->cf[0];
238  r2 = quadrat->cf[2] / t;
239 
240  if (r1 < r2) {
241  roots[0].re = r1;
242  roots[1].re = r2;
243  } else {
244  roots[0].re = r2;
245  roots[1].re = r1;
246  }
247  }
248  roots[1].im = roots[0].im = 0.0;
249  } else if (NEAR_ZERO(discrim, small)) {
250  roots[1].re = roots[0].re = -quadrat->cf[1] * denom;
251  roots[1].im = roots[0].im = 0.0;
252  } else {
253  roots[1].re = roots[0].re = -quadrat->cf[1] * denom;
254  roots[1].im = -(roots[0].im = sqrt(-discrim) * denom);
255  }
256  return 1; /* OK */
257 }
258 
259 
260 int
261 bn_poly_cubic_roots(register struct bn_complex *roots, register const struct bn_poly *eqn)
262 {
263  fastf_t a, b, c1, c1_3rd, delta;
264  register int i;
265 
266  c1 = eqn->cf[1];
267  if (fabs(c1) > SQRT_MAX_FASTF) return 0; /* FAIL */
268 
269  c1_3rd = c1 * THIRD;
270  a = eqn->cf[2] - c1*c1_3rd;
271  if (fabs(a) > SQRT_MAX_FASTF) return 0; /* FAIL */
272  b = (2.0*c1*c1*c1 - 9.0*c1*eqn->cf[2] + 27.0*eqn->cf[3])*TWENTYSEVENTH;
273  if (fabs(b) > SQRT_MAX_FASTF) return 0; /* FAIL */
274 
275  if ((delta = a*a) > SQRT_MAX_FASTF) return 0; /* FAIL */
276  delta = b*b*0.25 + delta*a*TWENTYSEVENTH;
277 
278  if (delta > 0.0) {
279  fastf_t r_delta, A, B;
280 
281  r_delta = sqrt(delta);
282  A = B = -0.5 * b;
283  A += r_delta;
284  B -= r_delta;
285 
286  A = CUBEROOT(A);
287  B = CUBEROOT(B);
288 
289  roots[2].re = roots[1].re = -0.5 * (roots[0].re = A + B);
290 
291  roots[0].im = 0.0;
292  roots[2].im = -(roots[1].im = (A - B)*M_SQRT3*0.5);
293  } else if (ZERO(delta)) {
294  fastf_t b_2;
295  b_2 = -0.5 * b;
296 
297  roots[0].re = 2.0* CUBEROOT(b_2);
298  roots[2].re = roots[1].re = -0.5 * roots[0].re;
299  roots[2].im = roots[1].im = roots[0].im = 0.0;
300  } else {
301  fastf_t phi, fact;
302  fastf_t cs_phi, sn_phi_s3;
303 
304  if (a >= 0.0) {
305  fact = 0.0;
306  cs_phi = 1.0; /* cos(phi); */
307  sn_phi_s3 = 0.0; /* sin(phi) * M_SQRT3; */
308  } else {
309  register fastf_t f;
310  a *= -THIRD;
311  fact = sqrt(a);
312  if ((f = b * (-0.5) / (a*fact)) >= 1.0) {
313  cs_phi = 1.0; /* cos(phi); */
314  sn_phi_s3 = 0.0; /* sin(phi) * M_SQRT3; */
315  } else if (f <= -1.0) {
316  phi = M_PI_3;
317  cs_phi = cos(phi);
318  sn_phi_s3 = sin(phi) * M_SQRT3;
319  } else {
320  phi = acos(f) * THIRD;
321  cs_phi = cos(phi);
322  sn_phi_s3 = sin(phi) * M_SQRT3;
323  }
324  }
325 
326  roots[0].re = 2.0*fact*cs_phi;
327  roots[1].re = fact*(sn_phi_s3 - cs_phi);
328  roots[2].re = fact*(-sn_phi_s3 - cs_phi);
329  roots[2].im = roots[1].im = roots[0].im = 0.0;
330  }
331  for (i=0; i < 3; ++i)
332  roots[i].re -= c1_3rd;
333 
334  return 1; /* OK */
335 }
336 
337 
338 int
339 bn_poly_quartic_roots(register struct bn_complex *roots, register const struct bn_poly *eqn)
340 {
341  struct bn_poly cube, quad1, quad2;
342  bn_complex_t u[3];
343  fastf_t U, p, q, q1, q2;
344 
345  /* something considerably larger than squared floating point fuss */
346  const fastf_t small = 1.0e-8;
347 
348 #define Max3(a, b, c) ((c)>((a)>(b)?(a):(b)) ? (c) : ((a)>(b)?(a):(b)))
349 
350  cube.dgr = 3;
351  cube.cf[0] = 1.0;
352  cube.cf[1] = -eqn->cf[2];
353  cube.cf[2] = eqn->cf[3]*eqn->cf[1]
354  - 4*eqn->cf[4];
355  cube.cf[3] = -eqn->cf[3]*eqn->cf[3]
356  - eqn->cf[4]*eqn->cf[1]*eqn->cf[1]
357  + 4*eqn->cf[4]*eqn->cf[2];
358 
359  if (!bn_poly_cubic_roots(u, &cube)) {
360  return 0; /* FAIL */
361  }
362  if (!ZERO(u[1].im)) {
363  U = u[0].re;
364  } else {
365  U = Max3(u[0].re, u[1].re, u[2].re);
366  }
367 
368  p = eqn->cf[1]*eqn->cf[1]*0.25 + U - eqn->cf[2];
369  U *= 0.5;
370  q = U*U - eqn->cf[4];
371  if (p < 0) {
372  if (p < -small) {
373  return 0; /* FAIL */
374  }
375  p = 0;
376  } else {
377  p = sqrt(p);
378  }
379  if (q < 0) {
380  if (q < -small) {
381  return 0; /* FAIL */
382  }
383  q = 0;
384  } else {
385  q = sqrt(q);
386  }
387 
388  quad1.dgr = quad2.dgr = 2;
389  quad1.cf[0] = quad2.cf[0] = 1.0;
390  quad1.cf[1] = eqn->cf[1]*0.5;
391  quad2.cf[1] = quad1.cf[1] + p;
392  quad1.cf[1] -= p;
393 
394  q1 = U - q;
395  q2 = U + q;
396 
397  p = quad1.cf[1]*q2 + quad2.cf[1]*q1 - eqn->cf[3];
398  if (NEAR_ZERO(p, small)) {
399  quad1.cf[2] = q1;
400  quad2.cf[2] = q2;
401  } else {
402  q = quad1.cf[1]*q1 + quad2.cf[1]*q2 - eqn->cf[3];
403  if (NEAR_ZERO(q, small)) {
404  quad1.cf[2] = q2;
405  quad2.cf[2] = q1;
406  } else {
407  return 0; /* FAIL */
408  }
409  }
410 
411  bn_poly_quadratic_roots(&roots[0], &quad1);
412  bn_poly_quadratic_roots(&roots[2], &quad2);
413  return 1; /* SUCCESS */
414 }
415 
416 
417 void
418 bn_pr_poly(const char *title, register const struct bn_poly *eqn)
419 {
420  register size_t n;
421  register size_t exponent;
422  struct bu_vls str = BU_VLS_INIT_ZERO;
423  char buf[48];
424 
425  bu_vls_extend(&str, 196);
426  bu_vls_strcat(&str, title);
427  snprintf(buf, 48, " polynomial, degree = %lu\n", (long unsigned int)eqn->dgr);
428  bu_vls_strcat(&str, buf);
429 
430  exponent = eqn->dgr;
431  for (n=0; n<=eqn->dgr; n++, exponent--) {
432  register double coeff = eqn->cf[n];
433  if (n > 0) {
434  if (coeff < 0) {
435  bu_vls_strcat(&str, " - ");
436  coeff = -coeff;
437  } else {
438  bu_vls_strcat(&str, " + ");
439  }
440  }
441  bu_vls_printf(&str, "%g", coeff);
442  if (exponent > 1) {
443  bu_vls_printf(&str, " *X^%zu", exponent);
444  } else if (exponent == 1) {
445 
446  bu_vls_strcat(&str, " *X");
447  } else {
448  /* For constant term, add nothing */
449  }
450  }
451  bu_vls_strcat(&str, "\n");
452  bu_log("%s", bu_vls_addr(&str));
453  bu_vls_free(&str);
454 }
455 
456 
457 void
458 bn_pr_roots(const char *title, const struct bn_complex *roots, int n)
459 {
460  register int i;
461 
462  bu_log("%s: %d roots:\n", title, n);
463  for (i=0; i<n; i++) {
464  bu_log("%4d %e + i * %e\n", i, roots[i].re, roots[i].im);
465  }
466 }
467 
468 
469 /** @} */
470 /*
471  * Local Variables:
472  * mode: C
473  * tab-width: 8
474  * indent-tabs-mode: t
475  * c-file-style: "stroustrup"
476  * End:
477  * ex: shiftwidth=4 tabstop=8
478  */
fastf_t cf[BN_MAX_POLY_DEGREE+1]
Definition: poly.h:50
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
struct bn_poly * bn_poly_scale(register struct bn_poly *eqn, double factor)
Definition: poly.c:99
int bn_poly_quartic_roots(register struct bn_complex *roots, register const struct bn_poly *eqn)
Definition: poly.c:339
#define BN_MAX_POLY_DEGREE
Library for dealing with polynomials.
Definition: poly.h:41
void bu_vls_strcat(struct bu_vls *vp, const char *s)
Definition: vls.c:368
struct bn_poly * bn_poly_sub(register struct bn_poly *diff, register const struct bn_poly *poly1, register const struct bn_poly *poly2)
Definition: poly.c:140
void bn_poly_synthetic_division(register struct bn_poly *quo, register struct bn_poly *rem, register const struct bn_poly *dvdend, register const struct bn_poly *dvsor)
Definition: poly.c:171
#define SMALL_FASTF
Definition: defines.h:342
Header file for the BRL-CAD common definitions.
Definition: poly.h:47
struct bn_poly * bn_poly_add(register struct bn_poly *sum, register const struct bn_poly *poly1, register const struct bn_poly *poly2)
Definition: poly.c:111
#define BN_POLY_MAGIC
Definition: magic.h:70
int bn_poly_quadratic_roots(register struct bn_complex *roots, register const struct bn_poly *quadrat)
Definition: poly.c:201
double im
imaginary part
Definition: complex.h:41
struct bn_poly * bn_poly_mul(register struct bn_poly *product, register const struct bn_poly *m1, register const struct bn_poly *m2)
Definition: poly.c:50
#define CUBEROOT(a)
Definition: poly.c:41
double re
real part
Definition: complex.h:40
#define Max3(a, b, c)
void bu_vls_free(struct bu_vls *vp)
Definition: vls.c:248
size_t dgr
Definition: poly.h:49
bn_poly_t quo[1]
void bn_pr_poly(const char *title, register const struct bn_poly *eqn)
Definition: poly.c:418
void bn_pr_roots(const char *title, const struct bn_complex *roots, int n)
Definition: poly.c:458
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#define BN_POLY_NULL
Definition: poly.h:54
char * bu_vls_addr(const struct bu_vls *vp)
Definition: vls.c:111
#define ZERO(val)
Definition: units.c:38
void bu_vls_printf(struct bu_vls *vls, const char *fmt,...) _BU_ATTR_PRINTF23
Definition: vls.c:694
void bu_vls_extend(struct bu_vls *vp, size_t extra)
Definition: vls.c:136
int bn_poly_cubic_roots(register struct bn_complex *roots, register const struct bn_poly *eqn)
Definition: poly.c:261
#define A
Definition: msr.c:51
#define SQRT_MAX_FASTF
Definition: defines.h:341
Complex numbers.
Definition: complex.h:39
#define BU_VLS_INIT_ZERO
Definition: vls.h:84
bn_poly_t rem[1]
Definition: vls.h:56
HIDDEN const point_t delta
Definition: sh_prj.c:618
double fastf_t
Definition: defines.h:300