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)");

   }