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> |
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); |
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 |
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); |
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]; |
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; |
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] = |