18 int CubicSolve2 (
double b,
double c,
double d,
double z[]);
83 double sQ, arg, theta, cs, sn, p;
84 const double one_3 = 1.0/3.0;
85 const double one_27 = 1.0/27.0;
86 const double sqrt3 = sqrt(3.0);
97 Q = b2*one_3 - c*(1.0 - 1.e-16);
98 R = b*(2.0*b2 - 9.0*
c)*one_27 + d;
124 arg =
MAX(-1.0, arg);
125 arg =
MIN( 1.0, arg);
128 theta = acos(arg)*one_3;
131 sn = sqrt3*sin(theta);
134 z[0] = -sQ*(cs + sn) + p;
135 z[1] = -sQ*(cs - sn) + p;
136 z[2] = 2.0*sQ*cs + p;
141 print (
"===========================================================\n");
142 print (
"> Resolvent cubic:\n");
143 print (
" g(x) = %18.12e + x*(%18.12e + x*(%18.12e + x))\n", d, c, b);
144 print (
" Q = %8.3e\n",Q);
145 print (
" arg-1 = %8.3e\n", -1.5*R/(Q*sQ)-1.0);
147 print (
"> Cubic roots = %8.3e %8.3e %8.3e\n",z[0],z[1],z[2]);
148 for (l = 0; l < 3; l++){
151 f = d + x*(c + x*(b + x));
152 print (
" verify: g(x[%d]) = %8.3e\n",l,f);
155 print (
"===========================================================\n");
183 double a2, a1, a0, u[4];
185 const double three_256 = 3.0/256.0;
186 const double one_64 = 1.0/64.0;
187 double sz1, sz2, sz3, sz4;
198 q = d + b2*b*0.125 - b*c*0.5;
199 r = e - 3.0*b2*b2/256.0 + b2*c/16.0 - 0.25*b*d;
206 if (ifail != 0)
return 1;
208 u[0] =
MAX(u[0],0.0);
209 u[1] =
MAX(u[1],0.0);
210 u[2] =
MAX(u[2],0.0);
213 if (u[0] != u[0] || u[1] != u[1] || u[2] != u[2])
return 1;
220 z[0] = -0.25*b + sz1 + sz2 + sz3;
221 z[1] = -0.25*b + sz1 - sz2 - sz3;
222 z[2] = -0.25*b - sz1 + sz2 - sz3;
223 z[3] = -0.25*b - sz1 - sz2 + sz3;
227 print (
"Quartic roots = %f %f %f %f; q = %8.3e\n",z[0],z[1],z[2],z[3],q);
241 double x, f,
fp, fm, fmax;
245 x = 1.0; fp = e + x*(d + x*(c + x*(b + x)));
246 x = -1.0; fm = e + x*(d + x*(c + x*(b + x)));
247 for (n = 0; n < 4; n++){
249 f = e + x*(d + x*(c + x*(b + x)));
250 print (
"> CheckSolutions(): f(z[%d]) = %8.3e\n",n,f);
271 print (
"! QuadraticSolve(): del = %8.3e, resetting to 0\n",del);
276 print (
"! QuadraticSolve(): complex roots, del = %8.3e\n",del);
285 sb = (b > 0.0 ? 1.0:-1.0);
286 q = -0.5*(b + sb*sqrt(del));
305 print (
" f(x) = %18.12e + x*(%18.12e + x*(%18.12e ",
307 print (
" + x*(%18.12e + x*%18.12e)))\n", b, 1.0);
309 print (
" b = %18.14e;\n",b);
310 print (
" c = %18.14e;\n",c);
311 print (
" d = %18.14e;\n",d);
312 print (
" e = %18.14e;\n",e);
315 double a1 = c*c + b*d - 4.0*e;
316 double a0 = -(b*c*d - b*b*e - d*d);
318 print (
" Resolvent cubic:\n");
319 print (
" g(x) = %18.12e + x*(%18.12e + x*(%18.12e + x))\n", a0, a1, a2);
320 double Y = 0.25*b*b + 2.0*d/(b+1.e-40) - c;
321 print (
" b^2/4 = %f\n",0.25*b*b);
322 print (
" Y = %8.3e\n",Y);
334 for (j = 1; j < 4; j++){
337 while(n >= 0 && z[n] > f) {
351 print (
"z = %f %f %f %f\n",z[0],z[1],z[2],z[3]);
358 double Q, R, R_Q, B, A, sQ, th, tn, sn, cn;
360 Q = (b*b - 3.0*
c)/9.0;
361 R = (2.0*b*b*b - 9.0*b*c + 27.0*d)/54.0;
364 if (
debug_print)
print (
"R = %8.3e, Q = %8.3e, R^2/Q^3 - 1.0 = %8.3e\n",R,Q,R_Q*R_Q/Q-1.0);
370 cn = 1.0/sqrt(1.0 + tn*tn);
389 z[0] = -2.0*sQ*cos(th/3.0) - 1.0/3.0;
390 z[1] = -2.0*sQ*cos((th - 2.0*
CONST_PI)/3.0) - b/3.0;
391 z[2] = -2.0*sQ*cos((th + 2.0*
CONST_PI)/3.0) - b/3.0;
394 print(
"! CubicSolve2(): 1 root only, arg = %8.3e !!\n", R_Q*R_Q/Q);
397 A = fabs(R)*(1.0 + sqrt(1.0 - A*A*Q));
398 A = -
DSIGN(R)*pow(A, 1.0/3.0);
401 z[0] = (A + B) - b/3.0;
409 print (
"===========================================================\n");
410 print (
"> Resolvent cubic:\n");
411 print (
" g(x) = %18.12e + x*(%18.12e + x*(%18.12e + x))\n", d, c, b);
412 print (
"> Cubic roots = %8.3e %8.3e %8.3e\n",z[0],z[1],z[2]);
413 for (l = 0; l < 3; l++){
416 f = d + x*(c + x*(b + x));
417 print (
" verify: g(x[%d]) = %8.3e\n",l,f);
420 print (
"===========================================================\n");
int QuarticSolveNew(double b, double c, double d, double e, double *z)
int CubicSolve2(double b, double c, double d, double z[])
void CheckSolutions(double b, double c, double d, double e, double *z)
void QuarticPrintCoeffs(double b, double c, double d, double e)
int CubicSolve(double b, double c, double d, double z[])
void print(const char *fmt,...)
void PrintSolution(double *z)
int QuadraticSolve(double a, double b, double c, double *x)
int QuarticSolve(double b, double c, double d, double e, double *z)
void SortArray(double *z, int n)
#define QUIT_PLUTO(e_code)