RTXI  2.1
cmpxpoly.cpp
Go to the documentation of this file.
1 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 //
3 // File = cmpxpoly.cpp
4 //
5 // class that implements a polynomial with
6 // complex-valued coefficients
7 //
8
9 #include "cmpxpoly.h"
10 #include "laguerre.h"
11 #include "pause.h"
12 #include "stdlib.h"
13 #include <math.h>
14 #include <unistd.h>
15
16 //======================================================
17 // default constructor
18
20 {
21  Degree = 0;
22  Coeff = new complex[1];
23  Coeff[0] = complex(0.0, 0.0);
24  RemCoeff = NULL;
25  Root = NULL;
26  return;
27 };
28
29 //===================================================================
30 // copy constructor
31
33 {
34  int i;
35  Degree = original.Degree;
36  Coeff = new complex[Degree + 1];
37  Root = new complex[Degree];
38
39  for (i = 0; i <= Degree; i++)
40  Coeff[i] = original.Coeff[i];
41
42  for (i = 0; i < Degree; i++)
43  Root[i] = original.Root[i];
44
45  return;
46 };
47
48 //===================================================================
49 // constructor for initializing a binomial
50
51 CmplxPolynomial::CmplxPolynomial(const complex coeff_1, const complex coeff_0)
52 {
53  Degree = 1;
54  Coeff = new complex[2];
55  RemCoeff = NULL;
56  Root = NULL;
57
58  Coeff[0] = coeff_0;
59  Coeff[1] = coeff_1;
60
61  return;
62 }
63
64 //===================================================================
65 // initializing constructor
66
67 CmplxPolynomial::CmplxPolynomial(const complex* coeff, const int degree)
68 {
69  Degree = degree;
70  Coeff = new complex[degree + 1];
71  RemCoeff = NULL;
72  Root = NULL;
73
74  for (int i = 0; i <= degree; i++)
75  Coeff[i] = coeff[i];
76  for (int j = 0; j <= degree; j++) {
77  std::cout << "poly_coeff[" << j << "] = " << coeff[j] << std::endl;
78  }
79
80  return;
81 }
82
83 CmplxPolynomial::CmplxPolynomial(const double* coeff, const int degree)
84 {
85  Degree = degree;
86  Coeff = new complex[degree + 1];
87  RemCoeff = NULL;
88  Root = NULL;
89
90  for (int i = 0; i <= degree; i++)
91  Coeff[i] = complex(coeff[i], 0.0);
92
93  return;
94 }
95
96 //=================================================================
97 // assignment operator
98
101 {
102  if (Coeff != right.Coeff) {
103  //-------------------------------------------------------------
104  // Get rid of old coefficient array to make way for a new one
105  // of the correct length for the new polynomial being assigned
106
107  delete[] Coeff;
108  delete[] Root;
109
110  Degree = right.Degree;
111  Coeff = new complex[Degree + 1];
112
113  for (int i = 0; i <= Degree; i++)
114  Coeff[i] = right.Coeff[i];
115  }
116  return *this;
117 }
118
119 //===================================================================
120 // multiply assign operator
121
124 {
125  //-----------------------------------------------------
126  // save pointer to original coefficient array so that
127  // this array can be deleted once no longer needed
128
129  complex* orig_coeff = Coeff;
130  int orig_degree = Degree;
131
132  //-------------------------------------------------------
133  // create new longer array to hold the new coefficients
134
135  Degree += right.Degree;
136  Coeff = new complex[Degree + 1];
137
138  for (int i = 0; i <= Degree; i++)
139  Coeff[i] = complex(0.0, 0.0);
140
141  //---------------------------------
142  // perform multiplication
143
144  for (int rgt_indx = 0; rgt_indx <= right.Degree; rgt_indx++) {
145  for (int orig_indx = 0; orig_indx <= orig_degree; orig_indx++) {
146  Coeff[orig_indx + rgt_indx] +=
147  (orig_coeff[orig_indx] * right.Coeff[rgt_indx]);
148  }
149  }
150
151  return *this;
152 }
153
154 //===================================================================
155 // divide assign operator
156
159 {
160  //----------------------------------------------------
161  // In general, polynomial division will produce a
162  // quotient and a remainder. This routine returns the
163  // quotient as its result. The remainder will be
164  // stored in a member variable so that it can be
165  // checked or retrived by subsequent calls to the
166  // appropriate member functions.
167  //-----------------------------------------------------
168  // save pointer to original coefficient array so that
169  // this array can be deleted once no longer needed
170
171  // complex *orig_coeff = Coeff;
172  int dvdnd_deg, dvsr_deg, j, k;
173
174  //-------------------------------------------------------
175  // create new array to hold the new coefficients
176
177  if (RemCoeff == NULL)
178  RemCoeff = new complex[Degree + 1];
179  dvdnd_deg = Degree;
180  dvsr_deg = divisor.Degree;
181  Degree -= dvsr_deg;
182
183  //---------------------------------
184  // perform division
185
186  for (j = 0; j <= dvdnd_deg; j++) {
187  RemCoeff[j] = Coeff[j];
188  Coeff[j] = complex(0.0, 0.0);
189  }
190  for (k = dvdnd_deg - dvsr_deg; k >= 0; k--) {
191  Coeff[k] = RemCoeff[dvsr_deg + k] / divisor.Coeff[dvsr_deg];
192  for (j = dvsr_deg + k - 1; j >= k; j--)
193  RemCoeff[j] -= Coeff[k] * divisor.Coeff[j - k];
194  }
195  for (j = dvsr_deg; j <= dvdnd_deg; j++)
196  RemCoeff[j] = complex(0.0, 0.0);
197  return *this;
198 }
199 //========================================================
200 // Find roots of polynomial
201
202 void
204 {
205  complex* root;
206  int status, i;
207  CmplxPolynomial root_factor;
208  // root = new complex[Degree]; // What's this supposed to do? -ansel
209  CmplxPolynomial work_poly;
210  double epsilon = 0.0000001;
211  double epsilon2 = 1.0e-10;
212  int max_iter = 12;
213
214  if (Root == NULL)
215  Root = new complex[Degree];
216  //------------------------------------------------
217  // find coarse locations for roots
218
219  work_poly = CmplxPolynomial(Coeff, Degree);
220
221  for (i = 0; i < Degree - 1; i++) {
222  Root[i] = complex(0.0, 0.0);
223  status =
224  LaguerreMethod(&work_poly, &(Root[i]), epsilon, epsilon2, max_iter);
225  if (status < 0) {
226  std::cout << "Laguerre method did not converge" << std::endl;
227  exit(55);
228  }
229  std::cout << "Root[" << i << "] = " << Root[i] << " (" << status << ")"
230  << std::endl;
231  root_factor = CmplxPolynomial(complex(1.0, 0.0), -Root[i]);
232  work_poly /= root_factor;
233  work_poly.DumpToStream(&std::cout);
234  pausewait();
235  }
236  Root[Degree - 1] = -(work_poly.GetCoeff(0));
237  std::cout << "Root[" << Degree - 1 << "] = " << Root[Degree - 1] << std::endl;
238
239  //------------------------------------------------
240  // polish the roots
241  work_poly = CmplxPolynomial(Coeff, Degree);
242  for (i = 0; i < Degree; i++) {
243  status =
244  LaguerreMethod(&work_poly, &(Root[i]), epsilon, epsilon2, max_iter);
245  if (status < 0) {
246  std::cout << "Laguerre method did not converge" << std::endl;
247  exit(55);
248  }
249  std::cout << "Polished Root[" << i << "] = " << Root[i] << " (" << status
250  << ")" << std::endl;
251  pause();
252  }
253  return;
254 }
255
256 //========================================================
257 // Get array of polynomial root values
258
259 complex*
261 {
262  complex* root;
263  int i;
264  root = new complex[Degree];
265  if (Root == NULL)
266  this->FindRoots();
267
268  for (i = 0; i < Degree; i++)
269  root[i] = Root[i];
270  return (root);
271 }
272
273 //=========================================================
274 // reflect root across the unit circle
275
276 void
278 {
279  // (1) we must first find the roots (with more accuracy than
280  // Laguerre alone provides) then (2) replace the root of
281  // interest with its reciprocal and (3) reconstitute
282  // the sum-of-powers coefficients from the revised set
283  // of roots.
284  this->FindRoots();
285  Root[root_idx] = complex(1.0, 0.0) / Root[root_idx];
286  this->BuildFromRoots();
287 }
288 //=========================================================
289 // build coefficients from the roots
290
291 void
293 {
294  if (Root != NULL) // do nothing if roots not defined
295  {
296  }
297 }
298 //=========================================================
299 // dump polynomial to an output stream
300
301 void
302 CmplxPolynomial::DumpToStream(std::ostream* output_stream)
303 {
304  (*output_stream) << "Degree = " << Degree << std::endl;
305
306  for (int i = Degree; i >= 0; i--) {
307  (*output_stream) << "Coeff[" << i << "] = " << Coeff[i] << std::endl;
308  }
309  return;
310 }
311
312 //====================================================
313 //
314
315 int
317 {
318  return (Degree);
319 }
320
321 //==================================================
322 //
323
324 complex
326 {
327  return Coeff[k];
328 }
329
330 //==================================================
331 //
332
333 void
335 {
336  for (int i = 0; i <= Degree; i++)
337  coeff[i] = Coeff[i];
338  return;
339 }
340
341 //\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$\$
CmplxPolynomial & operator*=(const CmplxPolynomial &right)
Definition: cmpxpoly.cpp:123
CmplxPolynomial & operator/=(const CmplxPolynomial &divisor)
Definition: cmpxpoly.cpp:158
void FindRoots(void)
Definition: cmpxpoly.cpp:203
complex * Root
Definition: cmpxpoly.h:67
void ReflectRoot(int root_idx)
Definition: cmpxpoly.cpp:277
void BuildFromRoots(void)
Definition: cmpxpoly.cpp:292
void CopyCoeffs(complex *)
Definition: cmpxpoly.cpp:334
complex * Coeff
Definition: cmpxpoly.h:65
CmplxPolynomial & operator=(const CmplxPolynomial &right)
Definition: cmpxpoly.cpp:100
void DumpToStream(std::ostream *output_stream)
Definition: cmpxpoly.cpp:302
int LaguerreMethod(CmplxPolynomial *poly, complex *root_ptr, double epsilon, double epsilon2, int max_iter)
Definition: laguerre.cpp:20
void pausewait(void)
Definition: pause.cpp:9
int GetDegree(void)
Definition: cmpxpoly.cpp:316
complex * GetRoots(void)
Definition: cmpxpoly.cpp:260
complex GetCoeff(int k)
Definition: cmpxpoly.cpp:325
void pause(logical pause_enabled)
Definition: fs_util.cpp:24
complex * RemCoeff
Definition: cmpxpoly.h:66