15 int roots(
double *a,
int n,
double *wr,
double *wi)
26 if (fabs(disc)/(fabs(b2*b2)+fabs(c)) <=
DBL_EPSILON) disc = 0.0;
37 wr[m-2] = fabs(b2)+sq;
38 if (b2 < 0.0) wr[m-2] = -wr[m-2];
60 void deflate(
double *a,
int n,
double *b,
double *quad,
double *err)
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];
76 *err = fabs(b[n])+fabs(b[n-1]);
84 void find_quad(
double *a,
int n,
double *b,
double *quad,
double *err,
int *iter)
86 double *c,dn,dr,ds,drn,dsn,eps,r,s;
98 while ((fabs(dr)+fabs(ds)) > eps) {
100 if (((*iter) % 200) == 0) {
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];
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];
114 if (fabs(dn) < 1e-15) {
115 if (dn < 0.0) dn = -1e-8;
126 *err = fabs(ds)+fabs(dr);
139 b[i] = a[i]*((double)(n-i))/coef;
162 void recurse(
double *a,
int n,
double *b,
int m,
double *quad,
163 double *err,
int *iter)
165 double *c,*x,rs[2],tst;
167 if (fabs(b[m]) < 1e-16) m--;
175 c =
new double [m+1];
176 x =
new double [n+1];
182 tst = fabs(rs[0]-quad[0])+fabs(rs[1]-quad[1]);
188 if (((*iter > 5) && (tst < 1e-4)) || ((*iter > 20) && (tst < 1e-1))) {
190 recurse(a,n,c,m-1,rs,err,iter);
203 double *b,*z,err,tmp;
206 if ((tmp = a[0]) != 1.0) {
222 b =
new double [n+1];
223 z =
new double [n+1];
231 quad[0] = 3.14159e-1;
232 quad[1] = 2.78127e-1;
236 if ((err > 1e-7) || (iter >
maxiter)) {
239 recurse(z,m,b,m-1,quad,&err,&iter);
249 quad[0] = -2.71828e-1;
250 quad[1] = -3.14159e-1;
254 if (err > 1)
goto loop;
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'.
int roots(double *a, int n, double *wr, double *wi)
Extract individual real or complex roots from list of quadratic factors.
#define DBL_EPSILON
some other epsilon and stuff
void find_quad(double *a, int n, double *b, double *quad, double *err, int *iter)
void recurse(double *a, int n, double *b, int m, double *quad, double *err, int *iter)
#define maxiter
maximum number of iterations
void get_quads(double *a, int n, double *quad, double *x)
void deflate(double *a, int n, double *b, double *quad, double *err)