[Bug-gsl] [bug #39055] gsl_poly_complex_solve fails on a 15th degree polynomial

[Patrick Alken]

URL:
  <http://savannah.gnu.org/bugs/?39055>

                 Summary: gsl_poly_complex_solve fails on a 15th degree
polynomial
                 Project: GNU Scientific Library
            Submitted by: psa
            Submitted on: Thu 23 May 2013 08:14:21 PM GMT
                Category: Runtime error
                Severity: 3 - Normal
        Operating System: all
                  Status: Confirmed
             Assigned to: None
             Open/Closed: Open
                 Release: 
         Discussion Lock: Any

    _______________________________________________________

Details:

The following code demonstrates a failure of gsl_poly_complex_solve(). I have
commented this out from poly/test.c since it will cause that test to fail
until the problem is fixed.

The QR reduction of the polynomial companion matrix does not converge. The
eigenvalues of the companion matrix can be successfully computed with
gsl_eigen_nonsymm(), but this is not a good long-term solution - the qr
algorithm should be fixed somehow. Until it is I think this bug report is best
left here rather than in poly/test.c where it causes 'make check' to fail.

  {
    /* 15-th order polynomial y = (x + 1) * (x^10 + x^9 + 2 * x^7 + 4 * x^2 +

       4 * x + 8) * (x - 1)^2 * (x - 2)^2

       Problem reported by Munagala Ramanath 
    */

    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.607810742347235857748423812, 0,
                              -1.306698248492076918607989103, 0,
                              -1.000000000000000000000000000, 0,
                              -0.6589385617524094731042279248,
0.8345975728742668637102804993,
                              -0.6589385617524094731042279248,
-0.8345975728742668637102804993,
                              -0.07089111740334127620785158739,
1.135924908758779018494317485,
                              -0.07089111740334127620785158739,
-1.135924908758779018494317485,
                              0.5728474783941085170880376851,
1.198780898828970478006562723,
                              0.5728474783941085170880376851,
-1.198780898828970478006562723,
                              1.114236696181298620402248284,
0.4808398120338998150567661114,
                              1.114236696181298620402248284,
-0.4808398120338998150567661114,
                              2.000000000000000000000000000, 0,
                              2.000000000000000000000000000, 0,
                              1.000000000000000000000000000, 0,
                              1.000000000000000000000000000, 0};


    int i;

    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (16);

    int status = gsl_poly_complex_solve (a, 16, w, z);

    gsl_poly_complex_workspace_free (w);

    gsl_test (status, "gsl_poly_complex_solve, 8th-order polynomial");
    
    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);
    }
  }





    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?39055>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/