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

Diff of /gsl/poly/zsolve_quartic.c

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

revision 1.1.2.1 by andybuckley, Thu Oct 9 14:27:30 2003 UTC revision 1.1.2.2 by andybuckley, Thu Oct 16 12:41:39 2003 UTC
# Line 46  gsl_poly_complex_solve_quartic (double a Line 46  gsl_poly_complex_solve_quartic (double a
46    double aa, pp, qq, rr, rc, sc, tc, q, h;    double aa, pp, qq, rr, rc, sc, tc, q, h;
47    int k1 = 0, k2 = 0, mt;    int k1 = 0, k2 = 0, mt;
48    
   gsl_complex tmp[3];  
   
49    GSL_SET_COMPLEX (&i, 0.0, 1.0);    GSL_SET_COMPLEX (&i, 0.0, 1.0);
50    GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);    GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);
51    GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);    GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);
# Line 59  gsl_poly_complex_solve_quartic (double a Line 57  gsl_poly_complex_solve_quartic (double a
57    
58    /* Deal easily with the cases where the quartic is degenerate. The    /* Deal easily with the cases where the quartic is degenerate. The
59     * ordering of solutions is done explicitly. */     * ordering of solutions is done explicitly. */
60    if (b == 0 && c == 0)    if (0 == b && 0 == c)
61      {      {
62        if (d == 0)        if (0 == d)
63          {          {
64            if (a > 0)            if (a > 0)
65              {              {
# Line 79  gsl_poly_complex_solve_quartic (double a Line 77  gsl_poly_complex_solve_quartic (double a
77              }              }
78            return 4;            return 4;
79          }          }
80        else if (a == 0)        else if (0 == a)
81          {          {
82            if (d > 0)            if (d > 0)
83              {              {
84                *z3 = gsl_complex_sqrt (gsl_complex_mul_real (i, sqrt (d)));                double sqrt_d = sqrt (d);
85                *z2 = gsl_complex_mul (gsl_complex_conjugate (i), *z3);                gsl_complex i_sqrt_d = gsl_complex_mul_real (i, sqrt_d);
86                  gsl_complex minus_i = gsl_complex_conjugate (i);
87                  *z3 = gsl_complex_sqrt (i_sqrt_d);
88                  *z2 = gsl_complex_mul (minus_i, *z3);
89                *z1 = gsl_complex_negative (*z2);                *z1 = gsl_complex_negative (*z2);
90                *z0 = gsl_complex_negative (*z3);                *z0 = gsl_complex_negative (*z3);
91              }              }
92            else            else
93              {              {
94                *z3 = gsl_complex_sqrt_real (sqrt (-d));                double sqrt_abs_d = sqrt (-d);
95                  *z3 = gsl_complex_sqrt_real (sqrt_abs_d);
96                *z2 = gsl_complex_mul (i, *z3);                *z2 = gsl_complex_mul (i, *z3);
97                *z1 = gsl_complex_negative (*z2);                *z1 = gsl_complex_negative (*z2);
98                *z0 = gsl_complex_negative (*z3);                *z0 = gsl_complex_negative (*z3);
# Line 99  gsl_poly_complex_solve_quartic (double a Line 101  gsl_poly_complex_solve_quartic (double a
101          }          }
102      }      }
103    
104    if (c == 0.0 && d == 0.0)    if (0.0 == c && 0.0 == d)
105      {      {
106        disc = a * a - 4.0 * b;        disc = (a * a - 4.0 * b);
107        if (disc < 0.0)        if (disc < 0.0)
108          {          {
109            mt = 3;            mt = 3;
# Line 149  gsl_poly_complex_solve_quartic (double a Line 151  gsl_poly_complex_solve_quartic (double a
151    
152          disc = (CR2 - CQ3) / 2125764.0;          disc = (CR2 - CQ3) / 2125764.0;
153    
154          if (R == 0 && Q == 0)          if (0 == R && 0 == Q)
155            {            {
156              u[0] = -rc / 3;              u[0] = -rc / 3;
157              u[1] = -rc / 3;              u[1] = -rc / 3;
# Line 187  gsl_poly_complex_solve_quartic (double a Line 189  gsl_poly_complex_solve_quartic (double a
189          else          else
190            {            {
191              double sgnR = (R >= 0 ? 1 : -1);              double sgnR = (R >= 0 ? 1 : -1);
192              double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0 / 3.0);              double modR = fabs (R);
193                double sqrt_disc = sqrt (R2 - Q3);
194                double A = -sgnR * pow (modR + sqrt_disc, 1.0 / 3.0);
195              double B = Q / A;              double B = Q / A;
196    
197                double mod_diffAB = fabs (A - B);
198              u[0] = A + B - rc / 3;              u[0] = A + B - rc / 3;
199              u[1] = -0.5 * (A + B) - rc / 3;              u[1] = -0.5 * (A + B) - rc / 3;
200              u[2] = -(sqrt (3.0) / 2.0) * fabs (A - B);              u[2] = -(sqrt (3.0) / 2.0) * mod_diffAB;
201            }            }
202        }        }
203        /* End of solution to resolvent cubic */        /* End of solution to resolvent cubic */
# Line 204  gsl_poly_complex_solve_quartic (double a Line 209  gsl_poly_complex_solve_quartic (double a
209         * mt=2 : 0 real roots         * mt=2 : 0 real roots
210         * mt=3 : 2 real roots         * mt=3 : 2 real roots
211         */         */
212        if (disc == 0)        if (0 == disc)
213          {          {
214            u[2] = u[1];            u[2] = u[1];
215          }          }
216        if (disc <= 0)        if (0 => disc)
217          {          {
218            mt = 2;            mt = 2;
219            v[0] = fabs (u[0]);            v[0] = fabs (u[0]);
220            v[1] = fabs (u[1]);            v[1] = fabs (u[1]);
221            v[2] = fabs (u[2]);            v[2] = fabs (u[2]);
222    
223            v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]);            v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]);
224            if (v1 == v[0])            if (v1 == v[0])
225              {              {
# Line 257  gsl_poly_complex_solve_quartic (double a Line 263  gsl_poly_complex_solve_quartic (double a
263        /* Solve the quadratic in order to obtain the roots        /* Solve the quadratic in order to obtain the roots
264         * to the quartic */         * to the quartic */
265        q = qq;        q = qq;
266        if (gsl_complex_abs (gsl_complex_mul (w1, w2)) != 0.0)        gsl_complex prod_w = gsl_complex_mul (w1, w2);
267          gsl_complex mod_prod_w = gsl_complex_abs (prod_w);
268          if (0.0 != mod_prod_w)
269          {          {
270            w3 =            gsl_complex inv_prod_w = gsl_complex_inverse (prod_w);
271              gsl_complex_mul_real (gsl_complex_inverse            w3 = gsl_complex_mul_real (inv_prod_w, -q / 8.0);
                                   (gsl_complex_mul (w1, w2)), -q / 8.0);  
272          }          }
273    
274        h = r4 * a;        h = r4 * a;
275        zarr[0] =        gsl_complex sum_w12 = gsl_complex_add (w1, w2);
276          gsl_complex_add_real (gsl_complex_add (gsl_complex_add (w1, w2), w3),        gsl_complex neg_sum_w12 = gsl_complex_negative (sum_w12);
277                                -h);        gsl_complex sum_w123 = gsl_complex_add (sum_w12, w3);
278        zarr[1] =        gsl_complex neg_sum_w123 = gsl_complex_add (neg_sum_w12, w3);
279          gsl_complex_add_real (gsl_complex_add  
280                                (gsl_complex_negative        gsl_complex diff_w12 = gsl_complex_sub (w2, w1);
281                                 (gsl_complex_add (w1, w2)), w3), -h);        gsl_complex neg_diff_w12 = gsl_complex_negative (diff_w12);
282        zarr[2] =        gsl_complex diff_w123 = gsl_complex_sub (diff_w12, w3);
283          gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w2, w1), w3),        gsl_complex neg_diff_w123 = gsl_complex_sub (neg_diff_w12, w3);
284                                -h);  
285        zarr[3] =        zarr[0] = gsl_complex_add_real (sum_w123, -h);
286          gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w1, w2), w3),        zarr[1] = gsl_complex_add_real (neg_sum_w123, -h);
287                                -h);        zarr[2] = gsl_complex_add_real (diff_w123, -h);
288          zarr[3] = gsl_complex_add_real (neg_diff_w123, -h);
289    
290        /* Arrange the roots into the variables z0, z1, z2, z3 */        /* Arrange the roots into the variables z0, z1, z2, z3 */
291        if (mt == 2)        if (2 == mt)
292          {          {
293            if (u[k1] >= 0 && u[k2] >= 0)            if (u[k1] >= 0 && u[k2] >= 0)
294              {              {
# Line 311  gsl_poly_complex_solve_quartic (double a Line 320  gsl_poly_complex_solve_quartic (double a
320                *z3 = zarr[2];                *z3 = zarr[2];
321              }              }
322          }          }
323        else if (mt == 3)        else if (3 == mt)
324          {          {
325            GSL_SET_COMPLEX (z0, GSL_REAL (zarr[0]), 0.0);            GSL_SET_COMPLEX (z0, GSL_REAL (zarr[0]), 0.0);
326            GSL_SET_COMPLEX (z1, GSL_REAL (zarr[1]), 0.0);            GSL_SET_COMPLEX (z1, GSL_REAL (zarr[1]), 0.0);
# Line 325  gsl_poly_complex_solve_quartic (double a Line 334  gsl_poly_complex_solve_quartic (double a
334     * sorting by ascending imaginary part     * sorting by ascending imaginary part
335     */     */
336    
337    if (mt == 1)    if (1 == mt)
338      {      {
339        /* Roots are all real, sort them by the real part */        /* Roots are all real, sort them by the real part */
340        if (GSL_REAL (*z0) > GSL_REAL (*z1))        if (GSL_REAL (*z0) > GSL_REAL (*z1)) SWAP (*z0, *z1);
341          SWAP (*z0, *z1);        if (GSL_REAL (*z0) > GSL_REAL (*z2)) SWAP (*z0, *z2);
342        if (GSL_REAL (*z0) > GSL_REAL (*z2))        if (GSL_REAL (*z0) > GSL_REAL (*z3)) SWAP (*z0, *z3);
         SWAP (*z0, *z2);  
       if (GSL_REAL (*z0) > GSL_REAL (*z3))  
         SWAP (*z0, *z3);  
343    
344        if (GSL_REAL (*z1) > GSL_REAL (*z2))        if (GSL_REAL (*z1) > GSL_REAL (*z2)) SWAP (*z1, *z2);
         SWAP (*z1, *z2);  
345        if (GSL_REAL (*z2) > GSL_REAL (*z3))        if (GSL_REAL (*z2) > GSL_REAL (*z3))
346          {          {
347            SWAP (*z2, *z3);            SWAP (*z2, *z3);
348            if (GSL_REAL (*z1) > GSL_REAL (*z2))            if (GSL_REAL (*z1) > GSL_REAL (*z2)) SWAP (*z1, *z2);
             SWAP (*z1, *z2);  
349          }          }
350      }      }
351    else if (mt == 2)    else if (2 == mt)
352      {      {
353        /* Roots are all complex. z0 and z1 are conjugates        /* Roots are all complex. z0 and z1 are conjugates
354         * and z2 and z3 are conjugates. Sort the real parts first */         * and z2 and z3 are conjugates. Sort the real parts first */
# Line 354  gsl_poly_complex_solve_quartic (double a Line 358  gsl_poly_complex_solve_quartic (double a
358            SWAP (*z1, *z3);            SWAP (*z1, *z3);
359          }          }
360        /* Then sort by the imaginary parts */        /* Then sort by the imaginary parts */
361        if (GSL_IMAG (*z0) > GSL_IMAG (*z1))        if (GSL_IMAG (*z0) > GSL_IMAG (*z1)) SWAP (*z0, *z1);
362          SWAP (*z0, *z1);        if (GSL_IMAG (*z2) > GSL_IMAG (*z3)) SWAP (*z2, *z3);
       if (GSL_IMAG (*z2) > GSL_IMAG (*z3))  
         SWAP (*z2, *z3);  
363      }      }
364    else    else
365      {      {
366        /* 2 real roots. z2 and z3 are conjugates. */        /* 2 real roots. z2 and z3 are conjugates. */
367    
368        /* Swap complex roots */        /* Swap complex roots */
369        if (GSL_IMAG (*z2) > GSL_IMAG (*z3))        if (GSL_IMAG (*z2) > GSL_IMAG (*z3)) SWAP (*z2, *z3);
         SWAP (*z2, *z3);  
370    
371        /* Sort real parts */        /* Sort real parts */
372        if (GSL_REAL (*z0) > GSL_REAL (*z1))        if (GSL_REAL (*z0) > GSL_REAL (*z1)) SWAP (*z0, *z1);
         SWAP (*z0, *z1);  
373        if (GSL_REAL (*z1) > GSL_REAL (*z2))        if (GSL_REAL (*z1) > GSL_REAL (*z2))
374          {          {
375            if (GSL_REAL (*z0) > GSL_REAL (*z2))            if (GSL_REAL (*z0) > GSL_REAL (*z2))
# Line 383  gsl_poly_complex_solve_quartic (double a Line 383  gsl_poly_complex_solve_quartic (double a
383                SWAP (*z2, *z3);                SWAP (*z2, *z3);
384              }              }
385          }          }
   
386      }      }
387    
388    return 4;    return 4;

Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.1.2.2

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