PLUTO
Macros | Functions
quartic.old.c File Reference

Solve quartic and cubic equations-. More...

#include "pluto.h"
#include "complex.h"
Include dependency graph for quartic.old.c:

Go to the source code of this file.

Macros

#define swap(x, y)   f = x; x = y; y = f;
 
#define ZQUARTIC(x)    (ez + (x)*(dz + (x)*(cz + (x)*(bz + (x)))))
 

Functions

int QuarticSolve (double b, double c, double d, double e, double *z)
 
int CubicSolve (double b, double c, double d, double z[])
 
int QuarticSolve2 (double b, double c, double d, double e, double *z)
 

Detailed Description

Solve quartic and cubic equations-.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Sept 10, 2012

Definition in file quartic.old.c.

Macro Definition Documentation

#define swap (   x,
  y 
)    f = x; x = y; y = f;

Definition at line 18 of file quartic.old.c.

#define ZQUARTIC (   x)    (ez + (x)*(dz + (x)*(cz + (x)*(bz + (x)))))

Function Documentation

int CubicSolve ( double  b,
double  c,
double  d,
double  z[] 
)

Solve a cubic equation in the form

\[ z^3 + bz^2 + cz + d = 0 \]

For its purpose, it is assumed that ALL roots are double. This makes things faster.

Parameters
[in]bcoefficient of the cubic
[in]ccoefficient of the cubic
[in]dcoefficient of the cubic
[out]zvector containing the roots of the cubic. Roots should be sorted in increasing order.

Reference: http://www.1728.com/cubic2.htm

Returns
Return 0 on success.

Definition at line 147 of file quartic.old.c.

167 {
168  double b2, g2;
169  double f, g, h;
170  double i, i2, j, k, m, n, p;
171  static double one_3 = 1.0/3.0, one_27=1.0/27.0;
172 
173  b2 = b*b;
174 
175 /* ----------------------------------------------
176  the expression for f should be
177  f = c - b*b/3.0; however, to avoid negative
178  round-off making h > 0.0 or g^2/4 - h < 0.0
179  we let c --> c(1- 1.1e-16)
180  ---------------------------------------------- */
181 
182  f = c*(1.0 - 1.e-16) - b2*one_3;
183  g = b*(2.0*b2 - 9.0*c)*one_27 + d;
184  g2 = g*g;
185  i2 = -f*f*f*one_27;
186  h = g2*0.25 - i2;
187 
188 /* --------------------------------------------
189  double roots are possible only when
190 
191  h <= 0
192  -------------------------------------------- */
193 
194  if (h > 1.e-12){
195 // printf ("Only one double root (%12.6e)!\n", h);
196  }
197  if (i2 < 0.0){
198 /*
199  printf ("i2 < 0.0 %12.6e\n",i2);
200  return(1);
201 */
202  i2 = 0.0;
203  }
204 
205 /* --------------------------------------
206  i^2 must be >= g2*0.25
207  -------------------------------------- */
208 
209  i = sqrt(i2); /* > 0 */
210  j = pow(i, one_3); /* > 0 */
211  k = -0.5*g/i;
212 
213 /* this is to prevent unpleseant situation
214  where both g and i are close to zero */
215 
216  k = (k < -1.0 ? -1.0:k);
217  k = (k > 1.0 ? 1.0:k);
218 
219  k = acos(k)*one_3; /* pi/3 < k < 0 */
220 
221  m = cos(k); /* > 0 */
222  n = sqrt(3.0)*sin(k); /* > 0 */
223  p = -b*one_3;
224 
225  z[0] = -j*(m + n) + p;
226  z[1] = -j*(m - n) + p;
227  z[2] = 2.0*j*m + p;
228 
229 /* ------------------------------------------------------
230  Since j, m, n > 0, it should follow that from
231 
232  z0 = -jm - jn + p
233  z1 = -jm + jn + p
234  z2 = 2jm + p
235 
236  z2 is the greatest of the roots, while z0 is the
237  smallest one.
238  ------------------------------------------------------ */
239 
240  return(0);
241 }
n
static int n
Definition: analysis.c:3
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
menu.c
tuple c
Definition: menu.py:375
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

int QuarticSolve ( double  b,
double  c,
double  d,
double  e,
double *  z 
)

Solve a quartic equation in the form

\[ z^4 + bz^3 + cz^2 + dz + e = 0 \]

For its purpose, it is assumed that ALL roots are double. This makes things faster.

Parameters
[in]bcoefficient of the quartic
[in]ccoefficient of the quartic
[in]dcoefficient of the quartic
[in]ecoefficient of the quartic
[out]zvector containing the (double) roots of the quartic

Reference:

http://www.1728.com/quartic2.htm

Returns
Return 0 on success, 1 if cubic solver fails, 2 if NaN has been found and 3 if quartic is not satisfied.

Definition at line 21 of file quartic.old.c.

45 {
46  int n, j, ifail;
47  double b2, f, g, h;
48  double a2, a1, a0, u[4];
49  double p, q, r, s;
50  static double three_256 = 3.0/256.0;
51  static double one_64 = 1.0/64.0;
52 
53  b2 = b*b;
54 
55  f = c - b2*0.375;
56  g = d + b2*b*0.125 - b*c*0.5;
57  h = e - b2*b2*three_256 + 0.0625*b2*c - 0.25*b*d;
58 
59  a2 = 0.5*f;
60  a1 = (f*f - 4.0*h)*0.0625;
61  a0 = -g*g*one_64;
62 
63  ifail = CubicSolve(a2, a1, a0, u);
64 
65  if (ifail) return 1;
66 
67  if (u[1] < 1.e-14){
68 
69  p = sqrt(u[2]);
70  s = 0.25*b;
71  z[0] = z[2] = - p - s;
72  z[1] = z[3] = + p - s;
73 
74  }else{
75 
76  p = sqrt(u[1]);
77  q = sqrt(u[2]);
78 
79  r = -0.125*g/(p*q);
80  s = 0.25*b;
81 
82  z[0] = - p - q + r - s;
83  z[1] = p - q - r - s;
84  z[2] = - p + q - r - s;
85  z[3] = p + q + r - s;
86 
87  }
88 
89 /* ----------------------------------------
90  Sort roots in ascending order.
91  Since z[0] < z[3], we first order
92 
93  z[0] < z[1] < z[3]
94 
95  and then put z[2] into the sequence.
96  ---------------------------------------- */
97 
98  if (z[1] < z[0]) {swap(z[1],z[0]);}
99  else if (z[1] > z[3]) {swap(z[1],z[3]);}
100 
101 
102  if (z[2] < z[0]) {
103  swap(z[2],z[0]);
104  swap(z[2],z[1]);
105  }else if (z[2] < z[1]){
106  swap(z[2],z[1]);
107  }else if (z[2] > z[3]){
108  swap(z[2],z[3]);
109  }
110 
111 /*
112  for (j = 1; j < 4; j++){
113  f = z[j];
114  n = j - 1;
115  while(n >= 0 && z[n] > f) {
116  z[n+1] = z[n];
117  n--;
118  }
119  z[n+1] = f;
120  }
121 */
122 
123  /* ----------------------------------------------
124  verify that cmax and cmin satisfy original
125  equation
126  ---------------------------------------------- */
127 
128  for (n = 0; n < 4; n++){
129  s = e + z[n]*(d + z[n]*(c + z[n]*(b + z[n])));
130  if (s != s) {
131 // print ("! QuarticSolve: NaN found.\n");
132  return 2;
133  }
134  if (fabs(s) > 1.e-6) {
135 // print ("! QuarticSolve: solution does not satisfy f(z) = 0; f(z) = %12.6e\n",s);
136  return 3;
137  }
138  }
139 
140  return(0);
141 /*
142  printf (" z: %f ; %f ; %f ; %f\n",z[0], z[1], z[2], z[3]);
143  */
144 }
n
static int n
Definition: analysis.c:3
CubicSolve
int CubicSolve(double b, double c, double d, double z[])
Definition: quartic.old.c:147
swap
#define swap(x, y)
Definition: quartic.old.c:18
j
int j
Definition: analysis.c:2
menu.c
tuple c
Definition: menu.py:375
s
#define s
q
static Runtime q
Definition: runtime_setup.c:490

Here is the call graph for this function:

Here is the caller graph for this function:

int QuarticSolve2 ( double  b,
double  c,
double  d,
double  e,
double *  z 
)

Solve a quartic equation in the form

\[ z^4 + bz^3 + cz^2 + dz + e = 0 \]

using Durand-Kerner method.

Parameters
[in]bcoefficient of the quartic
[in]ccoefficient of the quartic
[in]dcoefficient of the quartic
[in]ecoefficient of the quartic
[out]zvector containing the (double) roots of the quartic

Reference:

Returns
Return 0 on success, 1 if cubic solver fails, 2 if NaN has been found and 3 if quartic is not satisfied.

Definition at line 247 of file quartic.old.c.

270 {
271  int k;
272  double complex bz = b;
273  double complex cz = c;
274  double complex dz = d;
275  double complex ez = e;
276  double complex p, q, r, s;
277  double complex dp, dq, dr, ds;
278  double acc = 1.e-5, err;
279 
280  p = -1.0 + 0.000012*_Complex_I;
281  q = -0.25 + 0.00002*_Complex_I;
282  r = 0.25 + 0.00002*_Complex_I;
283  s = 1.0 - 0.000011*_Complex_I;
284 
285  for (k = 0; k < 4096; k++){
286 
287  dp = ZQUARTIC(p)/( (p - q)*(p - r)*(p - s) );
288  p -= dp;
289 
290  dq = ZQUARTIC(q)/( (q - p)*(q - r)*(q - s) );
291  q -= dq;
292 
293  dr = ZQUARTIC(r)/( (r - p)*(r - q)*(r - s) );
294  r -= dr;
295 
296  ds = ZQUARTIC(s)/( (s - p)*(s - q)*(s - r) );
297  s -= ds;
298 
299  err = MAX(cabs(dp), cabs(dq));
300  err = MAX(err, cabs(dr));
301  err = MAX(err, cabs(ds));
302 
303  if (err < acc) break;
304 
305  }
306 /*
307  printf ("Solution = %12.6e %12.6e %12.6e %12.6e\n",
308  creal(p), creal(q), creal(r), creal(s));
309  printf ("Number of iterations = %d\n",k);
310  printf ("Residual: %8.3e, %8.3e, %8.3e, %8.3e\n",cabs(dp),cabs(dq),cabs(dr),cabs(ds));
311 */
312 
313  z[0] = creal(p);
314  z[1] = creal(q);
315  z[2] = creal(r);
316  z[3] = creal(s);
317 
318 double zmin, zmax;
319 
320  zmax = MAX(z[0],z[1]);
321  zmax = MAX(zmax,z[2]);
322  zmax = MAX(zmax,z[3]);
323 
324  zmin = MIN(z[0],z[1]);
325  zmin = MIN(zmin,z[2]);
326  zmin = MIN(zmin,z[3]);
327  z[0] = zmin;
328  z[3] = zmax;
329 
330 }
MAX
#define MAX(a, b)
Definition: macros.h:101
MIN
#define MIN(a, b)
Definition: macros.h:104
k
int k
Definition: analysis.c:2
menu.c
tuple c
Definition: menu.py:375
s
#define s
ZQUARTIC
#define ZQUARTIC(x)
q
static Runtime q
Definition: runtime_setup.c:490