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; |
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 |
{ |
{ |
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 |
} |
} |
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) |