From 89f2d5667fd1ace51ab3e05ad614adbbcff8dab8 Mon Sep 17 00:00:00 2001
From: Manuel Yguel <manuel.yguel@gmail.com>
Date: Thu, 25 Mar 2010 03:25:47 +0100
Subject: [PATCH] Add the possibility to use the polynomial solver of the gsl.

---
 test/gsl_helper.h | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/test/gsl_helper.h b/test/gsl_helper.h
index f8f7a5d79..eefffe792 100644
--- a/test/gsl_helper.h
+++ b/test/gsl_helper.h
@@ -33,6 +33,7 @@
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_complex.h>
 #include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_poly.h>
 
 namespace Eigen {
 
@@ -69,6 +70,27 @@ template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex> struct
     gsl_eigen_gensymmv_free(w);
     free(a);
   }
+
+  template<class EIGEN_VECTOR, class EIGEN_ROOTS>
+  static void eigen_poly_solve(const EIGEN_VECTOR& poly, EIGEN_ROOTS& roots )
+  {
+    const int deg = poly.size()-1;
+    double *z = new double[2*deg];
+    double *a = new double[poly.size()];
+    for( int i=0; i<poly.size(); ++i ){ a[i] = poly[i]; }
+
+    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (poly.size());
+    gsl_poly_complex_solve(a, poly.size(), w, z);
+    gsl_poly_complex_workspace_free (w);
+    for( int i=0; i<deg; ++i )
+    {
+      roots[i].real() = z[2*i];
+      roots[i].imag() = z[2*i+1];
+    }
+
+    delete[] a;
+    delete[] z;
+  }
 };
 
 template<typename Scalar> struct GslTraits<Scalar,true>