CbmRoot
PolynomComplexRoots.h
Go to the documentation of this file.
1 #include <iomanip>
2 #include <iostream>
3 #include <math.h>
4 #include <stdlib.h>
5 
6 #define maxiter 500
7 
8 //
9 // Extract individual real or complex roots from list of quadratic factors
10 //
11 int roots(float* a, int n, float* wr, float* wi) {
12  float sq, b2, c, disc;
13  int m, numroots;
14 
15  m = n;
16  numroots = 0;
17  while (m > 1) {
18  b2 = -0.5 * a[m - 2];
19  c = a[m - 1];
20  disc = b2 * b2 - c;
21  if (disc < 0.0) { // complex roots
22  sq = sqrt(-disc);
23  wr[m - 2] = b2;
24  wi[m - 2] = sq;
25  wr[m - 1] = b2;
26  wi[m - 1] = -sq;
27  numroots += 2;
28  } else { // real roots
29  sq = sqrt(disc);
30  wr[m - 2] = fabs(b2) + sq;
31  if (b2 < 0.0) wr[m - 2] = -wr[m - 2];
32  if (wr[m - 2] == 0)
33  wr[m - 1] = 0;
34  else {
35  wr[m - 1] = c / wr[m - 2];
36  numroots += 2;
37  }
38  wi[m - 2] = 0.0;
39  wi[m - 1] = 0.0;
40  }
41  m -= 2;
42  }
43  if (m == 1) {
44  wr[0] = -a[0];
45  wi[0] = 0.0;
46  numroots++;
47  }
48  return numroots;
49 }
50 //
51 // Deflate polynomial 'a' by dividing out 'quad'. Return quotient
52 // polynomial in 'b' and error metric based on remainder in 'err'.
53 //
54 void deflate(float* a, int n, float* b, float* quad, float* err) {
55  float r, s;
56  int i;
57 
58  r = quad[1];
59  s = quad[0];
60 
61  b[1] = a[1] - r;
62 
63  for (i = 2; i <= n; i++) {
64  b[i] = a[i] - r * b[i - 1] - s * b[i - 2];
65  }
66  *err = fabs(b[n]) + fabs(b[n - 1]);
67 }
68 //
69 // Find quadratic factor using Bairstow's method (quadratic Newton method).
70 // A number of ad hoc safeguards are incorporated to prevent stalls due
71 // to common difficulties, such as zero slope at iteration point, and
72 // convergence problems.
73 //
74 // Bairstow's method is sensitive to the starting estimate. It is possible
75 // for convergence to fail or for 'wild' values to trigger an overflow.
76 //
77 // It is advisable to institute traps for these problems. (To do!)
78 //
79 void find_quad(float* a, int n, float* b, float* quad, float* err, int* iter) {
80  float *c, dn, dr, ds, drn, dsn, eps, r, s;
81  int i;
82 
83  c = new float[n + 1];
84  c[0] = 1.0;
85  r = quad[1];
86  s = quad[0];
87  eps = 1e-15;
88  *iter = 1;
89 
90  do {
91  if (*iter > maxiter) break;
92  if (((*iter) % 200) == 0) { eps *= 10.0; }
93  b[1] = a[1] - r;
94  c[1] = b[1] - r;
95 
96  for (i = 2; i <= n; i++) {
97  b[i] = a[i] - r * b[i - 1] - s * b[i - 2];
98  c[i] = b[i] - r * c[i - 1] - s * c[i - 2];
99  }
100  dn = c[n - 1] * c[n - 3] - c[n - 2] * c[n - 2];
101  drn = b[n] * c[n - 3] - b[n - 1] * c[n - 2];
102  dsn = b[n - 1] * c[n - 1] - b[n] * c[n - 2];
103 
104  if (fabs(dn) < 1e-10) {
105  if (dn < 0.0)
106  dn = -1e-8;
107  else
108  dn = 1e-8;
109  }
110  dr = drn / dn;
111  ds = dsn / dn;
112  r += dr;
113  s += ds;
114  (*iter)++;
115  } while ((fabs(dr) + fabs(ds)) > eps);
116  quad[0] = s;
117  quad[1] = r;
118  *err = fabs(ds) + fabs(dr);
119  delete[] c;
120 }
121 //
122 // Differentiate polynomial 'a' returning result in 'b'.
123 //
124 void diff_poly(float* a, int n, float* b) {
125  float coef;
126  int i;
127 
128  coef = (float) n;
129  b[0] = 1.0;
130  for (i = 1; i < n; i++) {
131  b[i] = a[i] * ((float) (n - i)) / coef;
132  }
133 }
134 //
135 // Attempt to find a reliable estimate of a quadratic factor using modified
136 // Bairstow's method with provisions for 'digging out' factors associated
137 // with multiple roots.
138 //
139 // This resursive routine operates on the principal that differentiation of
140 // a polynomial reduces the order of all multiple roots by one, and has no
141 // other roots in common with it. If a root of the differentiated polynomial
142 // is a root of the original polynomial, there must be multiple roots at
143 // that location. The differentiated polynomial, however, has lower order
144 // and is easier to solve.
145 //
146 // When the original polynomial exhibits convergence problems in the
147 // neighborhood of some potential root, a best guess is obtained and tried
148 // on the differentiated polynomial. The new best guess is applied
149 // recursively on continually differentiated polynomials until failure
150 // occurs. At this point, the previous polynomial is accepted as that with
151 // the least number of roots at this location, and its estimate is
152 // accepted as the root.
153 //
154 void recurse(float* a,
155  int n,
156  float* b,
157  int m,
158  float* quad,
159  float* err,
160  int* iter) {
161  float *c, *x, rs[2], tst;
162 
163  if (fabs(b[m]) < 1e-16) m--; // this bypasses roots at zero
164  if (m == 2) {
165  quad[0] = b[2];
166  quad[1] = b[1];
167  *err = 0;
168  *iter = 0;
169  return;
170  }
171  c = new float[m + 1];
172  x = new float[n + 1];
173  c[0] = x[0] = 1.0;
174  rs[0] = quad[0];
175  rs[1] = quad[1];
176  *iter = 0;
177  find_quad(b, m, c, rs, err, iter);
178  tst = fabs(rs[0] - quad[0]) + fabs(rs[1] - quad[1]);
179  if (*err < 1e-12) {
180  quad[0] = rs[0];
181  quad[1] = rs[1];
182  }
183  // tst will be 'large' if we converge to wrong root
184  if (((*iter > 5) && (tst < 1e-4)) || ((*iter > 20) && (tst < 1e-1))) {
185  diff_poly(b, m, c);
186  recurse(a, n, c, m - 1, rs, err, iter);
187  quad[0] = rs[0];
188  quad[1] = rs[1];
189  }
190  delete[] x;
191  delete[] c;
192 }
193 //
194 // Top level routine to manage the determination of all roots of the given
195 // polynomial 'a', returning the quadratic factors (and possibly one linear
196 // factor) in 'x'.
197 //
198 void get_quads(float* a, int n, float* quad, float* x) {
199  float *b, *z, err, tmp;
200  int iter, i, m;
201 
202  if ((tmp = a[0]) != 1.0) {
203  a[0] = 1.0;
204  for (i = 1; i <= n; i++) {
205  a[i] /= tmp;
206  }
207  }
208  if (n == 2) {
209  x[0] = a[1];
210  x[1] = a[2];
211  return;
212  } else if (n == 1) {
213  x[0] = a[1];
214  return;
215  }
216  m = n;
217  b = new float[n + 1];
218  z = new float[n + 1];
219  b[0] = 1.0;
220  for (i = 0; i <= n; i++) {
221  z[i] = a[i];
222  x[i] = 0.0;
223  }
224  do {
225  if (n > m) {
226  quad[0] = 3.14159e-1;
227  quad[1] = 2.78127e-1;
228  }
229  do { // This loop tries to assure convergence
230  for (i = 0; i < 5; i++) {
231  find_quad(z, m, b, quad, &err, &iter);
232  if ((err > 1e-7) || (iter > maxiter)) {
233  diff_poly(z, m, b);
234  recurse(z, m, b, m - 1, quad, &err, &iter);
235  }
236  deflate(z, m, b, quad, &err);
237  if (err < 0.001) break;
238  quad[0] = 9999.;
239  quad[1] = 9999.;
240  }
241  if (err > 0.01) {
242  printf("Error! Convergence failure in quadratic x^2 + r*x + s.\n");
243  exit(1);
244  }
245  } while (err > 0.01);
246  x[m - 2] = quad[1];
247  x[m - 1] = quad[0];
248  m -= 2;
249  for (i = 0; i <= m; i++) {
250  z[i] = b[i];
251  }
252  } while (m > 2);
253  if (m == 2) {
254  x[0] = b[1];
255  x[1] = b[2];
256  } else
257  x[0] = b[1];
258  delete[] z;
259  delete[] b;
260 }
261 
262 void polynomComplexRoots(float* wr, float* wi, int n, float* a, int& numr) {
263  float* quad = new float[2];
264  float* x = new float[n];
265 
266  // initialize estimate for 1st root pair
267  quad[0] = 2.71828e-1;
268  quad[1] = 3.14159e-1;
269 
270  // get roots
271  get_quads(a, n, quad, x);
272  numr = roots(x, n, wr, wi);
273 
274  delete[] quad;
275  delete[] x;
276 }
diff_poly
void diff_poly(float *a, int n, float *b)
Definition: PolynomComplexRoots.h:124
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
find_quad
void find_quad(float *a, int n, float *b, float *quad, float *err, int *iter)
Definition: PolynomComplexRoots.h:79
maxiter
#define maxiter
Definition: PolynomComplexRoots.h:6
polynomComplexRoots
void polynomComplexRoots(float *wr, float *wi, int n, float *a, int &numr)
Definition: PolynomComplexRoots.h:262
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
deflate
void deflate(float *a, int n, float *b, float *quad, float *err)
Definition: PolynomComplexRoots.h:54
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
recurse
void recurse(float *a, int n, float *b, int m, float *quad, float *err, int *iter)
Definition: PolynomComplexRoots.h:154
roots
int roots(float *a, int n, float *wr, float *wi)
Definition: PolynomComplexRoots.h:11
get_quads
void get_quads(float *a, int n, float *quad, float *x)
Definition: PolynomComplexRoots.h:198