RTXI  2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
laguerre.cpp
Go to the documentation of this file.
1 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 //
3 // File = laguerre.cpp
4 //
5 // Laguerre method for finding polynomial roots
6 //
7 
8 #include "laguerre.h"
9 #include "complex.h"
10 #include <fstream>
11 #include <iostream>
12 #include <math.h>
13 #include <stdlib.h>
14 
15 #ifdef _DEBUG
16 extern std::ofstream DebugFile;
17 #endif
18 
19 int
20 LaguerreMethod(CmplxPolynomial* poly, complex* root_ptr, double epsilon,
21  double epsilon2, int max_iter)
22 {
23  int iter, j, order;
24  complex p_eval, p_prime, p_double_prime;
25  complex root, f, f_sqrd, g, radical;
26  complex f_plus_rad, f_minus_rad, delta_root;
27  double old_delta_mag, root_mag, error;
28  complex* coeff;
29  order = poly->GetDegree();
30  coeff = new complex[order + 1];
31  poly->CopyCoeffs(coeff);
32  root = *root_ptr;
33  old_delta_mag = 1000.0;
34 
35  for (iter = 1; iter <= max_iter; iter++) {
36  p_double_prime = complex(0.0, 0.0);
37  p_prime = complex(0.0, 0.0);
38  p_eval = coeff[order];
39  error = mag(p_eval);
40  root_mag = mag(root);
41 
42  for (j = order - 1; j >= 0; j--) {
43  p_double_prime = p_prime + root * p_double_prime;
44  p_prime = p_eval + root * p_prime;
45  p_eval = coeff[j] + root * p_eval;
46  error = mag(p_eval) + root_mag * error;
47  }
48  error = epsilon2 * error;
49  p_double_prime = 2.0 * p_double_prime;
50 
51  if (mag(p_eval) < error) {
52  std::cout << "mag(p_eval) = " << mag(p_eval) << " error = " << error
53  << std::endl;
54  *root_ptr = root;
55  delete[] coeff;
56  return (1);
57  }
58  f = p_prime / p_eval;
59  f_sqrd = f * f;
60  g = f_sqrd - p_double_prime / p_eval;
61  radical = (order - 1) * (order * g - f_sqrd);
62  radical = sqrt(radical);
63  f_plus_rad = f + radical;
64  f_minus_rad = f - radical;
65  if (mag(f_plus_rad) > mag(f_minus_rad)) {
66  delta_root = complex(double(order), 0.0) / f_plus_rad;
67  } else {
68  delta_root = complex(double(order), 0.0) / f_minus_rad;
69  }
70  root = root - delta_root;
71  if ((iter > 6) && (mag(delta_root) > old_delta_mag)) {
72  *root_ptr = root;
73  delete[] coeff;
74  return (2);
75  }
76  if (mag(delta_root) < (epsilon * mag(root))) {
77  *root_ptr = root;
78  delete[] coeff;
79  return (3);
80  }
81  old_delta_mag = mag(delta_root);
82  }
83 #ifdef _DEBUG
84  DebugFile << "Laguerre method failed to converge" << std::endl;
85 #endif
86  delete[] coeff;
87  return (-1);
88 }
double mag(const complex _z)
Definition: complex.h:123
void CopyCoeffs(complex *)
Definition: cmpxpoly.cpp:334
complex sqrt(const complex _z)
Definition: complex.h:148
std::ofstream DebugFile
int LaguerreMethod(CmplxPolynomial *poly, complex *root_ptr, double epsilon, double epsilon2, int max_iter)
Definition: laguerre.cpp:20
int GetDegree(void)
Definition: cmpxpoly.cpp:316