VERB_code_2.3
rroots.cpp
Go to the documentation of this file.
1 
12 #include "rroots.h"
13 
15 int roots(double *a,int n,double *wr,double *wi)
16 {
17  double sq,b2,c,disc;
18  int m,numroots;
19 
20  m = n;
21  numroots = 0;
22  while (m > 1) {
23  b2 = -0.5*a[m-2];
24  c = a[m-1];
25  disc = b2*b2-c;
26  if (fabs(disc)/(fabs(b2*b2)+fabs(c)) <= DBL_EPSILON) disc = 0.0;
27  if (disc < 0.0) {
28  sq = sqrt(-disc);
29  wr[m-2] = b2;
30  wi[m-2] = sq;
31  wr[m-1] = b2;
32  wi[m-1] = -sq;
33  numroots+=2;
34  }
35  else {
36  sq = sqrt(disc);
37  wr[m-2] = fabs(b2)+sq;
38  if (b2 < 0.0) wr[m-2] = -wr[m-2];
39  if (wr[m-2] == 0)
40  wr[m-1] = 0;
41  else {
42  wr[m-1] = c/wr[m-2];
43  numroots+=2;
44  }
45  wi[m-2] = 0.0;
46  wi[m-1] = 0.0;
47  }
48  m -= 2;
49  }
50  if (m == 1) {
51  wr[0] = -a[0];
52  wi[0] = 0.0;
53  numroots++;
54  }
55  return numroots;
56 }
57 
60 void deflate(double *a,int n,double *b,double *quad,double *err)
61 {
62  double *c,r,s;
63  int i;
64 
65  c = new double [n+1];
66  r = quad[1];
67  s = quad[0];
68 
69  b[1] = a[1] - r;
70  c[1] = b[1] - r;
71 
72  for (i=2;i<=n;i++){
73  b[i] = a[i] - r * b[i-1] - s * b[i-2];
74  c[i] = b[i] - r * c[i-1] - s * c[i-2];
75  }
76  *err = fabs(b[n])+fabs(b[n-1]);
77  delete [] c;
78 }
79 
84 void find_quad(double *a,int n,double *b,double *quad,double *err, int *iter)
85 {
86  double *c,dn,dr,ds,drn,dsn,eps,r,s;
87  int i;
88 
89  c = new double [n+1];
90  c[0] = 1.0;
91  r = quad[1];
92  s = quad[0];
93  dr = 1.0;
94  ds = 0;
95  eps = 1e-15;
96  *iter = 1;
97 
98  while ((fabs(dr)+fabs(ds)) > eps) {
99  if (*iter > maxiter) break;
100  if (((*iter) % 200) == 0) {
101  eps*=10.0;
102  }
103  b[1] = a[1] - r;
104  c[1] = b[1] - r;
105 
106  for (i=2;i<=n;i++){
107  b[i] = a[i] - r * b[i-1] - s * b[i-2];
108  c[i] = b[i] - r * c[i-1] - s * c[i-2];
109  }
110  dn=c[n-1] * c[n-3] - c[n-2] * c[n-2];
111  drn=b[n] * c[n-3] - b[n-1] * c[n-2];
112  dsn=b[n-1] * c[n-1] - b[n] * c[n-2];
113 
114  if (fabs(dn) < 1e-15) {
115  if (dn < 0.0) dn = -1e-8;
116  else dn = 1e-8;
117  }
118  dr = drn / dn;
119  ds = dsn / dn;
120  r += dr;
121  s += ds;
122  (*iter)++;
123  }
124  quad[0] = s;
125  quad[1] = r;
126  *err = fabs(ds)+fabs(dr);
127  delete [] c;
128 }
129 
131 void diff_poly(double *a,int n,double *b)
132 {
133  double coef;
134  int i;
135 
136  coef = (double)n;
137  b[0] = 1.0;
138  for (i=1;i<n;i++) {
139  b[i] = a[i]*((double)(n-i))/coef;
140  }
141 }
142 
162 void recurse(double *a,int n,double *b,int m,double *quad,
163  double *err,int *iter)
164 {
165  double *c,*x,rs[2],tst;
166 
167  if (fabs(b[m]) < 1e-16) m--; // this bypasses roots at zero
168  if (m == 2) {
169  quad[0] = b[2];
170  quad[1] = b[1];
171  *err = 0;
172  *iter = 0;
173  return;
174  }
175  c = new double [m+1];
176  x = new double [n+1];
177  c[0] = x[0] = 1.0;
178  rs[0] = quad[0];
179  rs[1] = quad[1];
180  *iter = 0;
181  find_quad(b,m,c,rs,err,iter);
182  tst = fabs(rs[0]-quad[0])+fabs(rs[1]-quad[1]);
183  if (*err < 1e-12) {
184  quad[0] = rs[0];
185  quad[1] = rs[1];
186  }
187 // tst will be 'large' if we converge to wrong root
188  if (((*iter > 5) && (tst < 1e-4)) || ((*iter > 20) && (tst < 1e-1))) {
189  diff_poly(b,m,c);
190  recurse(a,n,c,m-1,rs,err,iter);
191  quad[0] = rs[0];
192  quad[1] = rs[1];
193  }
194  delete [] x;
195  delete [] c;
196 }
197 
201 void get_quads(double *a,int n,double *quad,double *x)
202 {
203  double *b,*z,err,tmp;
204  int iter,i,m;
205 
206  if ((tmp = a[0]) != 1.0) {
207  a[0] = 1.0;
208  for (i=1;i<=n;i++) {
209  a[i] /= tmp;
210  }
211  }
212  if (n == 2) {
213  x[0] = a[1];
214  x[1] = a[2];
215  return;
216  }
217  else if (n == 1) {
218  x[0] = a[1];
219  return;
220  }
221  m = n;
222  b = new double [n+1];
223  z = new double [n+1];
224  b[0] = 1.0;
225  for (i=0;i<=n;i++) {
226  z[i] = a[i];
227  x[i] = 0.0;
228  }
229  do {
230  if (n > m) {
231  quad[0] = 3.14159e-1;
232  quad[1] = 2.78127e-1;
233  }
234 loop:
235  find_quad(z,m,b,quad,&err,&iter);
236  if ((err > 1e-7) || (iter > maxiter)) {
237  diff_poly(z,m,b);
238  iter = 0;
239  recurse(z,m,b,m-1,quad,&err,&iter);
240  }
241  deflate(z,m,b,quad,&err);
242  if (err > 0.01) {
244  //cou << "Excessive error!" << endl;
245  //cou << "Enter new trial 'r': ";
246  //cin >> quad[1];
247  //cou << "Enter new trial 's': ";
248  //cin >> quad[0];
249  quad[0] = -2.71828e-1;
250  quad[1] = -3.14159e-1;
251  //cou << "Error finding roots" << endl;
252  //exit(0);
253  }
254  if (err > 1) goto loop;
255  x[m-2] = quad[1];
256  x[m-1] = quad[0];
257  m -= 2;
258  for (i=0;i<=m;i++) {
259  z[i] = b[i];
260  }
261  } while (m > 2);
262  if (m == 2) {
263  x[0] = b[1];
264  x[1] = b[2];
265  }
266  else x[0] = b[1];
267  delete [] z;
268  delete [] b;
269 }
static const double m
mass of electron in grams
void diff_poly(double *a, int n, double *b)
Differentiate polynomial 'a' returning result in 'b'.
Definition: rroots.cpp:131
int roots(double *a, int n, double *wr, double *wi)
Extract individual real or complex roots from list of quadratic factors.
Definition: rroots.cpp:15
#define DBL_EPSILON
some other epsilon and stuff
Definition: rroots.h:26
void find_quad(double *a, int n, double *b, double *quad, double *err, int *iter)
Definition: rroots.cpp:84
void recurse(double *a, int n, double *b, int m, double *quad, double *err, int *iter)
Definition: rroots.cpp:162
#define maxiter
maximum number of iterations
Definition: rroots.h:23
void get_quads(double *a, int n, double *quad, double *x)
Definition: rroots.cpp:201
void deflate(double *a, int n, double *b, double *quad, double *err)
Definition: rroots.cpp:60