Hi, Great work on the GSL, I love using it in my research.
I'm using the GSL in my research to solve arbitrary polynomials of order up to 4. I've written a little wrapper around the GSL to use gsl_poly_complex_solve_quadratic() for quadratics and gsl_poly_complex_solve_cubic() for cubics, and I use gsl_poly_complex_solve() to solve quartics. I also implemented my own analytical quartic solver based on Ferrari's method (http://mathworld.wolfram.com/QuarticEquation.html) using the GSL complex math operations such as gsl_complex_mul etc. to manipulate the complex numbers involved in Ferrari's solution. When comparing the speed of the GSL gsl_poly_complex_solve() solution of a quartic to my analytical version, I find that the analytic routine runs just a touch faster: *********** timing GSL version... GSL version took 1.078 seconds for 100000 runs; averaging 1.078e-05 seconds per run. *********** Analytic version... Analytic version took 0.797 seconds for 100000 runs; averaging 7.97e-06 seconds per run. (I'm running a 3.2 GHz Pentium IV with 1GB RAM, Cygwin on Windows XP for the timing tests, using the cstdlib rand() function to generate "random" coefficients between 0.0 and 10.0 for the quartics for each run) When comparing the "precision" of the two quartic solvers, I find a very small (machine-level?) difference between the roots that they find, typically 10^-30, using gsl_complex_sub() to calculate the difference between roots. I also find that sometimes the GSL gsl_poly_complex_solve() routine returns a non-zero status int meaning that it didn't find roots or it found unstable/overflow/underflow conditions, I presume. So, my question is, what is my advantage in using the GSL gsl_poly_complex_solve() routine to solve my quartics? I have a feeling the GSL gsl_poly_complex_solve() routine is more robust to large and small coefficients in the polynomial, because I haven't taken too much care in looking out for those in my analytic code. I'm sure the GSL developers looked at Ferrari's method when considering whether or not to implement a gsl_poly_complex_solve_quartic() routine, since an analytic solution is acheivable, and I'm wondering why they didn't implement one? Is there some obvious flaw when implementing Ferrari's method computationally? What are the other advantages of solving quartics using the iterative method in gsl_poly_complex_solve() instead? If you'd like to see any of the code I've used, either for Ferrari's method, or for testing, I'd be happy to share it. Thanks again for the great library. Cheers, Will _______________________________________________ Help-gsl mailing list Help-gsl@gnu.org http://lists.gnu.org/mailman/listinfo/help-gsl