47 #define MAT(m,i,j,n) ((m)[(i)*(n) + (j)]) 50 #define FMAT(m,i,j,n) ((m)[((i)-1)*(n) + ((j)-1)]) 53 #define GSL_DBL_EPSILON 2.2204460492503131e-16 56 #define RADIX2 (RADIX*RADIX) 59 #define GSL_FAILURE -1 63 #define GSL_SET_COMPLEX_PACKED(zp,n,x,y) do {*((zp)+2*(n))=(x); *((zp)+(2*(n)+1))=(y);} while(0) 77 gsl_poly_complex_workspace_alloc (
size_t n);
83 set_companion_matrix (
const double *a,
size_t n,
double *m);
89 balance_companion_matrix (
double *m,
size_t n);
92 gsl_poly_complex_solve (
const double *a,
size_t n,
108 assert (roots.
nrows() == a - 1);
109 assert (roots.
ncols() == 2);
110 assert (coeffs[a - 1] != 0);
116 c = (
double *)malloc (a *
sizeof (
double));
117 for (
Index i = 0; i < a; i++)
118 c[i] = (
double)coeffs[i];
120 s = (
double *)malloc ((a-1) * 2 *
sizeof (double));
125 int status = gsl_poly_complex_solve (c, a, w, s);
130 for (
Index i = 0; i < a - 1; i++)
132 roots (i, 0) = (
Numeric)s[i * 2];
133 roots (i, 1) = (
Numeric)s[i * 2 + 1];
141 gsl_poly_complex_workspace_free (w);
148 gsl_poly_complex_workspace_alloc (
size_t n)
156 cerr <<
"matrix size n must be positive integer" << endl;
166 cerr <<
"failed to allocate space for struct" << endl;
175 w->
matrix = (
double *) malloc (nc * nc *
sizeof(
double));
181 cerr <<
"failed to allocate space for workspace matrix" << endl;
199 balance_companion_matrix (
double *m,
size_t nc)
201 int not_converged = 1;
206 while (not_converged)
213 for (i = 0; i < nc; i++)
219 col_norm = fabs (
MAT (m, i + 1, i, nc));
225 for (j = 0; j < nc - 1; j++)
227 col_norm += fabs (
MAT (m, j, nc - 1, nc));
235 row_norm = fabs (
MAT (m, 0, nc - 1, nc));
237 else if (i == nc - 1)
239 row_norm = fabs (
MAT (m, i, i - 1, nc));
243 row_norm = (fabs (
MAT (m, i, i - 1, nc))
244 + fabs (
MAT (m, i, nc - 1, nc)));
247 if (col_norm == 0 || row_norm == 0)
252 g = row_norm /
RADIX;
254 s = col_norm + row_norm;
262 g = row_norm *
RADIX;
270 if ((row_norm + col_norm) < 0.95 * s * f)
278 MAT (m, 0, nc - 1, nc) *= g;
282 MAT (m, i, i - 1, nc) *= g;
283 MAT (m, i, nc - 1, nc) *= g;
288 for (j = 0; j < nc; j++)
290 MAT (m, j, i, nc) *= f;
295 MAT (m, i + 1, i, nc) *= f;
308 size_t iterations, e, i, j, k, m;
310 double w, x, y, s, z;
312 double p = 0,
q = 0, r = 0;
332 for (e = n; e >= 2; e--)
334 double a1 = fabs (
FMAT (h, e, e - 1, nc));
335 double a2 = fabs (
FMAT (h, e - 1, e - 1, nc));
336 double a3 = fabs (
FMAT (h, e, e, nc));
342 x =
FMAT (h, n, n, nc);
352 y =
FMAT (h, n - 1, n - 1, nc);
353 w =
FMAT (h, n - 1, n, nc) *
FMAT (h, n, n - 1, nc);
385 if (iterations == 60)
392 if (iterations % 10 == 0 && iterations > 0)
397 for (i = 1; i <= n; i++)
399 FMAT (h, i, i, nc) -= x;
402 s = fabs (
FMAT (h, n, n - 1, nc)) + fabs (
FMAT (h, n - 1, n - 2, nc));
410 for (m = n - 2; m >= e; m--)
414 z =
FMAT (h, m, m, nc);
417 p =
FMAT (h, m, m + 1, nc) + (r * s -
w) /
FMAT (h, m + 1, m, nc);
418 q =
FMAT (h, m + 1, m + 1, nc) - z - r - s;
419 r =
FMAT (h, m + 2, m + 1, nc);
420 s = fabs (p) + fabs (
q) + fabs (r);
428 a1 = fabs (
FMAT (h, m, m - 1, nc));
429 a2 = fabs (
FMAT (h, m - 1, m - 1, nc));
430 a3 = fabs (
FMAT (h, m + 1, m + 1, nc));
436 for (i = m + 2; i <= n; i++)
438 FMAT (h, i, i - 2, nc) = 0;
441 for (i = m + 3; i <= n; i++)
443 FMAT (h, i, i - 3, nc) = 0;
448 for (k = m; k <= n - 1; k++)
450 notlast = (k != n - 1);
454 p =
FMAT (h, k, k - 1, nc);
455 q =
FMAT (h, k + 1, k - 1, nc);
456 r = notlast ?
FMAT (h, k + 2, k - 1, nc) : 0.0;
458 x = fabs (p) + fabs (
q) + fabs (r);
468 s = sqrt (p * p +
q *
q + r * r);
475 FMAT (h, k, k - 1, nc) = -s * x;
479 FMAT (h, k, k - 1, nc) *= -1;
491 for (j = k; j <= n; j++)
493 p =
FMAT (h, k, j, nc) +
q *
FMAT (h, k + 1, j, nc);
497 p += r *
FMAT (h, k + 2, j, nc);
498 FMAT (h, k + 2, j, nc) -= p * z;
501 FMAT (h, k + 1, j, nc) -= p * y;
502 FMAT (h, k, j, nc) -= p * x;
505 j = (k + 3 < n) ? (k + 3) : n;
509 for (i = e; i <= j; i++)
511 p = x *
FMAT (h, i, k, nc) + y *
FMAT (h, i, k + 1, nc);
515 p += z *
FMAT (h, i, k + 2, nc);
516 FMAT (h, i, k + 2, nc) -= p * r;
518 FMAT (h, i, k + 1, nc) -= p *
q;
519 FMAT (h, i, k, nc) -= p;
528 set_companion_matrix (
const double *a,
size_t nc,
double *m)
532 for (i = 0; i < nc; i++)
533 for (j = 0; j < nc; j++)
534 MAT (m, i, j, nc) = 0.0;
536 for (i = 1; i < nc; i++)
537 MAT (m, i, i - 1, nc) = 1.0;
539 for (i = 0; i < nc; i++)
540 MAT (m, i, nc - 1, nc) = -a[i] / a[nc];
545 gsl_poly_complex_solve (
const double *a,
size_t n,
554 cerr <<
"number of terms must be a positive integer" << endl;
561 cerr <<
"cannot solve for only one term" << endl;
568 cerr <<
"leading term of polynomial must be non-zero" << endl;
575 cerr <<
"size of workspace does not match polynomial" << endl;
583 set_companion_matrix (a, n - 1, m);
585 balance_companion_matrix (m, n - 1);
587 status = qr_companion (m, n - 1, z);
INDEX Index
The type to use for all integer numbers and indices.
Numeric * mdata
Pointer to the plain C array that holds the data.
Contains the code to determine roots of polynomials.
#define GSL_SET_COMPLEX_PACKED(zp, n, x, y)
cmplx FADDEEVA() w(cmplx z, double relerr)
Index nelem() const
Returns the number of elements.
Index ncols() const
Returns the number of columns.
double * gsl_complex_packed_ptr
NUMERIC Numeric
The type to use for all floating point numbers.
Numeric * mdata
Pointer to the plain C array that holds the data.
int poly_root_solve(Matrix &roots, Vector &coeffs)
Index nrows() const
Returns the number of rows.