ROOT   6.10/09 Reference Guide
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
12  * it under the terms of the GNU General Public License as published by
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 <math.h>
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
38 int
39 gsl_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  gsl_complex i, zarr[4], w1, w2, w3;
44  double r4 = 1.0 / 4.0;
45  double q2 = 1.0 / 2.0, q4 = 1.0 / 4.0, q8 = 1.0 / 8.0;
46  double q1 = 3.0 / 8.0, q3 = 3.0 / 16.0;
47  double u[3], v[3], v1, v2, disc;
48  double aa, pp, qq, rr, rc, sc, tc, q, h;
49  int k1 = 0, k2 = 0, mt;
50
51  GSL_SET_COMPLEX (&i, 0.0, 1.0);
52  GSL_SET_COMPLEX (&zarr[0], 0.0, 0.0);
53  GSL_SET_COMPLEX (&zarr[1], 0.0, 0.0);
54  GSL_SET_COMPLEX (&zarr[2], 0.0, 0.0);
55  GSL_SET_COMPLEX (&zarr[3], 0.0, 0.0);
56  GSL_SET_COMPLEX (&w1, 0.0, 0.0);
57  GSL_SET_COMPLEX (&w2, 0.0, 0.0);
58  GSL_SET_COMPLEX (&w3, 0.0, 0.0);
59
60  /* Deal easily with the cases where the quartic is degenerate. The
61  * ordering of solutions is done explicitly. */
62  if (0 == b && 0 == c)
63  {
64  if (0 == d)
65  {
66  if (a > 0)
67  {
68  GSL_SET_COMPLEX (z0, -a, 0.0);
69  GSL_SET_COMPLEX (z1, 0.0, 0.0);
70  GSL_SET_COMPLEX (z2, 0.0, 0.0);
71  GSL_SET_COMPLEX (z3, 0.0, 0.0);
72  }
73  else
74  {
75  GSL_SET_COMPLEX (z0, 0.0, 0.0);
76  GSL_SET_COMPLEX (z1, 0.0, 0.0);
77  GSL_SET_COMPLEX (z2, 0.0, 0.0);
78  GSL_SET_COMPLEX (z3, -a, 0.0);
79  }
80  return 4;
81  }
82  else if (0 == a)
83  {
84  if (d > 0)
85  {
86  double sqrt_d = sqrt (d);
87  gsl_complex i_sqrt_d = gsl_complex_mul_real (i, sqrt_d);
88  gsl_complex minus_i = gsl_complex_conjugate (i);
89  *z3 = gsl_complex_sqrt (i_sqrt_d);
90  *z2 = gsl_complex_mul (minus_i, *z3);
91  *z1 = gsl_complex_negative (*z2);
92  *z0 = gsl_complex_negative (*z3);
93  }
94  else
95  {
96  double sqrt_abs_d = sqrt (-d);
97  *z3 = gsl_complex_sqrt_real (sqrt_abs_d);
98  *z2 = gsl_complex_mul (i, *z3);
99  *z1 = gsl_complex_negative (*z2);
100  *z0 = gsl_complex_negative (*z3);
101  }
102  return 4;
103  }
104  }
105
106  if (0.0 == c && 0.0 == d)
107  {
108  disc = (a * a - 4.0 * b);
109  if (disc < 0.0)
110  {
111  mt = 3;
112  }
113  else
114  {
115  mt = 1;
116  }
117  *z0 = zarr[0];
118  *z1 = zarr[0];
119  gsl_poly_complex_solve_quadratic (1.0, a, b, z2, z3);
120  }
121  else
122  {
123  /* For non-degenerate solutions, proceed by constructing and
124  * solving the resolvent cubic */
125  aa = a * a;
126  pp = b - q1 * aa;
127  qq = c - q2 * a * (b - q4 * aa);
128  rr = d - q4 * (a * c - q4 * aa * (b - q3 * aa));
129  rc = q2 * pp;
130  sc = q4 * (q4 * pp * pp - rr);
131  tc = -(q8 * qq * q8 * qq);
132
133  /* This code solves the resolvent cubic in a convenient fashion
134  * for this implementation of the quartic. If there are three real
135  * roots, then they are placed directly into u[]. If two are
136  * complex, then the real root is put into u[0] and the real
137  * and imaginary part of the complex roots are placed into
138  * u[1] and u[2], respectively. Additionally, this
139  * calculates the discriminant of the cubic and puts it into the
140  * variable disc. */
141  {
142  double qcub = (rc * rc - 3 * sc);
143  double rcub = (2 * rc * rc * rc - 9 * rc * sc + 27 * tc);
144
145  double Q = qcub / 9;
146  double R = rcub / 54;
147
148  double Q3 = Q * Q * Q;
149  double R2 = R * R;
150
151  disc = R2 - Q3;
152
153 // more numerical problems with this calculation of disc
154 // double CR2 = 729 * rcub * rcub;
155 // double CQ3 = 2916 * qcub * qcub * qcub;
156 // disc = (CR2 - CQ3) / 2125764.0;
157
158
159  if (0 == R && 0 == Q)
160  {
161  u[0] = -rc / 3;
162  u[1] = -rc / 3;
163  u[2] = -rc / 3;
164  }
165  else if (R2 == Q3)
166  {
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 (0 >= 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 }
401
static double B[]
TH1 * h
Definition: legend2.C:5
TArc * a
Definition: textangle.C:12
static double A[]
double cos(double)
double sqrt(double)
double acos(double)
double pow(double, double)
#define SWAP(a, b)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
#define M_PI
Definition: Rotated.cxx:105
SVector< double, 2 > v
Definition: Dict.h:5
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)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
float * q
Definition: THbookFile.cxx:87
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
TRandom3 R
a TMatrixD.
Definition: testIO.cxx:28
static double Q[]