BRL-CAD
bn_poly_quartic_roots.c
Go to the documentation of this file.
1 /* B N _ P O L Y _ Q U A R T 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 /**
43  * Initialises quartic equations storing a negative, positive and zero coefficients.
44  */
45 void
46 poly_init(void)
47 {
48 
49  /* initializes a zero equation */
50 
51  input[0] = bn_Zero_poly;
52  input[0].dgr = 4;
53  input[0].cf[0] = input[0].cf[1] = input[0].cf[2] = input[0].cf[3] = input[0].cf[4] = 0.0;
54 
55  rts[0].re = rts[0].im = 0.0;
56 
57 
58  /* initializes a negative quadratic eqn. */
59  input[1] = bn_Zero_poly;
60  input[1].dgr = 4;
61 
62 
63  input[1].cf[0] = -4, input[1].cf[1] = -3, input[1].cf[2] = -2, input[1].cf[3] = -25, input[1].cf[4] = -38;/* input coeff */
64 
65  /**
66  * The known output values used for these tests were generated from
67  * GNU Octave, version 3.4.3
68  */
69  rts[1].re = 0.2613656082942032, rts[1].im = 0.4284631324677022;
70  rts[2].re = 0.2613656082942032, rts[2].im = -0.4284631324677022;
71  rts[3].re = -0.5903129767152558, rts[3].im = 0.263475942656035;
72  rts[4].re = -0.5903129767152558, rts[4].im = -0.263475942656035;
73 
74  /* initializes a positive quadratic equation */
75  input[2] = bn_Zero_poly;
76  input[2].dgr = 4;
77  input[2].cf[0] = 5478, input[2].cf[1] = 5485, input[2].cf[2] = 458, input[2].cf[3] = 258564, input[2].cf[4] = 54785;/* input coeff */
78 
79  rts[5].re = 0.12889648467110737 , rts[5].im = 0.25711127015404556;
80  rts[6].re = 0.12889648467110737, rts[6].im = -0.25711127015404556;
81  rts[7].re = -0.25602234520349354, rts[7].im = 0.0;
82  rts[8].re = -4.721383656903164, rts[8].im = 0.0;
83 
84  return;
85 }
86 
87 
88 /* tests the polynomials to make sure bn_poly_mul() works properly. */
89 int
91 {
92  int val, val1[4], val2[4];/* variables get results for comparisons */
93  bn_complex_t r1, r2[4], r3[4];
94 
95  int i, j, ind1 = 1, ind2 = 5; /* creates array indexes for rts */
96 
97  bn_poly_quartic_roots(&r1, &input[0]);
98  bn_poly_quartic_roots(r2, &input[1]);
99  bn_poly_quartic_roots(r3, &input[2]);
100 
101  /* loop moves through arrays for roots and output comparing variables. */
102  for (i = 0; i < 4; i++) {
103 
104  /*checks r1 */
105  if (EQUAL(r1.re, rts[0].re) && EQUAL(r1.im, rts[0].im))
106  val = 0;
107 
108  /*checks r2 */
109  for (j = ind1; j< ind1 + 4; j++) {
110  if (EQUAL(r2[i].re, rts[j].re) && EQUAL(r2[i].im, rts[j].im)) {
111  val1[i] = 0;
112  continue;
113  }
114  }
115 
116  /*check r3 */
117  for (j = ind2; j< ind2 + 4; j++) {
118  if (EQUAL(r3[i].re, rts[j].re) && EQUAL(r3[i].im, rts[j].im)) {
119  val2[i] = 0;
120  continue;
121  }
122  }
123  }
124 
125  if (!EQUAL(val, 0))
126  return -1;
127 
128  for (i = 0; i < 4; i++) {
129  if (!EQUAL(val1[i], val2[i]) && !EQUAL(val1[i], 0))
130  return -1;
131  }
132 
133  return 0;
134 }
135 
136 
137 int
138 main(void)
139 {
140  int ret;
141  poly_init();
142  ret = test_bn_poly_qua_rts();
143 
144  if (ret == 0) {
145  bu_log("\nFunction computes correctly\n");
146  return ret;
147 
148  } else
149  exit(EXIT_FAILURE);
150 
151  return 0;
152 
153 }
154 
155 
156 /*
157  * Local Variables:
158  * mode: C
159  * tab-width: 8
160  * indent-tabs-mode: t
161  * c-file-style: "stroustrup"
162  * End:
163  * ex: shiftwidth=4 tabstop=8
164  */
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
void poly_init(void)
struct bn_poly bn_Zero_poly
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...
int test_bn_poly_qua_rts(void)
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]
bn_complex_t rts[9]
Complex numbers.
Definition: complex.h:39
int main(void)