Fixes bug #39055: gsl_poly_complex_solve fails on a 15th degree polynomial
http://savannah.gnu.org/bugs/?39055
Changes are already integrated upstream but not in a released version yet.
--- gsl-1.16/poly/test.c.orig 2015-03-31 10:20:01.997591810 -0700
+++ gsl-1.16/poly/test.c 2015-03-31 10:20:13.280773508 -0700
@@ -25,11 +25,21 @@
#include <gsl/gsl_poly.h>
#include <gsl/gsl_heapsort.h>
+/* sort by Re(z) then by Im(z) */
static int
cmp_cplx(const double *a, const double *b)
{
- double t = (a[0] * a[0] + a[1] * a[1]) - (b[0] * b[0] + b[1] * b[1]);
- return t < 0.0 ? -1 : t > 0.0 ? 1 : 0;
+ double r = a[0] - b[0];
+
+ if (r == 0.0)
+ {
+ double t = a[1] - b[1];
+ return t < 0.0 ? -1 : t > 0.0 ? 1 : 0;
+ }
+ else if (r < 0.0)
+ return -1;
+ else
+ return 1;
}
int
@@ -534,25 +544,26 @@
Problem reported by Munagala Ramanath (bug #39055)
*/
- double a[16] = { 32, -48, -8, 28, -8, 16, -16, 12, -16, 6, 10, -17, 10, 2, -4, 1 };
+ 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.0000000000000000, 0.00000000000000000,
- 1.0000000000000000, 0.00000000000000000,
- -1.0000000000000000, 0.00000000000000000,
- -0.65893856175240950, 0.83459757287426684,
- -0.65893856175240950, -0.83459757287426684,
- -0.070891117403341281, -1.1359249087587791,
- -0.070891117403341281, 1.1359249087587791,
- 1.1142366961812986, -0.48083981203389980,
- 1.1142366961812986, 0.48083981203389980,
- -1.3066982484920768, 0.00000000000000000,
- 0.57284747839410854, 1.1987808988289705,
- 0.57284747839410854, -1.1987808988289705,
- -1.6078107423472359, 0.00000000000000000,
- 2.0000000000000000, 0.00000000000000000,
- 2.0000000000000000, 0.00000000000000000 };
+ double expected[16*2] = {
+ -1.6078107423472359, 0.00000000000000000,
+ -1.3066982484920768, 0.00000000000000000,
+ -1.0000000000000000, 0.00000000000000000,
+ -0.65893856175240950, -0.83459757287426684,
+ -0.65893856175240950, 0.83459757287426684,
+ -0.070891117403341281, -1.1359249087587791,
+ -0.070891117403341281, 1.1359249087587791,
+ 0.57284747839410854, -1.1987808988289705,
+ 0.57284747839410854, 1.1987808988289705,
+ 1.0000000000000000, 0.00000000000000000,
+ 1.0000000000000000, 0.00000000000000000,
+ 1.1142366961812986, -0.48083981203389980,
+ 1.1142366961812986, 0.48083981203389980,
+ 2.0000000000000000, 0.00000000000000000,
+ 2.0000000000000000, 0.00000000000000000 };
int i;
@@ -568,8 +579,8 @@
for (i = 0; i<15; i++)
{
- gsl_test_abs (z[2*i], expected[2*i], 1e-7, "z%d.real, 15th-order polynomial", i);
- gsl_test_abs (z[2*i+1], expected[2*i+1], 1e-7, "z%d.imag, 15th-order polynomial", i);
+ gsl_test_rel (z[2*i], expected[2*i], 1e-7, "z%d.real, 15th-order polynomial", i);
+ gsl_test_rel (z[2*i+1], expected[2*i+1], 1e-7, "z%d.imag, 15th-order polynomial", i);
}
}
@@ -654,12 +665,12 @@
x = -0.5;
gsl_poly_eval_derivs(c, 6, x, dc, 6);
- gsl_test_rel (dc[0], c[0] + c[1]*x + c[2]*x*x + c[3]*x*x*x + c[4]*x*x*x*x + c[5]*x*x*x*x*x , eps, "gsl_poly_eval_dp({+1, -2, +3, -4, +5, -6}, 3.75)");
- gsl_test_rel (dc[1], c[1] + 2.0*c[2]*x + 3.0*c[3]*x*x + 4.0*c[4]*x*x*x + 5.0*c[5]*x*x*x*x , eps, "gsl_poly_eval_dp({+1, -2, +3, -4, +5, -6} deriv 1, -12.375)");
- gsl_test_rel (dc[2], 2.0*c[2] + 3.0*2.0*c[3]*x + 4.0*3.0*c[4]*x*x + 5.0*4.0*c[5]*x*x*x , eps, "gsl_poly_eval_dp({+1, -2, +3, -4, +5, -6} deriv 2, +48.0)");
- gsl_test_rel (dc[3], 3.0*2.0*c[3] + 4.0*3.0*2.0*c[4]*x + 5.0*4.0*3.0*c[5]*x*x , eps,"gsl_poly_eval_dp({+1, -2, +3, -4, +5, -6} deriv 3, -174.0)");
- gsl_test_rel (dc[4], 4.0*3.0*2.0*c[4] + 5.0*4.0*3.0*2.0*c[5]*x, eps, "gsl_poly_eval_dp({+1, -2, +3, -4, +5, -6} deriv 4, +480.0)");
- gsl_test_rel (dc[5], 5.0*4.0*3.0*2.0*c[5] , eps, "gsl_poly_eval_dp({+1, -2, +3, -4, +5, -6} deriv 5, -720.0)");
+ gsl_test_rel (dc[0], c[0] + c[1]*x + c[2]*x*x + c[3]*x*x*x + c[4]*x*x*x*x + c[5]*x*x*x*x*x , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6}, 3.75)");
+ gsl_test_rel (dc[1], c[1] + 2.0*c[2]*x + 3.0*c[3]*x*x + 4.0*c[4]*x*x*x + 5.0*c[5]*x*x*x*x , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 1, -12.375)");
+ gsl_test_rel (dc[2], 2.0*c[2] + 3.0*2.0*c[3]*x + 4.0*3.0*c[4]*x*x + 5.0*4.0*c[5]*x*x*x , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 2, +48.0)");
+ gsl_test_rel (dc[3], 3.0*2.0*c[3] + 4.0*3.0*2.0*c[4]*x + 5.0*4.0*3.0*c[5]*x*x , eps,"gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 3, -174.0)");
+ gsl_test_rel (dc[4], 4.0*3.0*2.0*c[4] + 5.0*4.0*3.0*2.0*c[5]*x, eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 4, +480.0)");
+ gsl_test_rel (dc[5], 5.0*4.0*3.0*2.0*c[5] , eps, "gsl_poly_eval_derivs({+1, -2, +3, -4, +5, -6} deriv 5, -720.0)");
}