/[gsl]/gsl/poly/solve_quartic.c
ViewVC logotype

Diff of /gsl/poly/solve_quartic.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.3 by andybuckley, Thu Oct 16 14:49:37 2003 UTC revision 1.1.2.4 by andybuckley, Thu Oct 23 13:09:59 2003 UTC
# Line 42  int Line 42  int
42  gsl_poly_solve_quartic (double a, double b, double c, double d,  gsl_poly_solve_quartic (double a, double b, double c, double d,
43                          double *x0, double *x1, double *x2, double *x3)                          double *x0, double *x1, double *x2, double *x3)
44  {  {
45    double u[3], v[3], v1, v2;    double u[3];
46    
47    /* remaining complex variables are zarr[4] and w{1,2,3} */    /* remaining complex variables are zarr[4] and w{1,2,3} */
48    gsl_complex zarr[4], w1, w2, w3;    gsl_complex zarr[4], w1, w2, w3;
# Line 183  gsl_poly_solve_quartic (double a, double Line 183  gsl_poly_solve_quartic (double a, double
183     * mt=3 : 2 real roots (disc > 0)     * mt=3 : 2 real roots (disc > 0)
184     */     */
185    double mt;    double mt;
186    if (0 == disc)  
     {  
       u[2] = u[1];  
     }  
187    if (0 >= disc)    if (0 >= disc)
188      {      {
189        mt = 2;        if (0 == disc)
190            {
191              mt = 1;
192              u[2] = u[1]; /* *** ? */
193            }
194          else
195            {
196              mt = 2; /* *** so should we return 0 here? */
197            }
198    
199          double v[3], v1, v2;
200        v[0] = fabs (u[0]);        v[0] = fabs (u[0]);
201        v[1] = fabs (u[1]);        v[1] = fabs (u[1]);
202        v[2] = fabs (u[2]);        v[2] = fabs (u[2]);
203        v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]);        v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]);
204    
205          /* Work out which two roots have the largest moduli */
206        int k1 = 0, k2 = 0;        int k1 = 0, k2 = 0;
207        if (v1 == v[0])        if (v1 == v[0])
208          {          {
# Line 224  gsl_poly_solve_quartic (double a, double Line 232  gsl_poly_solve_quartic (double a, double
232          {          {
233            k2 = 2;            k2 = 2;
234          }          }
235    
236        w1 = gsl_complex_sqrt_real (u[k1]);        w1 = gsl_complex_sqrt_real (u[k1]);
237        w2 = gsl_complex_sqrt_real (u[k2]);        w2 = gsl_complex_sqrt_real (u[k2]);
238      }      }
# Line 238  gsl_poly_solve_quartic (double a, double Line 247  gsl_poly_solve_quartic (double a, double
247    
248    /* Solve the quadratic to obtain the roots to the quartic */    /* Solve the quadratic to obtain the roots to the quartic */
249    double q = qq; /* ! */    double q = qq; /* ! */
250    if (0.0 != gsl_complex_abs (gsl_complex_mul (w1, w2)))    if (0.0 != gsl_complex_abs (gsl_complex_mul (w1, w2))) {
251      {      w3 = gsl_complex_mul_real (gsl_complex_inverse (gsl_complex_mul (w1, w2)), -q / 8.0);
252        w3 =    }
         gsl_complex_mul_real (gsl_complex_inverse (gsl_complex_mul (w1, w2)),  
                               -q / 8.0);  
     }  
253    double h = a / 4.0;    double h = a / 4.0;
254    zarr[0] =    zarr[0] = gsl_complex_add_real (gsl_complex_add (gsl_complex_add (w1, w2), w3), -h);
255      gsl_complex_add_real (gsl_complex_add (gsl_complex_add (w1, w2), w3), -h);    zarr[1] = gsl_complex_add_real (gsl_complex_add (gsl_complex_negative (gsl_complex_add (w1, w2)), w3), -h);
256    zarr[1] =    zarr[2] = gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w2, w1), w3), -h);
257      gsl_complex_add_real (gsl_complex_add    zarr[3] = gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w1, w2), w3), -h);
258                            (gsl_complex_negative (gsl_complex_add (w1, w2)),  
259                             w3), -h);    /* *** the mt==1 situation is never handled...nothing done for assigning to x2 and x3 */
   zarr[2] =  
     gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w2, w1), w3), -h);  
   zarr[3] =  
     gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w1, w2), w3), -h);  
260    
261    /* Arrange the roots into the variables z0, z1, z2, z3 */    /* Arrange the roots into the variables z0, z1, z2, z3 */
262    if (2 == mt)    if (2 == mt)
263      {      {
264        /* No real roots. Return 0; */        /* No real roots. Return 0; */
265          /* *** this can be done much earlier */
266        return 0;        return 0;
267      }      }
268    else if (3 == mt)    else if (3 == mt)

Legend:
Removed from v.1.1.2.3  
changed lines
  Added in v.1.1.2.4

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26