/[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.1 by andybuckley, Thu Oct 9 14:27:30 2003 UTC revision 1.1.2.2 by andybuckley, Thu Oct 16 12:52:49 2003 UTC
# Line 50  gsl_poly_solve_quartic (double a, double Line 50  gsl_poly_solve_quartic (double a, double
50    double aa, pp, qq, rr, rc, sc, tc, q, h;    double aa, pp, qq, rr, rc, sc, tc, q, h;
51    int k1 = 0, k2 = 0, mt;    int k1 = 0, k2 = 0, mt;
52    
   gsl_complex tmp[3];  
   
53    GSL_SET_COMPLEX (&i, 0.0, 1.0);    GSL_SET_COMPLEX (&i, 0.0, 1.0);
54    GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);    GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);
55    GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);    GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);
# Line 63  gsl_poly_solve_quartic (double a, double Line 61  gsl_poly_solve_quartic (double a, double
61    
62    /* Deal easily with the cases where the quartic is degenerate. The    /* Deal easily with the cases where the quartic is degenerate. The
63     * ordering of solutions is done explicitly. */     * ordering of solutions is done explicitly. */
64    if (b == 0 && c == 0)    if (0 == b && 0 == c)
65      {      {
66        if (d == 0)        if (0 == d)
67          {          {
68            if (a > 0)            if (a > 0)
69              {              {
# Line 83  gsl_poly_solve_quartic (double a, double Line 81  gsl_poly_solve_quartic (double a, double
81              }              }
82            return 4;            return 4;
83          }          }
84        else if (a == 0)        else if (0 == a)
85          {          {
86            if (d > 0)            if (d > 0)
87              {              {
# Line 131  gsl_poly_solve_quartic (double a, double Line 129  gsl_poly_solve_quartic (double a, double
129    
130      disc = (CR2 - CQ3) / 2125764.0;      disc = (CR2 - CQ3) / 2125764.0;
131    
132      if (R == 0 && Q == 0)      if (0 == R && 0 == Q)
133        {        {
134          u[0] = -rc / 3;          u[0] = -rc / 3;
135          u[1] = -rc / 3;          u[1] = -rc / 3;
# Line 158  gsl_poly_solve_quartic (double a, double Line 156  gsl_poly_solve_quartic (double a, double
156          double sqrtQ = sqrt (Q);          double sqrtQ = sqrt (Q);
157          double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;          double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
158          double theta = acos (R / sqrtQ3);          double theta = acos (R / sqrtQ3);
159          if (R / sqrtQ3 >= 1.0)          if (R / sqrtQ3 >= 1.0) theta = 0.0;
           theta = 0.0;  
160          double norm = -2 * sqrtQ;          double norm = -2 * sqrtQ;
161    
162          u[0] = norm * cos (theta / 3) - rc / 3;          u[0] = norm * cos (theta / 3) - rc / 3;
# Line 169  gsl_poly_solve_quartic (double a, double Line 166  gsl_poly_solve_quartic (double a, double
166      else      else
167        {        {
168          double sgnR = (R >= 0 ? 1 : -1);          double sgnR = (R >= 0 ? 1 : -1);
169          double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0 / 3.0);          double modR = fabs (R);
170            double sqrt_disc = sqrt (R2 - Q3);
171            double A = -sgnR * pow (modR + sqrt_disc, 1.0 / 3.0);
172          double B = Q / A;          double B = Q / A;
173            double mod_diffAB = fabs (A - B);
174    
175          u[0] = A + B - rc / 3;          u[0] = A + B - rc / 3;
176          u[1] = -0.5 * (A + B) - rc / 3;          u[1] = -0.5 * (A + B) - rc / 3;
177          u[2] = -(sqrt (3.0) / 2.0) * fabs (A - B);          u[2] = -(sqrt (3.0) / 2.0) * mod_diffAB;
178        }        }
179    }    }
180    /* End of solution to resolvent cubic */    /* End of solution to resolvent cubic */
# Line 186  gsl_poly_solve_quartic (double a, double Line 186  gsl_poly_solve_quartic (double a, double
186     * mt=2 : 0 real roots (disc < 0)     * mt=2 : 0 real roots (disc < 0)
187     * mt=3 : 2 real roots (disc > 0)     * mt=3 : 2 real roots (disc > 0)
188     */     */
189    if (disc == 0)    if (0 == disc)
190      {      {
191        u[2] = u[1];        u[2] = u[1];
192      }      }
193    if (disc <= 0)    if (0 >= disc)
194      {      {
195        mt = 2;        mt = 2;
196        v[0] = fabs (u[0]);        v[0] = fabs (u[0]);
# Line 237  gsl_poly_solve_quartic (double a, double Line 237  gsl_poly_solve_quartic (double a, double
237        w2 = gsl_complex_sqrt (w2);        w2 = gsl_complex_sqrt (w2);
238      }      }
239    
240    /* Solve the quadratic in order to obtain the roots    /* Solve the quadratic to obtain the roots to the quartic */
    * to the quartic */  
241    q = qq;    q = qq;
242    if (gsl_complex_abs (gsl_complex_mul (w1, w2)) != 0.0)    if (0.0 != gsl_complex_abs (gsl_complex_mul (w1, w2)))
243      {      {
244        w3 =        w3 =
245          gsl_complex_mul_real (gsl_complex_inverse (gsl_complex_mul (w1, w2)),          gsl_complex_mul_real (gsl_complex_inverse (gsl_complex_mul (w1, w2)),
# Line 259  gsl_poly_solve_quartic (double a, double Line 258  gsl_poly_solve_quartic (double a, double
258      gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w1, w2), w3), -h);      gsl_complex_add_real (gsl_complex_sub (gsl_complex_sub (w1, w2), w3), -h);
259    
260    /* Arrange the roots into the variables z0, z1, z2, z3 */    /* Arrange the roots into the variables z0, z1, z2, z3 */
261    if (mt == 2)    if (2 == mt)
262      {      {
263        /* No real roots. Return 0; */        /* No real roots. Return 0; */
264        return 0;        return 0;
265      }      }
266    else if (mt == 3)    else if (3 == mt)
267      {      {
268        *x0 = GSL_REAL (zarr[0]);        *x0 = GSL_REAL (zarr[0]);
269        *x1 = GSL_REAL (zarr[1]);        *x1 = GSL_REAL (zarr[1]);
270      }      }
271    
272    /*    /* Sort the roots as usual */
273     * Sort the roots as usual: main sorting by ascending real part, secondary    if (1 == mt)
    * sorting by ascending imaginary part  
    */  
   
   if (mt == 1)  
274      {      {
275        /* Roots are all real, sort them by the real part */        /* Roots are all real, sort them by the real part */
276        if (*x0 > *x1)        if (*x0 > *x1)

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