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); |
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 |
{ |
{ |
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); |
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; |
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; |
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 */ |
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 |
{ |
{ |
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 |
{ |
{ |
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); |
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 */ |
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)) |
383 |
SWAP (*z2, *z3); |
SWAP (*z2, *z3); |
384 |
} |
} |
385 |
} |
} |
|
|
|
386 |
} |
} |
387 |
|
|
388 |
return 4; |
return 4; |