/[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.2 by andybuckley, Thu Oct 16 12:52:49 2003 UTC revision 1.1.2.3 by andybuckley, Thu Oct 16 14:49:37 2003 UTC
# Line 31  Line 31 
31    
32  /* Work needed to remove the complex numbers from this  /* Work needed to remove the complex numbers from this
33   * algorithm: when done, compilation will work without these   * algorithm: when done, compilation will work without these
34   * includes   * #includes
35   */   */
36  #include <gsl/gsl_complex.h>  #include <gsl/gsl_complex.h>
37  #include <gsl/gsl_complex_math.h>  #include <gsl/gsl_complex_math.h>
# 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    gsl_complex i, zarr[4], w1, w2, w3;    double u[3], v[3], v1, v2;
   double r4 = 1.0 / 4.0;  
   double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0;  
   double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0;  
   double u[3], v[3], v1, v2, disc;  
   double aa, pp, qq, rr, rc, sc, tc, q, h;  
   int k1 = 0, k2 = 0, mt;  
46    
47    GSL_SET_COMPLEX (&i, 0.0, 1.0);    /* remaining complex variables are zarr[4] and w{1,2,3} */
48      gsl_complex zarr[4], w1, w2, w3;
49    GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);    GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);
50    GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);    GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);
51    GSL_SET_COMPLEX (&zarr[2], 0.0, 0.0);    GSL_SET_COMPLEX (&zarr[2], 0.0, 0.0);
# Line 98  gsl_poly_solve_quartic (double a, double Line 93  gsl_poly_solve_quartic (double a, double
93    
94    /* For non-degenerate solutions, proceed by constructing and    /* For non-degenerate solutions, proceed by constructing and
95     * solving the resolvent cubic */     * solving the resolvent cubic */
96    aa = a * a;    double aa = a * a;
97    pp = b - q1 * aa;    double pp = b - (3.0/8.0) * aa;
98    qq = c - q2 * a * (b - q4 * aa);    double qq = c - (1.0/2.0) * a * (b - (1.0/4.0) * aa);
99    rr = d - q4 * (a * c - q4 * aa * (b - q3 * aa));    double rr = d - (1.0/4.0) * (a * c - (1.0/4.0) * aa * (b - (3.0/16.0) * aa));
100    rc = q2 * pp;    double rc = (1.0/2.0) * pp;
101    sc = q4 * (q4 * pp * pp - rr);    double sc = (1.0/4.0) * ((1.0/4.0) * pp * pp - rr);
102    tc = -(q8 * qq * q8 * qq);    double tc = -((1.0/8.0) * qq * (1.0/8.0) * qq);
103    
104    /* This code solves the resolvent cubic in a convenient fashion    /* This code solves the resolvent cubic in a convenient fashion
105     * for this implementation of the quartic. If there are three real     * for this implementation of the quartic. If there are three real
# Line 114  gsl_poly_solve_quartic (double a, double Line 109  gsl_poly_solve_quartic (double a, double
109     * u[1] and u[2], respectively. Additionally, this     * u[1] and u[2], respectively. Additionally, this
110     * calculates the discriminant of the cubic and puts it into the     * calculates the discriminant of the cubic and puts it into the
111     * variable disc. */     * variable disc. */
112      double disc;
113    {    {
114      double qcub = (rc * rc - 3 * sc);      double qcub = (rc * rc - 3 * sc);
115      double rcub = (2 * rc * rc * rc - 9 * rc * sc + 27 * tc);      double rcub = (2 * rc * rc * rc - 9 * rc * sc + 27 * tc);
# Line 186  gsl_poly_solve_quartic (double a, double Line 182  gsl_poly_solve_quartic (double a, double
182     * mt=2 : 0 real roots (disc < 0)     * mt=2 : 0 real roots (disc < 0)
183     * mt=3 : 2 real roots (disc > 0)     * mt=3 : 2 real roots (disc > 0)
184     */     */
185      double mt;
186    if (0 == disc)    if (0 == disc)
187      {      {
188        u[2] = u[1];        u[2] = u[1];
# Line 197  gsl_poly_solve_quartic (double a, double Line 194  gsl_poly_solve_quartic (double a, double
194        v[1] = fabs (u[1]);        v[1] = fabs (u[1]);
195        v[2] = fabs (u[2]);        v[2] = fabs (u[2]);
196        v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]);        v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]);
197    
198          int k1 = 0, k2 = 0;
199        if (v1 == v[0])        if (v1 == v[0])
200          {          {
201            k1 = 0;            k1 = 0;
# Line 238  gsl_poly_solve_quartic (double a, double Line 237  gsl_poly_solve_quartic (double a, double
237      }      }
238    
239    /* Solve the quadratic to obtain the roots to the quartic */    /* Solve the quadratic to obtain the roots to the quartic */
240    q = qq;    double q = qq; /* ! */
241    if (0.0 != gsl_complex_abs (gsl_complex_mul (w1, w2)))    if (0.0 != gsl_complex_abs (gsl_complex_mul (w1, w2)))
242      {      {
243        w3 =        w3 =
244          gsl_complex_mul_real (gsl_complex_inverse (gsl_complex_mul (w1, w2)),          gsl_complex_mul_real (gsl_complex_inverse (gsl_complex_mul (w1, w2)),
245                                -q / 8.0);                                -q / 8.0);
246      }      }
247    h = r4 * a;    double h = a / 4.0;
248    zarr[0] =    zarr[0] =
249      gsl_complex_add_real (gsl_complex_add (gsl_complex_add (w1, w2), w3), -h);      gsl_complex_add_real (gsl_complex_add (gsl_complex_add (w1, w2), w3), -h);
250    zarr[1] =    zarr[1] =

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

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