BRL-CAD
bn_poly_cubic_roots.c
Go to the documentation of this file.
1 /* B N _ P O L Y _ C U B I C _ R O O T S. C
2  * BRL-CAD
3  *
4  * Copyright (c) 2013-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 <stdio.h>
25 #include <stdlib.h>
26 #include <math.h>
27 #include <string.h>
28 #include <signal.h>
29 
30 #include "bu.h"
31 #include "vmath.h"
32 #include "bn.h"
33 
34 
35 /* holds three polynomials to be used in test. */
38 
39 struct bn_poly bn_Zero_poly = { BN_POLY_MAGIC, 0, {0.0} };
40 
41 /**
42  * Initialises cubic equations storing a negative, positive and zero coefficients.
43  */
44 void
45 poly_init(void)
46 {
47 
48  /* initializes a zero equation */
49 
50  input[0] = bn_Zero_poly;
51  input[0].dgr = 3;
52  input[0].cf[0] = input[0].cf[1] = input[0].cf[2] = input[0].cf[3] = 0.0;
53 
54  rts[0].re = rts[0].im = 0.0;
55 
56 
57  /* initializes a negative cubic eqn. */
58  input[1] = bn_Zero_poly;
59  input[1].dgr = 3;
60 
61 
62  input[1].cf[0] = -4, input[1].cf[1] = -3, input[1].cf[2] = -2, input[1].cf[3] = -25; /* input coeff */
63 
64  /**
65  * The known output values used for these tests were generated from
66  * GNU Octave, version 3.4.3
67  */
68  rts[1].re = -0.49359, rts[1].im = 0.0;
69  rts[2].re = 0.20679876865588492, rts[2].im = 0.5304573452575734;
70  rts[3].re = 0.20679876865588492 , rts[3].im = -0.5304573452575734;
71 
72  /* initializes a positive cubic equation */
73  input[2] = bn_Zero_poly;
74  input[2].dgr = 3;
75  input[2].cf[0] = 5478, input[2].cf[1] = 5485, input[2].cf[2] = 458, input[2].cf[3] = 786;/* input coeff */
76 
77  rts[4].re = -0.9509931181746001, rts[4].im = 0.0;
78  rts[5].re = 0.18414795857839417 , rts[5].im = 2.700871695081346;
79  rts[6].re = 0.18414795857839417 , rts[6].im = -2.700871695081346;
80 
81  return;
82 }
83 
84 
85 /* tests the polynomials to make sure bn_poly_mul() works properly. */
86 int
88 {
89  int val, val1[3], val2[3];/* variables get results for comparisons */
90  bn_complex_t r1[3], r2[3], r3[3];
91  int i, j, ind1 = 1, ind2 = 4; /* creates indexes for rts array. */
92 
93  bn_poly_cubic_roots(r1, &input[0]);
94  bn_poly_cubic_roots(r2, &input[1]);
95  bn_poly_cubic_roots(r3, &input[2]);
96 
97  /* loop moves through arrays for roots and output comparing variables. */
98  for (i = 0; i < 3; i++) {
99 
100  /*checks r1 */
101  if (EQUAL(r1[i].re, rts[0].re) && EQUAL(r1[i].im, rts[0].im))
102  val = 0;
103 
104  /*checks r2 */
105  for (j = ind1; j< ind1 + 3; j++) {
106  if (EQUAL(r2[i].re, rts[j].re) && EQUAL(r2[i].im, rts[j].im)) {
107  val1[i] = 0;
108  continue;
109  }
110  }
111 
112  /*check r3 */
113  for (j = ind2; j< ind2 + 3; j++) {
114  if (EQUAL(r3[i].re, rts[j].re) && EQUAL(r3[i].im, rts[j].im)) {
115  val2[i] = 0;
116  continue;
117  }
118  }
119  }
120 
121  if (!EQUAL(val, 0))
122  return -1;
123 
124  for (i = 0; i < 3; i++) {
125  if (!EQUAL(val1[i], val2[i]) && !EQUAL(val1[i], 0))
126  return -1;
127  }
128 
129  return 0;
130 }
131 
132 
133 int
134 main(void)
135 {
136  int ret;
137  poly_init();
138  ret = test_bn_poly_cubic_rts();
139 
140  if (ret == 0) {
141  bu_log("\nFunction computes correctly\n");
142  return ret;
143  } else
144  exit(EXIT_FAILURE);
145 
146  return 0;
147 }
148 
149 
150 /*
151  * Local Variables:
152  * mode: C
153  * tab-width: 8
154  * indent-tabs-mode: t
155  * c-file-style: "stroustrup"
156  * End:
157  * ex: shiftwidth=4 tabstop=8
158  */
Definition: db_flip.c:35
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_Zero_poly
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
double re
real part
Definition: complex.h:40
size_t dgr
Definition: poly.h:49
bn_poly_t input[3]
void poly_init(void)
int main(void)
int test_bn_poly_cubic_rts(void)
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...
bn_complex_t rts[7]
Complex numbers.
Definition: complex.h:39