Thu 23 May 2013 08:14:21 PM UTC, original submission:
The following code demonstrates a failure of gsl_poly_complex_solve(). I have commented this out from poly/test.c since it will cause that test to fail until the problem is fixed.
The QR reduction of the polynomial companion matrix does not converge. The eigenvalues of the companion matrix can be successfully computed with gsl_eigen_nonsymm(), but this is not a good long-term solution - the qr algorithm should be fixed somehow. Until it is I think this bug report is best left here rather than in poly/test.c where it causes 'make check' to fail.
{
/* 15-th order polynomial y = (x + 1) * (x^10 + x^9 + 2 * x^7 + 4 * x^2 +
4 * x + 8) * (x - 1)^2 * (x - 2)^2
Problem reported by Munagala Ramanath
*/
double a[16] = { 32, -48, -8, 28, -8, 16, -16, 12, -16, 6, 10, -17, 10, 2, -4, 1 };
double z[16*2];
double expected[16*20] = {-1.607810742347235857748423812, 0,
-1.306698248492076918607989103, 0,
-1.000000000000000000000000000, 0,
-0.6589385617524094731042279248, 0.8345975728742668637102804993,
-0.6589385617524094731042279248, -0.8345975728742668637102804993,
-0.07089111740334127620785158739, 1.135924908758779018494317485,
-0.07089111740334127620785158739, -1.135924908758779018494317485,
0.5728474783941085170880376851, 1.198780898828970478006562723,
0.5728474783941085170880376851, -1.198780898828970478006562723,
1.114236696181298620402248284, 0.4808398120338998150567661114,
1.114236696181298620402248284, -0.4808398120338998150567661114,
2.000000000000000000000000000, 0,
2.000000000000000000000000000, 0,
1.000000000000000000000000000, 0,
1.000000000000000000000000000, 0};
int i;
gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (16);
int status = gsl_poly_complex_solve (a, 16, w, z);
gsl_poly_complex_workspace_free (w);
gsl_test (status, "gsl_poly_complex_solve, 8th-order polynomial");
for (i = 0; i<15; i++) {
gsl_test_abs (z[2i], expected[2i], 1e-7, "z%d.real, 15th-order polynomial", i);
gsl_test_abs (z[2i+1], expected[2i+1], 1e-7, "z%d.imag, 15th-order polynomial", i);
}
}
|