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); |
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 |
{ |
{ |
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 |
{ |
{ |
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; |
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; |
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 */ |
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]); |
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)), |
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) |