ROOT   Reference Guide
Searching...
No Matches
complex_quartic.h
Go to the documentation of this file.
1// @(#)root/mathmore:$Id$
2// Authors: L. Moneta, A. Zsenei 08/2005
3
4/* poly/zsolve_quartic.c
5 *
6 * Copyright (C) 2003 CERN and K.S. K\"{o}lbig
7 *
8 * Converted to C and implemented into the GSL Library - Sept. 2003
9 * by Andrew W. Steiner and Andy Buckley
10 *
11 * This program is free software; you can redistribute it and/or modify
13 * the Free Software Foundation; either version 2 of the License, or (at
14 * your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful, but
17 * WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program; if not, write to the Free Software
23 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 */
25
26/* zsolve_quartic.c - finds the complex roots of
27 * x^4 + a x^3 + b x^2 + c x + d = 0
28 */
29
30#include <cmath>
31#include <gsl/gsl_math.h>
32#include <gsl/gsl_complex.h>
33#include <gsl/gsl_complex_math.h>
34#include <gsl/gsl_poly.h>
35
36#define SWAP(a,b) do { gsl_complex tmp = b ; b = a ; a = tmp ; } while(0)
37
38int
39gsl_poly_complex_solve_quartic (double a, double b, double c, double d,
40 gsl_complex * z0, gsl_complex * z1,
41 gsl_complex * z2, gsl_complex * z3)
42{
43 /* Relative tolerance below which disc is treated as zero. */
44 static constexpr double ZeroTolerance = 1e-15;
45
46 gsl_complex i, zarr[4], w1, w2, w3;
47 double r4 = 1.0 / 4.0;
48 double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0;
49 double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0;
50 double u[3], v[3], v1, v2, disc, disc_tol;
51 double aa, pp, qq, rr, rc, sc, tc, q, h;
52 int k1 = 0, k2 = 0, mt;
53
54 GSL_SET_COMPLEX (&i, 0.0, 1.0);
55 GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);
56 GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);
57 GSL_SET_COMPLEX (&zarr[2], 0.0, 0.0);
58 GSL_SET_COMPLEX (&zarr[3], 0.0, 0.0);
59 GSL_SET_COMPLEX (&w1, 0.0, 0.0);
60 GSL_SET_COMPLEX (&w2, 0.0, 0.0);
61 GSL_SET_COMPLEX (&w3, 0.0, 0.0);
62
63 /* Deal easily with the cases where the quartic is degenerate. The
64 * ordering of solutions is done explicitly. */
65 if (0 == b && 0 == c)
66 {
67 if (0 == d)
68 {
69 if (a > 0)
70 {
71 GSL_SET_COMPLEX (z0, -a, 0.0);
72 GSL_SET_COMPLEX (z1, 0.0, 0.0);
73 GSL_SET_COMPLEX (z2, 0.0, 0.0);
74 GSL_SET_COMPLEX (z3, 0.0, 0.0);
75 }
76 else
77 {
78 GSL_SET_COMPLEX (z0, 0.0, 0.0);
79 GSL_SET_COMPLEX (z1, 0.0, 0.0);
80 GSL_SET_COMPLEX (z2, 0.0, 0.0);
81 GSL_SET_COMPLEX (z3, -a, 0.0);
82 }
83 return 4;
84 }
85 else if (0 == a)
86 {
87 if (d > 0)
88 {
89 double sqrt_d = sqrt (d);
90 gsl_complex i_sqrt_d = gsl_complex_mul_real (i, sqrt_d);
91 gsl_complex minus_i = gsl_complex_conjugate (i);
92 *z3 = gsl_complex_sqrt (i_sqrt_d);
93 *z2 = gsl_complex_mul (minus_i, *z3);
94 *z1 = gsl_complex_negative (*z2);
95 *z0 = gsl_complex_negative (*z3);
96 }
97 else
98 {
99 double sqrt_abs_d = sqrt (-d);
100 *z3 = gsl_complex_sqrt_real (sqrt_abs_d);
101 *z2 = gsl_complex_mul (i, *z3);
102 *z1 = gsl_complex_negative (*z2);
103 *z0 = gsl_complex_negative (*z3);
104 }
105 return 4;
106 }
107 }
108
109 if (0.0 == c && 0.0 == d)
110 {
111 disc = (a * a - 4.0 * b);
112 if (disc < 0.0)
113 {
114 mt = 3;
115 }
116 else
117 {
118 mt = 1;
119 }
120 *z0 = zarr[0];
121 *z1 = zarr[0];
122 gsl_poly_complex_solve_quadratic (1.0, a, b, z2, z3);
123 }
124 else
125 {
126 /* For non-degenerate solutions, proceed by constructing and
127 * solving the resolvent cubic */
128 aa = a * a;
129 pp = b - q1 * aa;
130 qq = c - q2 * a * (b - q4 * aa);
131 rr = d - q4 * (a * c - q4 * aa * (b - q3 * aa));
132 rc = q2 * pp;
133 sc = q4 * (q4 * pp * pp - rr);
134 tc = -(q8 * qq * q8 * qq);
135
136 /* This code solves the resolvent cubic in a convenient fashion
137 * for this implementation of the quartic. If there are three real
138 * roots, then they are placed directly into u[]. If two are
139 * complex, then the real root is put into u[0] and the real
140 * and imaginary part of the complex roots are placed into
141 * u[1] and u[2], respectively. Additionally, this
142 * calculates the discriminant of the cubic and puts it into the
143 * variable disc. */
144 {
145 double qcub = (rc * rc - 3 * sc);
146 double rcub = (2 * rc * rc * rc - 9 * rc * sc + 27 * tc);
147
148 double Q = qcub / 9;
149 double R = rcub / 54;
150
151 double Q3 = Q * Q * Q;
152 double R2 = R * R;
153
154 disc = R2 - Q3;
155 disc_tol = std::max(ZeroTolerance, ZeroTolerance * (fabs(R2) + fabs(Q3)));
156
157
158 if (0 == R && 0 == Q)
159 {
160 u[0] = -rc / 3;
161 u[1] = -rc / 3;
162 u[2] = -rc / 3;
163 }
164 else if (fabs(disc) < disc_tol)
165 {
166 // R2 and Q3 are considered equal here
167 double sqrtQ = sqrt (Q);
168 if (R > 0)
169 {
170 u[0] = -2 * sqrtQ - rc / 3;
171 u[1] = sqrtQ - rc / 3;
172 u[2] = sqrtQ - rc / 3;
173 }
174 else
175 {
176 u[0] = -sqrtQ - rc / 3;
177 u[1] = -sqrtQ - rc / 3;
178 u[2] = 2 * sqrtQ - rc / 3;
179 }
180 }
181 else if ( R2 < Q3)
182 {
183 double sqrtQ = sqrt (Q);
184 double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
185 double ctheta = R / sqrtQ3;
186 double theta = 0;
187 // protect against numerical error can make this larger than one
188 if ( fabs(ctheta) < 1.0 )
189 theta = acos( ctheta);
190 else if ( ctheta <= -1.0)
191 theta = M_PI;
192
193 double norm = -2 * sqrtQ;
194
195 u[0] = norm * cos (theta / 3) - rc / 3;
196 u[1] = norm * cos ((theta + 2.0 * M_PI) / 3) - rc / 3;
197 u[2] = norm * cos ((theta - 2.0 * M_PI) / 3) - rc / 3;
198 }
199 else
200 {
201 double sgnR = (R >= 0 ? 1 : -1);
202 double modR = fabs (R);
203 double sqrt_disc = sqrt (disc);
204 double A = -sgnR * pow (modR + sqrt_disc, 1.0 / 3.0);
205 double B = Q / A;
206
207 double mod_diffAB = fabs (A - B);
208 u[0] = A + B - rc / 3;
209 u[1] = -0.5 * (A + B) - rc / 3;
210 u[2] = -(sqrt (3.0) / 2.0) * mod_diffAB;
211 }
212 }
213 /* End of solution to resolvent cubic */
214
215 /* Combine the square roots of the roots of the cubic
216 * resolvent appropriately. Also, calculate 'mt' which
217 * designates the nature of the roots:
218 * mt=1 : 4 real roots
219 * mt=2 : 0 real roots
220 * mt=3 : 2 real roots
221 */
222 // when disc == 0 2 roots are identicals
223 if (disc_tol >= disc)
224 {
225 mt = 2;
226 v[0] = fabs (u[0]);
227 v[1] = fabs (u[1]);
228 v[2] = fabs (u[2]);
229
230 v1 = GSL_MAX (GSL_MAX (v[0], v[1]), v[2]);
231 if (v1 == v[0])
232 {
233 k1 = 0;
234 v2 = GSL_MAX (v[1], v[2]);
235 }
236 else if (v1 == v[1])
237 {
238 k1 = 1;
239 v2 = GSL_MAX (v[0], v[2]);
240 }
241 else
242 {
243 k1 = 2;
244 v2 = GSL_MAX (v[0], v[1]);
245 }
246
247 if (v2 == v[0])
248 {
249 k2 = 0;
250 }
251 else if (v2 == v[1])
252 {
253 k2 = 1;
254 }
255 else
256 {
257 k2 = 2;
258 }
259 w1 = gsl_complex_sqrt_real (u[k1]);
260 w2 = gsl_complex_sqrt_real (u[k2]);
261 }
262 else
263 {
264 mt = 3;
265 GSL_SET_COMPLEX (&w1, u[1], u[2]);
266 GSL_SET_COMPLEX (&w2, u[1], -u[2]);
267 w1 = gsl_complex_sqrt (w1);
268 w2 = gsl_complex_sqrt (w2);
269 }
270 /* Solve the quadratic in order to obtain the roots
271 * to the quartic */
272 q = qq;
273 gsl_complex prod_w = gsl_complex_mul (w1, w2);
274 //gsl_complex mod_prod_w = gsl_complex_abs (prod_w);
275 /*
276 Changed from gsl_complex to double in order to make it compile.
277 */
278 double mod_prod_w = gsl_complex_abs (prod_w);
279 if (0.0 != mod_prod_w)
280 {
281 gsl_complex inv_prod_w = gsl_complex_inverse (prod_w);
282 w3 = gsl_complex_mul_real (inv_prod_w, -q / 8.0);
283 }
284
285 h = r4 * a;
286 gsl_complex sum_w12 = gsl_complex_add (w1, w2);
287 gsl_complex neg_sum_w12 = gsl_complex_negative (sum_w12);
288 gsl_complex sum_w123 = gsl_complex_add (sum_w12, w3);
289 gsl_complex neg_sum_w123 = gsl_complex_add (neg_sum_w12, w3);
290
291 gsl_complex diff_w12 = gsl_complex_sub (w2, w1);
292 gsl_complex neg_diff_w12 = gsl_complex_negative (diff_w12);
293 gsl_complex diff_w123 = gsl_complex_sub (diff_w12, w3);
294 gsl_complex neg_diff_w123 = gsl_complex_sub (neg_diff_w12, w3);
295
296 zarr[0] = gsl_complex_add_real (sum_w123, -h);
297 zarr[1] = gsl_complex_add_real (neg_sum_w123, -h);
298 zarr[2] = gsl_complex_add_real (diff_w123, -h);
299 zarr[3] = gsl_complex_add_real (neg_diff_w123, -h);
300
301 /* Arrange the roots into the variables z0, z1, z2, z3 */
302 if (2 == mt)
303 {
304 if (u[k1] >= 0 && u[k2] >= 0)
305 {
306 mt = 1;
307 GSL_SET_COMPLEX (z0, GSL_REAL (zarr[0]), 0.0);
308 GSL_SET_COMPLEX (z1, GSL_REAL (zarr[1]), 0.0);
309 GSL_SET_COMPLEX (z2, GSL_REAL (zarr[2]), 0.0);
310 GSL_SET_COMPLEX (z3, GSL_REAL (zarr[3]), 0.0);
311 }
312 else if (u[k1] >= 0 && u[k2] < 0)
313 {
314 *z0 = zarr[0];
315 *z1 = zarr[3];
316 *z2 = zarr[2];
317 *z3 = zarr[1];
318 }
319 else if (u[k1] < 0 && u[k2] >= 0)
320 {
321 *z0 = zarr[0];
322 *z1 = zarr[2];
323 *z2 = zarr[3];
324 *z3 = zarr[1];
325 }
326 else if (u[k1] < 0 && u[k2] < 0)
327 {
328 *z0 = zarr[0];
329 *z1 = zarr[1];
330 *z2 = zarr[3];
331 *z3 = zarr[2];
332 }
333 }
334 else if (3 == mt)
335 {
336 GSL_SET_COMPLEX (z0, GSL_REAL (zarr[0]), 0.0);
337 GSL_SET_COMPLEX (z1, GSL_REAL (zarr[1]), 0.0);
338 *z2 = zarr[3];
339 *z3 = zarr[2];
340 }
341 }
342
343 /*
344 * Sort the roots as usual: main sorting by ascending real part, secondary
345 * sorting by ascending imaginary part
346 */
347
348 if (1 == mt)
349 {
350 /* Roots are all real, sort them by the real part */
351 if (GSL_REAL (*z0) > GSL_REAL (*z1)) SWAP (*z0, *z1);
352 if (GSL_REAL (*z0) > GSL_REAL (*z2)) SWAP (*z0, *z2);
353 if (GSL_REAL (*z0) > GSL_REAL (*z3)) SWAP (*z0, *z3);
354
355 if (GSL_REAL (*z1) > GSL_REAL (*z2)) SWAP (*z1, *z2);
356 if (GSL_REAL (*z2) > GSL_REAL (*z3))
357 {
358 SWAP (*z2, *z3);
359 if (GSL_REAL (*z1) > GSL_REAL (*z2)) SWAP (*z1, *z2);
360 }
361 }
362 else if (2 == mt)
363 {
364 /* Roots are all complex. z0 and z1 are conjugates
365 * and z2 and z3 are conjugates. Sort the real parts first */
366 if (GSL_REAL (*z0) > GSL_REAL (*z2))
367 {
368 SWAP (*z0, *z2);
369 SWAP (*z1, *z3);
370 }
371 /* Then sort by the imaginary parts */
372 if (GSL_IMAG (*z0) > GSL_IMAG (*z1)) SWAP (*z0, *z1);
373 if (GSL_IMAG (*z2) > GSL_IMAG (*z3)) SWAP (*z2, *z3);
374 }
375 else
376 {
377 /* 2 real roots. z2 and z3 are conjugates. */
378
379 /* Swap complex roots */
380 if (GSL_IMAG (*z2) > GSL_IMAG (*z3)) SWAP (*z2, *z3);
381
382 /* Sort real parts */
383 if (GSL_REAL (*z0) > GSL_REAL (*z1)) SWAP (*z0, *z1);
384 if (GSL_REAL (*z1) > GSL_REAL (*z2))
385 {
386 if (GSL_REAL (*z0) > GSL_REAL (*z2))
387 {
388 SWAP (*z0, *z2);
389 SWAP (*z1, *z3);
390 }
391 else
392 {
393 SWAP (*z1, *z2);
394 SWAP (*z2, *z3);
395 }
396 }
397 }
398
399 return 4;
400}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
#define R(a, b, c, d, e, f, g, h, i)
Definition RSha256.hxx:110
#define e(i)
Definition RSha256.hxx:103
#define GSL_MAX(a, b)
#define M_PI
Definition Rotated.cxx:105
float * q
int gsl_poly_complex_solve_quartic(double a, double b, double c, double d, gsl_complex *z0, gsl_complex *z1, gsl_complex *z2, gsl_complex *z3)
#define SWAP(a, b)
#define R2(v, w, x, y, z, i)
Definition sha1.inl:137