#include #include #include #define polyOrderPlus1 6 /* This is a test program to show failure of gsl_poly_complex_solve method. * Accordining to the documentation, the method should return failure code * but it does call default handler and exit the application instead. This * code relies on the GSL_EFAILED return to fudge the coefficients to find * roots. * * -- Suleyman G */ int main (void) { int i, status, arrayIndex, failCount = 0; const int polyOrder = polyOrderPlus1 - 1; double a[polyOrderPlus1] = { -0.058690412428659, 0.028123386142471, 0.008523354273975, -0.002639648481342, -0.000378023608125, 0.000099773554248 }; double z[2*polyOrder]; gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (polyOrderPlus1); while (status=gsl_poly_complex_solve(a, polyOrderPlus1, w, z) != GSL_SUCCESS) { printf("gsl_poly_complex_solve FAILED for coefficients...\n"); for (i = 0; i < polyOrderPlus1; i++) { printf (" a[%d]=%+.15lf\n", i, a[i]); } if (++failCount > 10) { break; } else { arrayIndex = (failCount-1) % polyOrderPlus1; a[arrayIndex] += 1e-15; } } gsl_poly_complex_workspace_free (w); if (status == GSL_SUCCESS) { printf("gsl_poly_complex_solve FOUND the roots...\n Coefficients:\n"); for (i = 0; i < polyOrderPlus1; i++) { printf (" a[%d] = %+.15lf\n", i, a[i]); } printf(" Roots:\n"); for (i = 0; i < polyOrder; i++) { printf (" z[%d] = %+.8f %+.8fi\n", i, z[2*i], z[2*i+1]); } } return 0; }