1584 lines
38 KiB
C++
1584 lines
38 KiB
C++
#include "Racine.h"
|
|
#include "MathUtil.h"
|
|
|
|
#include <stdio.h>
|
|
#include <cmath>
|
|
|
|
#define NCASES 40
|
|
|
|
// constructeur permettant d'initialiser les variables
|
|
Quartic::Quartic():
|
|
debug(6),iterate(false),
|
|
doub0((double)0),
|
|
doub1((double)1),
|
|
doub2((double)2),
|
|
doub3((double)3),
|
|
doub4((double)4),
|
|
doub5((double)5),
|
|
doub6((double)6),
|
|
doub8((double)8),
|
|
doub9((double)9),
|
|
doub12((double)12),
|
|
doub16((double)16),
|
|
doub24((double)24),
|
|
doub27((double)27),
|
|
doub64((double)64),
|
|
inv2(doub1/doub2),
|
|
inv3(doub1/doub3),
|
|
inv4(doub1/doub4),
|
|
inv8(doub1/doub8),
|
|
inv16(doub1/doub16),
|
|
inv32(doub1/(double)32),
|
|
inv64(doub1/doub64),
|
|
inv128(doub1/(double)128),
|
|
d3o8(doub3/doub8),
|
|
d3o256(doub3/(double)256),
|
|
rt3(sqrt(doub3))
|
|
{ doubtol = doub1;
|
|
for (int j = 1; doub1+doubtol > doub1; ++ j)
|
|
{
|
|
doubtol *= inv2;
|
|
}
|
|
doubtol = sqrt(doubtol);
|
|
|
|
doubmin = inv2;
|
|
for (int j = 1; j <= 100; ++ j)
|
|
{
|
|
doubmin=doubmin*doubmin;
|
|
if ((doubmin*doubmin) <= (doubmin*doubmin*inv2))
|
|
break;
|
|
}
|
|
doubmax=doub1/sqrt(doub2*doubmin);
|
|
};
|
|
|
|
|
|
/****************************************************/
|
|
/*
|
|
// find the errors
|
|
//
|
|
// called by quartictest, docoeff, compare,
|
|
// chris, descartes, ferrari, neumark, yacfraid.
|
|
/****************************************************/
|
|
|
|
double Quartic::errors(double a,double b,double c,double d,double rts[4],double rterr[4],int nrts)
|
|
{
|
|
double deriv,test,worst;
|
|
|
|
worst = doubmax;
|
|
if (nrts > 0)
|
|
{
|
|
worst = doub0;
|
|
for (int k = 0 ; k < nrts ; ++ k )
|
|
{
|
|
test = (((rts[k]+a)*rts[k]+b)*rts[k]+c)*rts[k]+d;
|
|
if (test == doub0) rterr[k] = doub0;
|
|
else
|
|
{
|
|
deriv =
|
|
((doub4*rts[k]+doub3*a)*rts[k]+doub2*b)*rts[k]+c;
|
|
if (deriv != doub0)
|
|
rterr[k] = Dabs(test/deriv);
|
|
else
|
|
{
|
|
deriv = (doub12*rts[k]+doub6*a)*rts[k]+doub2*b;
|
|
if (deriv != doub0)
|
|
rterr[k] = sqrt(Dabs(test/deriv));
|
|
else
|
|
{
|
|
deriv = doub24*rts[k]+doub6*a;
|
|
if (deriv != doub0)
|
|
rterr[k] = curoot(Dabs(test/deriv));
|
|
else
|
|
rterr[k] = sqrt(sqrt(Dabs(test)/doub24));
|
|
}
|
|
}
|
|
}
|
|
if (rts[k] != doub0) rterr[k] /= rts[k];
|
|
if (rterr[k] < doub0) rterr[k] = -rterr[k];
|
|
if (rterr[k] > worst) worst = rterr[k];
|
|
}
|
|
}
|
|
return(worst);
|
|
}; /* errors */
|
|
|
|
/**********************************************/
|
|
// find cos(acos(x)/3)
|
|
//
|
|
// 16 Jul 1981 Don Herbison-Evans
|
|
//
|
|
// called by cubic .
|
|
/**********************************************/
|
|
|
|
double Quartic::acos3(double x)
|
|
{
|
|
double value;
|
|
double y = acos(x);
|
|
|
|
// value = std::cos(acos(x)*inv3); pb de nan
|
|
// on traite dans le cas d'une erreur (qui ça arrive !!)
|
|
// comme ça ne marche pas pour l'instant sur la lo-lg2m-045 on spécifie les machines !!!
|
|
#ifdef SYSTEM_MAC_OS_X
|
|
switch(fpclassify(y))
|
|
{ case FP_NAN: case FP_INFINITE:
|
|
{ // cas d'un nombre infini on fait quelques tests
|
|
if (Dabs(x-1.) <= ConstMath::pasmalpetit)
|
|
{ // dans ce cas on a affaire à une erreur d'arrondi
|
|
y=0.;
|
|
}
|
|
else if (x > 1.)
|
|
{ // il y a un pb dans l'algo, on le signal mais on continu
|
|
cout << "\n erreur dans la fonction Qartic::acos3, l'argument est > 1 , on limite ! "
|
|
<< " x-1= " << (x-1.);
|
|
y=0.;
|
|
};
|
|
};
|
|
break;
|
|
};
|
|
#endif
|
|
#ifdef SYSTEM_MAC_OS_X_unix
|
|
switch(fpclassify(y))
|
|
{ case FP_NAN: case FP_INFINITE:
|
|
{ // cas d'un nombre infini on fait quelques tests
|
|
if (Dabs(x-1.) <= ConstMath::pasmalpetit)
|
|
{ // dans ce cas on a affaire à une erreur d'arrondi
|
|
y=0.;
|
|
}
|
|
else if (x > 1.)
|
|
{ // il y a un pb dans l'algo, on le signal mais on continu
|
|
cout << "\n erreur dans la fonction Qartic::acos3, l'argument est > 1 , on limite ! "
|
|
<< " x-1= " << (x-1.);
|
|
y=0.;
|
|
};
|
|
};
|
|
break;
|
|
};
|
|
#endif
|
|
// calcul de la valeur de retour
|
|
value = std::cos(y*inv3);
|
|
return(value);
|
|
}; /* acos3 */
|
|
|
|
|
|
/***************************************/
|
|
// find cube root of x.
|
|
//
|
|
// 30 Jan 1989 Don Herbison-Evans
|
|
//
|
|
// called by cubic .
|
|
/***************************************/
|
|
|
|
double Quartic::curoot(double x)
|
|
{
|
|
double value;
|
|
double absx;
|
|
int neg;
|
|
|
|
neg = 0;
|
|
absx = x;
|
|
if (x < doub0)
|
|
{
|
|
absx = -x;
|
|
neg = 1;
|
|
}
|
|
if (absx != doub0) value = exp( log(absx)*inv3 );
|
|
else value = doub0;
|
|
if (neg == 1) value = -value;
|
|
return(value);
|
|
}; /* curoot */
|
|
|
|
/****************************************************/
|
|
// solve the quadratic equation -
|
|
//
|
|
// x**2 + b*x + c = 0
|
|
//
|
|
// 14 Jan 2004 cut determinant in quadratic call
|
|
// 29 Nov 2003 improved
|
|
// 16 Jul 1981 Don Herbison-Evans
|
|
//
|
|
// called by cubic,quartic,chris,descartes,ferrari,neumark.
|
|
/****************************************************/
|
|
|
|
int Quartic::quadratic(double b,double c,double rts[4])
|
|
{
|
|
int nquad;
|
|
double dis,rtdis ;
|
|
|
|
dis = b*b - doub4*c;
|
|
rts[0] = doub0;
|
|
rts[1] = doub0;
|
|
if (b == doub0)
|
|
{
|
|
if (c == doub0)
|
|
{
|
|
nquad = 2;
|
|
++nqud[0];
|
|
}
|
|
else
|
|
{
|
|
if (c < doub0)
|
|
{
|
|
nquad = 2;
|
|
rts[0] = sqrt(-c);
|
|
rts[1] = -rts[0];
|
|
++nqud[1];
|
|
}
|
|
else
|
|
{
|
|
nquad = 0;
|
|
++nqud[2];
|
|
}
|
|
}
|
|
}
|
|
else
|
|
if (c == doub0)
|
|
{
|
|
nquad = 2;
|
|
rts[0] = -b;
|
|
++nqud[3];
|
|
}
|
|
else
|
|
if (dis >= doub0)
|
|
{
|
|
nquad = 2 ;
|
|
rtdis = sqrt(dis);
|
|
if (b > doub0)
|
|
{
|
|
rts[0] = ( -b - rtdis)*inv2;
|
|
++nqud[4];
|
|
}
|
|
else
|
|
{
|
|
rts[0] = ( -b + rtdis)*inv2;
|
|
++nqud[5];
|
|
}
|
|
if (rts[0] == doub0)
|
|
{
|
|
rts[1] = -b;
|
|
++nqud[6];
|
|
}
|
|
else
|
|
{
|
|
rts[1] = c/rts[0];
|
|
++nqud[7];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
nquad = 0;
|
|
++nqud[8];
|
|
}
|
|
if (debug < 1)
|
|
{
|
|
printf("quad b %g c %g dis %g\n",
|
|
b,c,dis);
|
|
printf(" %d %g %g\n",
|
|
nquad,rts[0],rts[1]);
|
|
}
|
|
return(nquad);
|
|
}; /* quadratic */
|
|
|
|
/**************************************************/
|
|
// find the real roots of the cubic -
|
|
// x**3 + p*x**2 + q*x + r = 0
|
|
//
|
|
// 12 Dec 2003 initialising n,m,po3
|
|
// 12 Dec 2003 allow return of 3 zero roots if p=q=r=0
|
|
// 2 Dec 2003 negating j if p>0
|
|
// 1 Dec 2003 changing v from (sinsqk > doub0) to (sinsqk >= doub0)
|
|
// 1 Dec 2003 test changing v from po3sq+po3sq to doub2*po3sq
|
|
// 16 Jul 1981 Don Herbison-Evans
|
|
//
|
|
// input parameters -
|
|
// p,q,r - coeffs of cubic equation.
|
|
//
|
|
// output-
|
|
// the number of real roots
|
|
// v3 - the roots.
|
|
//
|
|
// global constants -
|
|
// rt3 - sqrt(3)
|
|
// inv3 - 1/3
|
|
// doubmax - square root of largest number held by machine
|
|
//
|
|
// method -
|
|
// see D.E. Littlewood, "A University Algebra" pp.173 - 6
|
|
//
|
|
// 15 Nov 2003 output 3 real roots: Don Herbison-Evans
|
|
// Apr 1981 initial version: Charles Prineas
|
|
//
|
|
// called by cubictest,quartic,chris,yacfraid,neumark,descartes,ferrari.
|
|
// calls quadratic,acos3,curoot,cubnewton.
|
|
/**************************************************/
|
|
|
|
int Quartic::cubic(double p,double q,double r,double v3[4])
|
|
{
|
|
int j,n3;
|
|
double po3,po3sq,qo3,po3q;
|
|
double uo3,u2o3,uo3sq4,uo3cu4;
|
|
double v,vsq,wsq;
|
|
double m1,m2,mcube;
|
|
double muo3,s,scube,t,cosk,rt3sink,sinsqk;
|
|
|
|
m1=doub0; m2=doub0; po3=doub0;
|
|
v=doub0; uo3=doub0; cosk=doub0;
|
|
if (r == doub0)
|
|
{
|
|
++ncub[0];
|
|
n3 = quadratic(p,q,v3);
|
|
v3[n3++] = doub0;
|
|
goto done;
|
|
}
|
|
if ((p == doub0) && (q == doub0))
|
|
{
|
|
++ncub[1];
|
|
v3[0] = curoot(-r);
|
|
v3[1] = v3[0];
|
|
v3[2] = v3[0];
|
|
n3 = 3;
|
|
goto done;
|
|
}
|
|
n3 = 1;
|
|
if ((p > doubmax) || (p < -doubmax))
|
|
{
|
|
v3[0] = -p;
|
|
++ncub[2];
|
|
goto done;
|
|
}
|
|
if ((q > doubmax) || (q < -doubmax))
|
|
{
|
|
if (q > doub0)
|
|
{
|
|
v3[0] = -r/q;
|
|
++ncub[3];
|
|
goto done;
|
|
}
|
|
else
|
|
if (q < doub0)
|
|
{
|
|
v3[0] = -sqrt(-q);
|
|
++ncub[4];
|
|
goto done;
|
|
}
|
|
else
|
|
{
|
|
v3[0] = doub0;
|
|
++ncub[5];
|
|
goto done;
|
|
}
|
|
}
|
|
else
|
|
if ((r > doubmax)|| (r < -doubmax))
|
|
{
|
|
v3[0] = -curoot(r);
|
|
++ncub[6];
|
|
goto done;
|
|
}
|
|
else
|
|
{
|
|
po3 = p*inv3;
|
|
po3q = po3*q;
|
|
po3sq = po3*po3;
|
|
if (po3sq > doubmax)
|
|
{
|
|
v3[0] = -p;
|
|
++ncub[7];
|
|
goto done;
|
|
}
|
|
else
|
|
{
|
|
v = r + po3*(po3sq+po3sq - q);
|
|
if ((v > doubmax) || (v < -doubmax))
|
|
{
|
|
v3[0] = -p;
|
|
++ncub[8];
|
|
goto done;
|
|
}
|
|
else
|
|
{
|
|
vsq = v*v;
|
|
qo3 = q*inv3;
|
|
uo3 = qo3 - po3sq;
|
|
u2o3 = uo3 + uo3;
|
|
if ((u2o3 > doubmax) || (u2o3 < -doubmax))
|
|
{
|
|
if (p == doub0)
|
|
{
|
|
if (q > doub0)
|
|
{
|
|
v3[0] = -r/q;
|
|
++ncub[9];
|
|
goto done;
|
|
}
|
|
else
|
|
if (q < doub0)
|
|
{
|
|
v3[0] = -sqrt(-q);
|
|
++ncub[10];
|
|
goto done;
|
|
}
|
|
else
|
|
{
|
|
v3[0] = doub0;
|
|
++ncub[11];
|
|
goto done;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
v3[0] = -q/p;
|
|
++ncub[12];
|
|
goto done;
|
|
}
|
|
}
|
|
uo3sq4 = u2o3*u2o3;
|
|
if (uo3sq4 > doubmax)
|
|
{
|
|
if (p == doub0)
|
|
{
|
|
if (q > doub0)
|
|
{
|
|
v3[0] = -r/q;
|
|
++ncub[13];
|
|
goto done;
|
|
}
|
|
else
|
|
if (q < doub0)
|
|
{
|
|
v3[0] = -sqrt(-q);
|
|
++ncub[14];
|
|
goto done;
|
|
}
|
|
else
|
|
{
|
|
v3[0] = doub0;
|
|
++ncub[15];
|
|
goto done;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
v3[0] = -q/p;
|
|
++ncub[16];
|
|
goto done;
|
|
}
|
|
}
|
|
uo3cu4 = uo3sq4*uo3;
|
|
wsq = uo3cu4 + vsq;
|
|
if (wsq > doub0)
|
|
{
|
|
/*
|
|
cubic has one real root -
|
|
*/
|
|
if (v <= doub0)
|
|
{
|
|
mcube = ( -v + sqrt(wsq))*inv2;
|
|
++ncub[17];
|
|
}
|
|
else
|
|
{
|
|
mcube = ( -v - sqrt(wsq))*inv2;
|
|
++ncub[18];
|
|
}
|
|
m1 = curoot(mcube);
|
|
if (m1 != doub0)
|
|
{
|
|
m2 = -uo3/m1;
|
|
++ncub[19];
|
|
}
|
|
else
|
|
{
|
|
m2 = doub0;
|
|
++ncub[20];
|
|
}
|
|
v3[0] = m1 + m2 - po3;
|
|
}
|
|
else
|
|
{
|
|
/*
|
|
cubic has three real roots -
|
|
*/
|
|
if (uo3 < doub0)
|
|
{
|
|
muo3 = -uo3;
|
|
if (muo3 > doub0)
|
|
{
|
|
s = sqrt(muo3);
|
|
++ncub[21];
|
|
if (p > doub0)
|
|
{
|
|
s = -s;
|
|
++ncub[22];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
s = doub0;
|
|
++ncub[23];
|
|
}
|
|
scube = s*muo3;
|
|
if (scube == doub0)
|
|
{
|
|
v3[0] = m1 + m2 - po3;
|
|
n3 = 1;
|
|
++ncub[24];
|
|
}
|
|
else
|
|
{
|
|
t = -v/(scube+scube);
|
|
cosk = acos3(t);
|
|
v3[0] = (s+s)*cosk - po3;
|
|
n3 = 1 ;
|
|
sinsqk = doub1 - cosk*cosk;
|
|
if (sinsqk >= doub0)
|
|
{
|
|
rt3sink = rt3*sqrt(sinsqk);
|
|
v3[1] = s*(-cosk + rt3sink) - po3;
|
|
v3[2] = s*(-cosk - rt3sink) - po3;
|
|
n3 = 3;
|
|
++ncub[25];
|
|
}
|
|
else ++ncub[26];
|
|
}
|
|
}
|
|
else
|
|
/*
|
|
cubic has multiple root -
|
|
*/
|
|
{
|
|
++ncub[27];
|
|
v3[0] = curoot(v) - po3;
|
|
v3[1] = v3[0];
|
|
v3[2] = v3[0];
|
|
n3 = 3;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
done:
|
|
if (debug < 1)
|
|
{
|
|
printf("cubic %d %g %g %g\n",n3,p,q,r);
|
|
for (j = 0; j < n3; ++j)
|
|
printf(" %d %13g %13g\n",j,v3[j],r+v3[j]*(q+v3[j]*(p+v3[j])));
|
|
printf("v %g, uo3 %g, m1 %g, m2 %g, po3 %g, cosk %g\n",
|
|
v,uo3,m1,m2,po3,cosk);
|
|
for (j = 0; j < 28; ++j)
|
|
{
|
|
printf(" %d",ncub[j]);
|
|
if ((j%10) == 9) printf("\n");
|
|
}
|
|
printf("\n");
|
|
}
|
|
if (iterate == true) cubnewton(p,q,r,n3,v3);
|
|
return(n3) ;
|
|
}; /* cubic */
|
|
|
|
|
|
|
|
/****************************************************/
|
|
// improve roots of cubic by Newton-Raphson iteration
|
|
//
|
|
// 5 Jan 2004 Don Herbison-Evans
|
|
//
|
|
// called by cubic.
|
|
/****************************************************/
|
|
|
|
void Quartic::cubnewton(double p,double q,double r,int n3,double v3[4])
|
|
{
|
|
int j,k;
|
|
double corr,deriv,err,root;
|
|
|
|
if (debug < 2) printf("cubnewtona %d %g\n",
|
|
n3,v3[0]);
|
|
for (j = 0; j < n3; ++j)
|
|
{
|
|
for (k = 0; k < 4; ++k)
|
|
{
|
|
root = v3[j];
|
|
err = ((root + p)*root + q)*root + r;
|
|
deriv = (doub3*root + doub2*p)*root + q;
|
|
if (deriv != doub0) corr = err/deriv; else corr = doub0;
|
|
v3[j] -= corr;
|
|
if (debug < 1) printf("cubnewtonb %d %d %g %g %g %g %g\n",
|
|
j,k,root,err,deriv,corr,v3[j]);
|
|
}
|
|
}
|
|
}; /* cubnewton */
|
|
|
|
/****************************************************/
|
|
// Solve quartic equation using either
|
|
// quadratic, Ferrari's or Neumark's algorithm.
|
|
//
|
|
// called by
|
|
// calls descartes, ferrari, neumark, yacfraid.
|
|
//
|
|
// 15 Dec 2003 added yacfraid
|
|
// 10 Dec 2003 added descartes with neg coeffs
|
|
// 21 Jan 1989 Don Herbison-Evans
|
|
/****************************************************/
|
|
int Quartic::quartic(double a,double b,double c,double d,double rts[4])
|
|
{
|
|
int j,k,nq,nr;
|
|
double roots[4];
|
|
|
|
if (Dabs(a) > doubmax)
|
|
nr = yacfraid(a,b,c,d,rts);
|
|
else
|
|
if ((a == doub0) && (c == doub0))
|
|
{
|
|
nq = quadratic(b,d,roots);
|
|
nr = 0;
|
|
for (j = 0; j < nq; ++j)
|
|
{
|
|
if (roots[0] >= doub0)
|
|
{
|
|
rts[0] = sqrt(roots[0]);
|
|
rts[1] = -rts[0];
|
|
nr = 2;
|
|
}
|
|
if (roots[1] >= doub0)
|
|
{
|
|
rts[nr] = sqrt(roots[1]);
|
|
rts[nr+1] = -rts[nr];
|
|
nr += 2;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
k = 0;
|
|
if (a < doub0) k += 2;
|
|
if (b < doub0) k += 1;
|
|
if (c < doub0) k += 8;
|
|
if (d < doub0) k += 4;
|
|
switch (k)
|
|
{
|
|
case 0 : nr = neumark(a,b,c,d,rts); break;
|
|
case 1 : nr = neumark(a,b,c,d,rts); break;
|
|
case 2 : nr = neumark(a,b,c,d,rts); break;
|
|
case 3 : nr = ferrari(a,b,c,d,rts); break;
|
|
case 4 : nr = neumark(a,b,c,d,rts); break;
|
|
case 5 : nr = descartes(a,b,c,d,rts); break;
|
|
case 6 : nr = neumark(a,b,c,d,rts); break;
|
|
case 7 : nr = neumark(a,b,c,d,rts); break;
|
|
case 8 : nr = neumark(a,b,c,d,rts); break;
|
|
case 9 : nr = ferrari(a,b,c,d,rts); break;
|
|
case 10 : nr = neumark(a,b,c,d,rts); break;
|
|
case 11 : nr = neumark(a,b,c,d,rts); break;
|
|
case 12 : nr = neumark(a,b,c,d,rts); break;
|
|
case 13 : nr = neumark(a,b,c,d,rts); break;
|
|
case 14 : nr = neumark(a,b,c,d,rts); break;
|
|
case 15 : nr = descartes(-a,b,-c,d,rts); break;
|
|
}
|
|
if (k == 15)
|
|
for (j = 0; j < nr; ++j)
|
|
rts[j] = -rts[j];
|
|
}
|
|
return(nr);
|
|
}; /* quartic */
|
|
|
|
|
|
//*****************************************/
|
|
// Solve quartic equation using
|
|
// Descartes-Euler-Cardano algorithm
|
|
//
|
|
// called by quartic, compare, quartictest.
|
|
//
|
|
// Strong, T. "Elemementary and Higher Algebra"
|
|
// Pratt and Oakley, p. 469 (1859)
|
|
//
|
|
// 16 Jul 1981 Don Herbison-Evans
|
|
//*****************************************/
|
|
int Quartic::descartes(double a,double b,double c,double d,double rts[4])
|
|
{
|
|
int j;
|
|
double v1[4],v2[4],v3[4];
|
|
double k,y;
|
|
double p,q,r;
|
|
double e0,e1,e2;
|
|
double g,h;
|
|
double asq;
|
|
double ainv4;
|
|
double e1invk;
|
|
|
|
asq = a*a;
|
|
e2 = b - asq*d3o8;
|
|
e1 = c + a*(asq*inv8 - b*inv2);
|
|
e0 = d + asq*(b*inv16 - asq*d3o256) - a*c*inv4;
|
|
|
|
p = doub2*e2;
|
|
q = e2*e2 - doub4*e0;
|
|
r = -e1*e1;
|
|
|
|
n3 = cubic(p,q,r,v3);
|
|
for (j3 = 0; j3 < n3; ++j3)
|
|
{
|
|
y = v3[j3];
|
|
if (y <= doub0)
|
|
{
|
|
n4[j3] = 0;
|
|
++ndes[0];
|
|
} /* y<0 */
|
|
else
|
|
{
|
|
k = sqrt(y);
|
|
ainv4 = a*inv4;
|
|
e1invk = e1/k;
|
|
g = (y + e2 + e1invk)*inv2;
|
|
h = (y + e2 - e1invk)*inv2 ;
|
|
n1 = quadratic(-k, g, v1);
|
|
n2 = quadratic( k, h, v2);
|
|
qrts[0][j3] = v1[0] - ainv4;
|
|
qrts[1][j3] = v1[1] - ainv4;
|
|
qrts[n1][j3] = v2[0] - ainv4;
|
|
qrts[n1+1][j3] = v2[1] - ainv4;
|
|
n4[j3]= n1 + n2;
|
|
++ndes[1];
|
|
} /* y>=0 */
|
|
donej3:
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
worst3[j3] = errors(a,b,c,d,rts,rterd,n4[j3]);
|
|
} /* j3 loop */
|
|
done:
|
|
j3 = 0;
|
|
if (n3 != 1)
|
|
{
|
|
if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
|
|
++ndes[NCASES-2];
|
|
else
|
|
++ndes[NCASES-1];
|
|
if ((n4[1] > n4[j3]) ||
|
|
((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
|
|
if ((n4[2] > n4[j3]) ||
|
|
((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
|
|
}
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
if (debug < 1)
|
|
printf("descartes chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
|
|
++ndes[30+n4[j3]];
|
|
++ndes[35+j3];
|
|
return(n4[j3]);
|
|
}; /* descartes */
|
|
|
|
//****************************************************/
|
|
// solve the quartic equation -
|
|
//
|
|
// x**4 + a*x**3 + b*x**2 + c*x + d = 0
|
|
//
|
|
// calls cubic, quadratic.
|
|
//
|
|
// input -
|
|
// a,b,c,e - coeffs of equation.
|
|
//
|
|
// output -
|
|
// n4 - number of real roots.
|
|
// rts - array of root values.
|
|
//
|
|
// method : Ferrari - Lagrange
|
|
// Theory of Equations, H.W. Turnbull p. 140 (1947)
|
|
//
|
|
// 16 Jul 1981 Don Herbison-Evans
|
|
//
|
|
// calls cubic, quadratic
|
|
//****************************************************/
|
|
int Quartic::ferrari(double a,double b,double c,double d,double rts[4])
|
|
{
|
|
int j;
|
|
double asqinv4;
|
|
double ainv2;
|
|
double d4;
|
|
double yinv2;
|
|
double v1[4],v2[4],v3[4];
|
|
double p,q,r;
|
|
double y;
|
|
double e,f,esq,fsq,ef;
|
|
double g,gg,h,hh;
|
|
|
|
ainv2 = a*inv2;
|
|
asqinv4 = ainv2*ainv2;
|
|
d4 = d*doub4 ;
|
|
|
|
p = b;
|
|
q = a*c-d4;
|
|
r = (asqinv4 - b)*d4 + c*c;
|
|
n3 = cubic(p,q,r,v3);
|
|
for (j3 = 0; j3 < n3; ++j3)
|
|
{
|
|
y = v3[j3];
|
|
yinv2 = y*inv2;
|
|
esq = asqinv4 - b - y;
|
|
fsq = yinv2*yinv2 - d;
|
|
if ((esq < doub0) && (fsq < doub0))
|
|
{
|
|
n4[j3] = 0;
|
|
++nfer[0];
|
|
}
|
|
else
|
|
{
|
|
ef = -(inv4*a*y + inv2*c);
|
|
if ( ((a > doub0)&&(y > doub0)&&(c > doub0))
|
|
|| ((a > doub0)&&(y < doub0)&&(c < doub0))
|
|
|| ((a < doub0)&&(y > doub0)&&(c < doub0))
|
|
|| ((a < doub0)&&(y < doub0)&&(c > doub0))
|
|
|| (a == doub0)||(y == doub0)||(c == doub0))
|
|
/* use ef - */
|
|
{
|
|
if ((b < doub0)&&(y < doub0))
|
|
{
|
|
e = sqrt(esq);
|
|
f = ef/e;
|
|
++nfer[1];
|
|
}
|
|
else if (d < doub0)
|
|
{
|
|
f = sqrt(fsq);
|
|
e = ef/f;
|
|
++nfer[2];
|
|
}
|
|
else
|
|
{
|
|
if (esq > doub0)
|
|
{
|
|
e = sqrt(esq);
|
|
++nfer[3];
|
|
}
|
|
else
|
|
{
|
|
e = doub0;
|
|
++nfer[4];
|
|
}
|
|
if (fsq > doub0)
|
|
{
|
|
f = sqrt(fsq);
|
|
++nfer[5];
|
|
}
|
|
else
|
|
{
|
|
f = doub0;
|
|
++nfer[6];
|
|
}
|
|
if (ef < doub0)
|
|
{
|
|
f = -f;
|
|
++nfer[7];
|
|
}
|
|
}
|
|
}
|
|
else
|
|
/* use esq and fsq - */
|
|
{
|
|
if (esq > doub0)
|
|
{
|
|
e = sqrt(esq);
|
|
++nfer[8];
|
|
}
|
|
else
|
|
{
|
|
e = doub0;
|
|
++nfer[9];
|
|
}
|
|
if (fsq > doub0)
|
|
{
|
|
f = sqrt(fsq);
|
|
++nfer[10];
|
|
}
|
|
else
|
|
{
|
|
f = doub0;
|
|
++nfer[11];
|
|
}
|
|
if (ef < doub0)
|
|
{
|
|
f = -f;
|
|
++nfer[12];
|
|
}
|
|
}
|
|
/* note that e >= doub0 */
|
|
g = ainv2 - e;
|
|
gg = ainv2 + e;
|
|
if ( ((b > doub0)&&(y > doub0))
|
|
|| ((b < doub0)&&(y < doub0)) )
|
|
{
|
|
if ((a > doub0) && (e > doub0)
|
|
|| (a < doub0) && (e < doub0) )
|
|
{
|
|
g = (b + y)/gg;
|
|
++nfer[13];
|
|
}
|
|
else
|
|
if ((a > doub0) && (e < doub0)
|
|
|| (a < doub0) && (e > doub0) )
|
|
{
|
|
gg = (b + y)/g;
|
|
++nfer[14];
|
|
}
|
|
else
|
|
++nfer[15];
|
|
}
|
|
hh = -yinv2 + f;
|
|
h = -yinv2 - f;
|
|
if ( ((f > doub0)&&(y < doub0))
|
|
|| ((f < doub0)&&(y > doub0)) )
|
|
{
|
|
h = d/hh;
|
|
++nfer[16];
|
|
}
|
|
else
|
|
if ( ((f < doub0)&&(y < doub0))
|
|
|| ((f > doub0)&&(y > doub0)) )
|
|
{
|
|
hh = d/h;
|
|
++nfer[17];
|
|
}
|
|
else
|
|
++nfer[18];
|
|
|
|
n1 = quadratic(gg,hh,v1);
|
|
n2 = quadratic(g,h,v2);
|
|
n4[j3] = n1+n2;
|
|
qrts[0][j3] = v1[0];
|
|
qrts[1][j3] = v1[1];
|
|
qrts[n1+0][j3] = v2[0];
|
|
qrts[n1+1][j3] = v2[1];
|
|
}
|
|
donej3:
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
worst3[j3] = errors(a,b,c,d,rts,rterf,n4[j3]);
|
|
} /* j3 loop */
|
|
done:
|
|
j3 = 0;
|
|
if (n3 != 1)
|
|
{
|
|
if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
|
|
++nfer[NCASES-2];
|
|
else
|
|
++nfer[NCASES-1];
|
|
if ((n4[1] > n4[j3]) ||
|
|
((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
|
|
if ((n4[2] > n4[j3]) ||
|
|
((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
|
|
}
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
if (debug < 1)
|
|
printf("ferrari chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
|
|
++nfer[30+n4[j3]];
|
|
++nfer[35+j3];
|
|
return(n4[j3]);
|
|
}; /* ferrari */
|
|
|
|
|
|
//*****************************************/
|
|
// solve the quartic equation -
|
|
//
|
|
// x**4 + a*x**3 + b*x**2 + c*x + d = 0
|
|
//
|
|
// calls cubic, quadratic.
|
|
//
|
|
// input parameters -
|
|
// a,b,c,e - coeffs of equation.
|
|
//
|
|
// output parameters -
|
|
// n4 - number of real roots.
|
|
// rts - array of root values.
|
|
//
|
|
// method - S. Neumark
|
|
// "Solution of Cubic and Quartic Equations" - Pergamon 1965
|
|
//
|
|
// 1 Dec 1985 translated to C with help of Shawn Neely
|
|
// 16 Jul 1981 Don Herbison-Evans
|
|
//*****************************************/
|
|
int Quartic::neumark(double a,double b,double c,double d,double rts[4])
|
|
{
|
|
int j;
|
|
double y,g,gg,h,hh,gdis,gdisrt,hdis,hdisrt,g1,g2,h1,h2;
|
|
double bmy,gerr,herr,y4,bmysq;
|
|
double v1[4],v2[4],v3[4];
|
|
double asq;
|
|
double d4;
|
|
double p,q,r;
|
|
double hmax,gmax;
|
|
|
|
if (d == doub0)
|
|
{
|
|
n3 = 0;
|
|
n4[0] = cubic(a,b,c,rts);
|
|
for (j = 0; j < n4[0]; ++j)
|
|
qrts[j][0] = rts[j];
|
|
qrts[n4[0]++][0] = doub0;
|
|
goto done;
|
|
}
|
|
asq = a*a;
|
|
d4 = d*doub4;
|
|
p = -b*doub2;
|
|
q = b*b + a*c - d4;
|
|
r = (c - a*b)*c + asq*d;
|
|
if (debug < 3)
|
|
printf("neumarka %g %g %g %g, %g %g %g\n",
|
|
a,b,c,d,p,q,r);
|
|
n3 = cubic(p,q,r,v3);
|
|
for (j3 = 0; j3 < n3; ++j3)
|
|
{
|
|
y = v3[j3];
|
|
bmy = b - y;
|
|
y4 = y*doub4;
|
|
bmysq = bmy*bmy;
|
|
gdis = asq - y4;
|
|
hdis = bmysq - d4;
|
|
if (debug < 3)
|
|
printf("neumarkb %g %g\n",
|
|
gdis,hdis);
|
|
if ((gdis < doub0) || (hdis < doub0))
|
|
{
|
|
n4[j3] = 0;
|
|
++nneu[0];
|
|
}
|
|
else
|
|
{
|
|
g1 = a*inv2;
|
|
h1 = bmy*inv2;
|
|
gerr = asq + y4;
|
|
herr = hdis;
|
|
if (d > doub0)
|
|
{
|
|
herr = bmysq + d4;
|
|
++nneu[1];
|
|
}
|
|
if ((y < doub0) || (herr*gdis > gerr*hdis))
|
|
{
|
|
gdisrt = sqrt(gdis);
|
|
g2 = gdisrt*inv2;
|
|
if (gdisrt != doub0)
|
|
{
|
|
h2 = (a*h1 - c)/gdisrt;
|
|
++nneu[2];
|
|
}
|
|
else
|
|
{
|
|
h2 = doub0;
|
|
++nneu[3];
|
|
}
|
|
}
|
|
else
|
|
{
|
|
hdisrt = sqrt(hdis);
|
|
h2 = hdisrt*inv2;
|
|
if (hdisrt != doub0)
|
|
{
|
|
g2 = (a*h1 - c)/hdisrt;
|
|
++nneu[4];
|
|
}
|
|
else
|
|
{
|
|
g2 = doub0;
|
|
++nneu[5];
|
|
}
|
|
}
|
|
/*
|
|
note that in the following, the tests ensure non-zero
|
|
denominators -
|
|
*/
|
|
h = h1 - h2;
|
|
hh = h1 + h2;
|
|
hmax = hh;
|
|
if (hmax < doub0)
|
|
{
|
|
hmax = -hmax;
|
|
++nneu[6];
|
|
}
|
|
if (hmax < h)
|
|
{
|
|
hmax = h;
|
|
++nneu[7];
|
|
}
|
|
if (hmax < -h)
|
|
{
|
|
hmax = -h;
|
|
++nneu[8];
|
|
}
|
|
if ((h1 > doub0)&&(h2 > doub0))
|
|
{
|
|
h = d/hh;
|
|
++nneu[9];
|
|
}
|
|
if ((h1 < doub0)&&(h2 < doub0))
|
|
{
|
|
h = d/hh;
|
|
++nneu[10];
|
|
}
|
|
if ((h1 > doub0)&&(h2 < doub0))
|
|
{
|
|
hh = d/h;
|
|
++nneu[11];
|
|
}
|
|
if ((h1 < doub0)&&(h2 > doub0))
|
|
{
|
|
hh = d/h;
|
|
++nneu[12];
|
|
}
|
|
if (h > hmax)
|
|
{
|
|
h = hmax;
|
|
++nneu[13];
|
|
}
|
|
if (h < -hmax)
|
|
{
|
|
h = -hmax;
|
|
++nneu[14];
|
|
}
|
|
if (hh > hmax)
|
|
{
|
|
hh = hmax;
|
|
++nneu[15];
|
|
}
|
|
if (hh < -hmax)
|
|
{
|
|
hh = -hmax;
|
|
++nneu[16];
|
|
}
|
|
|
|
g = g1 - g2;
|
|
gg = g1 + g2;
|
|
gmax = gg;
|
|
if (gmax < doub0)
|
|
{
|
|
gmax = -gmax;
|
|
++nneu[17];
|
|
}
|
|
if (gmax < g)
|
|
{
|
|
gmax = g;
|
|
++nneu[18];
|
|
}
|
|
if (gmax < -g)
|
|
{
|
|
gmax = -g;
|
|
++nneu[19];
|
|
}
|
|
if ((g1 > doub0)&&(g2 > doub0))
|
|
{
|
|
g = y/gg;
|
|
++nneu[20];
|
|
}
|
|
if ((g1 < doub0)&&(g2 < doub0))
|
|
{
|
|
g = y/gg;
|
|
++nneu[21];
|
|
}
|
|
if ((g1 > doub0)&&(g2 < doub0))
|
|
{
|
|
gg = y/g;
|
|
++nneu[22];
|
|
}
|
|
if ((g1 < doub0)&&(g2 > doub0))
|
|
{
|
|
gg = y/g;
|
|
++nneu[23];
|
|
}
|
|
if (g > gmax)
|
|
{
|
|
g = gmax;
|
|
++nneu[24];
|
|
}
|
|
if (g < -gmax)
|
|
{
|
|
g = -gmax;
|
|
++nneu[25];
|
|
}
|
|
if (gg > gmax)
|
|
{
|
|
gg = gmax;
|
|
++nneu[26];
|
|
}
|
|
if (gg < -gmax)
|
|
{
|
|
gg = -gmax;
|
|
++nneu[27];
|
|
}
|
|
|
|
n1 = quadratic(gg,hh,v1);
|
|
n2 = quadratic(g,h,v2);
|
|
n4[j3] = n1 + n2;
|
|
qrts[0][j3] = v1[0];
|
|
qrts[1][j3] = v1[1];
|
|
qrts[n1+0][j3] = v2[0];
|
|
qrts[n1+1][j3] = v2[1];
|
|
}
|
|
donej3:
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
worst3[j3] = errors(a,b,c,d,rts,rtern,n4[j3]);
|
|
} /* j3 loop */
|
|
done:
|
|
j3 = 0;
|
|
if (n3 > 1)
|
|
{
|
|
if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
|
|
++nneu[NCASES-2];
|
|
else
|
|
++nneu[NCASES-1];
|
|
if ((n4[1] > n4[j3]) ||
|
|
((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
|
|
if ((n4[2] > n4[j3]) ||
|
|
((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
|
|
}
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
if (debug < 1)
|
|
printf("neumark chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
|
|
++nneu[30+n4[j3]];
|
|
++nneu[35+j3];
|
|
return(n4[j3]);
|
|
}; /* neumark */
|
|
|
|
|
|
|
|
//****************************************************/
|
|
// solve the quartic equation -
|
|
// x**4 + a*x**3 + b*x**2 + c*x + d = 0
|
|
//
|
|
// calls cubic, quadratic.
|
|
//
|
|
// input parameters -
|
|
// a,b,c,e - coeffs of equation.
|
|
//
|
|
// output parameters -
|
|
// n4 - number of real roots.
|
|
// rts - array of root values.
|
|
//
|
|
// method -
|
|
// K.S. Brown
|
|
// Reducing Quartics to Cubics,
|
|
// http://www.seanet.com/~ksbrown/kmath296.htm (1967)
|
|
//
|
|
// Michael Daoud Yacoub & Gustavo Fraidenraich
|
|
// "A new simple solution of the general quartic equation"
|
|
// Revised 16 Feb 2004
|
|
//
|
|
// 14 Nov 2003 Don Herbison-Evans
|
|
//*****************************************/
|
|
int Quartic::yacfraid(double a,double b,double c,double d,double rts[4])
|
|
{
|
|
int j;
|
|
double y;
|
|
double v3[4];
|
|
double asq,acu;
|
|
double b4;
|
|
double det0,det1,det2,det3;
|
|
double det0rt,det1rt,det2rt,det3rt;
|
|
double e,f,g,h,k;
|
|
double fsq,gsq,hsq,invk;
|
|
double P,Q,R,U;
|
|
|
|
if (d == doub0)
|
|
{
|
|
n3 = 0;
|
|
n4[0] = cubic(a,b,c,rts);
|
|
for (j = 0; j < n4[0]; ++j)
|
|
qrts[j][0] = rts[j];
|
|
qrts[n4[0]++][0] = doub0;
|
|
goto done;
|
|
}
|
|
asq = a*a;
|
|
acu = a*asq;
|
|
b4 = b*doub4;
|
|
n3 = 0;
|
|
|
|
P = asq*b - b4*b + doub2*a*c + doub16*d ;
|
|
Q = asq*c - b4*c + doub8*a*d;
|
|
R = asq*d - c*c ;
|
|
U = acu - b4*a + doub8*c;
|
|
n4[0] = 0;
|
|
if (U == doub0)
|
|
{
|
|
if (P == doub0)
|
|
{
|
|
det0 = doub3*asq - doub8*b;
|
|
if (det0 < doub0)
|
|
{
|
|
++nyac[0];
|
|
goto done;
|
|
}
|
|
det0rt = sqrt(det0);
|
|
qrts[0][0] = (-a + det0rt)*inv4;
|
|
qrts[1][0] = qrts[0][0];
|
|
qrts[2][0] = (-a - det0rt)*inv4;
|
|
qrts[3][0] = qrts[2][0];
|
|
++nyac[1];
|
|
n4[0] = 4;
|
|
goto done;
|
|
} /* P=0 */
|
|
else
|
|
{
|
|
det1 = asq*asq - doub8*asq*b + doub16*b*b - doub64*d;
|
|
if (det1 < doub0)
|
|
{
|
|
++nyac[2];
|
|
goto done;;
|
|
}
|
|
n4[0] = 0;
|
|
det1rt = sqrt(det1);
|
|
det2 = doub3*asq - doub8*b + doub2*det1rt;
|
|
if (det2 >= doub0)
|
|
{
|
|
det2rt = sqrt(det2);
|
|
qrts[0][0] = (-a + det2rt)*inv4;
|
|
qrts[1][0] = (-a - det2rt)*inv4;
|
|
n4[0] = 2;
|
|
++nyac[3];
|
|
}
|
|
det3 = doub3*asq - doub8*b - doub2*det1rt;
|
|
if (det3 >= doub0)
|
|
{
|
|
det3rt = sqrt(det3);
|
|
qrts[n4[0]++][0] = (-a + det3rt)*inv4;
|
|
qrts[n4[0]++][0] = (-a - det3rt)*inv4;
|
|
++nyac[5];
|
|
}
|
|
if (n4[0] == 0) ++nyac[6];
|
|
goto done;
|
|
} /* P<>0 */
|
|
}
|
|
|
|
n3 = cubic(P/U,Q/U,R/U,v3);
|
|
for (j3 = 0; j3 < n3; ++j3)
|
|
{
|
|
y = v3[j3];
|
|
j = 0;
|
|
k = a + doub4*y;
|
|
if (k == doub0)
|
|
{
|
|
++nyac[9];
|
|
goto donej3;
|
|
}
|
|
invk = doub1/k;
|
|
e = (acu - doub4*c - doub2*a*b + (doub6*asq - doub16*b)*y)*invk;
|
|
fsq = (acu + doub8*c - doub4*a*b)*invk;
|
|
if (fsq < doub0)
|
|
{
|
|
++nyac[10];
|
|
goto donej3;
|
|
}
|
|
f = sqrt(fsq);
|
|
gsq = doub2*(e + f*k);
|
|
hsq = doub2*(e - f*k);
|
|
if (gsq >= doub0)
|
|
{
|
|
++nyac[11];
|
|
g = sqrt(gsq);
|
|
qrts[j++][j3] = (-a - f - g)*inv4;
|
|
qrts[j++][j3] = (-a - f + g)*inv4;
|
|
}
|
|
if (hsq >= doub0)
|
|
{
|
|
++nyac[12];
|
|
h = sqrt(hsq);
|
|
qrts[j++][j3] = (-a + f - h)*inv4;
|
|
qrts[j++][j3] = (-a + f + h)*inv4;
|
|
}
|
|
if (debug < 1)
|
|
{
|
|
printf("j3 %d y %g k %g fsq %g gsq %g hsq %g\n",
|
|
j3,y,k,fsq,gsq,hsq);
|
|
printf("e %g f %g g %g h %g\n",e,f,g,h);
|
|
}
|
|
donej3:
|
|
n4[j3] = j;
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
worst3[j3] = errors(a,b,c,d,rts,rtery,n4[j3]);
|
|
} /* j3 loop */
|
|
done:
|
|
j3 = 0;
|
|
if (n3 > 1)
|
|
{
|
|
if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
|
|
++nyac[NCASES-2];
|
|
else
|
|
{
|
|
++nyac[NCASES-1];
|
|
if ((n4[0] != n4[1]) && (n4[0] != n4[2]) && (n4[1] != n4[2]))
|
|
printf("yace %d %d %d %g %g %g %g\n",
|
|
n4[0],n4[1],n4[2],a,b,c,d);
|
|
}
|
|
if ((n4[1] > n4[j3]) ||
|
|
((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
|
|
if ((n4[2] > n4[j3]) ||
|
|
((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
|
|
}
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
++nyac[30+n4[j3]];
|
|
++nyac[35+j3];
|
|
if (debug < 1)
|
|
printf("yacfraid chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
|
|
return(n4[j3]);
|
|
}; /* yacfraid */
|
|
|
|
|
|
//*****************************************/
|
|
// Solve quartic equation using
|
|
// Christianson's algorithm
|
|
// calls errors, quadratic, cubic.
|
|
//
|
|
// Bruce Christianson, Solving Quartics Using Palindromes,
|
|
// Mathematical Gazette, Vol. 75, pp. 327-328 (1991)
|
|
//
|
|
// 14 Jan 2004 Don Herbison-Evans
|
|
//*****************************************/
|
|
|
|
int Quartic::chris(double a,double b,double c,double d,double rts[4])
|
|
{
|
|
int j;
|
|
int n0,n1,n2;
|
|
double v0[4],v1[4],v2[4],v3[4];
|
|
double y,ysq,ycu;
|
|
double p,q,r;
|
|
double asq,acu,aqu,ao4;
|
|
double e0,e1,e2;
|
|
double k,ksq,kcu,kqu,kquinv,k1,k2;
|
|
double g,h,g1,Z1,Z2;
|
|
|
|
if (d == doub0)
|
|
{
|
|
n3 = 0;
|
|
n4[0] = cubic(a,b,c,rts);
|
|
for (j = 0; j < n4[0]; ++j)
|
|
qrts[j][0] = rts[j];
|
|
qrts[n4[0]++][0] = doub0;
|
|
goto done;
|
|
}
|
|
asq = a*a;
|
|
acu = asq*a;
|
|
aqu = acu*a;
|
|
ao4 = a*inv4;
|
|
e2 = b - d3o8*asq;
|
|
e1 = c - inv2*b*a + inv8*acu;
|
|
e0 = d - inv4*c*a + inv16*b*asq - d3o256*acu*a;
|
|
if (debug < 1)
|
|
printf("chrisa e0 %g e1 %g e2 %g\n",e0,e1,e2);
|
|
if (e1 == doub0)
|
|
{
|
|
n3 = 0;
|
|
n4[0] = 0;
|
|
++nchr[0];
|
|
n3 = 3;
|
|
v3[2] = doub0;
|
|
v3[1] = doub0;
|
|
v3[0] = -inv8*(doub16*e0 - doub4*e2*e2);
|
|
}
|
|
else
|
|
{
|
|
p = (inv2*b*asq - inv2*b*b - inv2*c*a + doub2*d - doub3*inv32*aqu)/e1;
|
|
q = doub3*inv16*asq - inv2*b;
|
|
r = inv16*b*a - inv8*c - inv64*acu;
|
|
if (debug < 1) printf("chrisb %g %g %g\n",p,q,r);
|
|
++nchr[1];
|
|
n3 = cubic(p,q,r,v3);
|
|
}
|
|
for (j3 = 0; j3 < n3; ++j3)
|
|
{
|
|
y = v3[j3];
|
|
n4[j3] = 0;
|
|
ysq = y*y;
|
|
ycu = y*ysq;
|
|
if ((y == doub0) ||
|
|
(( y < doub0) && (e1 <= doub0) && (e0 >= doub0) && (e2 >= doub0)))
|
|
{
|
|
kqu = y*ycu + e2*ysq + e1*y + e0;
|
|
if (kqu <= doub0)
|
|
{
|
|
ksq = doub0;
|
|
++nchr[2];
|
|
goto donej3;
|
|
}
|
|
ksq = sqrt(kqu);
|
|
}
|
|
else
|
|
{
|
|
ksq = ysq + inv2*e2 + inv4*e1/y;
|
|
if (ksq <= doub0)
|
|
{
|
|
kqu = doub0;
|
|
++nchr[3];
|
|
goto donej3;
|
|
}
|
|
kqu = ksq*ksq;
|
|
}
|
|
k = sqrt(ksq);
|
|
kcu = k*ksq;
|
|
kquinv = doub1/kqu;
|
|
g = doub4*y*kcu;
|
|
h = (doub6*ysq + e2)*ksq;
|
|
if (debug < 1)
|
|
{
|
|
k2 = sqrt(sqrt(Dabs(y*ycu + e2*ysq + e1*y + e0)));
|
|
if (y == doub0)
|
|
k1 = doub0;
|
|
else
|
|
k1 = sqrt(Dabs(ysq + inv2*e2 + inv4*e1/y));
|
|
g1 = (doub4*ycu + doub2* y*e2 + e1)*k;
|
|
printf("chrisc %g %g %g %g %g\n",k,k1,y,g,g1,h);
|
|
}
|
|
n0 = quadratic(g*kquinv,h*kquinv-doub2,v0);
|
|
if (n0 < 1)
|
|
{
|
|
++nchr[4];
|
|
goto donej3;
|
|
}
|
|
Z1 = v0[0];
|
|
Z2 = v0[1];
|
|
n1 = quadratic(-Z1,doub1,v1);
|
|
if (n1 > 0)
|
|
{
|
|
++nchr[5];
|
|
n4[j3] = n1;
|
|
qrts[0][j3] = y + k*v1[0] - ao4;
|
|
qrts[1][j3] = y + k*v1[1] - ao4;
|
|
}
|
|
n2 = quadratic(-Z2,doub1,v2);
|
|
if (n2 > 0)
|
|
{
|
|
++nchr[6];
|
|
n4[j3] += n2;
|
|
qrts[n1][j3] = y + k*v2[0] - ao4;
|
|
qrts[n1+1][j3] = y + k*v2[1] - ao4;
|
|
}
|
|
donej3:
|
|
if (debug < 1)
|
|
printf("chrisd %g %g\n",ksq,kqu);
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
worst3[j3] = errors(a,b,c,d,rts,rterc,n4[j3]);
|
|
} /* j3 loop */
|
|
done:
|
|
j3 = 0;
|
|
if (n3 > 1)
|
|
{
|
|
if ((n4[0] == n4[1]) && (n4[1] == n4[2]))
|
|
++nchr[NCASES-2];
|
|
else
|
|
{
|
|
++nchr[NCASES-1];
|
|
if (debug < 1)
|
|
printf("chrise %d %d %d %g %g %g %g\n",
|
|
n4[0],n4[1],n4[2],a,b,c,d);
|
|
}
|
|
if ((n4[1] > n4[j3]) ||
|
|
((worst3[1] < worst3[j3] ) && (n4[1] == n4[j3]))) j3 = 1;
|
|
if ((n4[2] > n4[j3]) ||
|
|
((worst3[2] < worst3[j3] ) && (n4[2] == n4[j3]))) j3 = 2;
|
|
}
|
|
for (j = 0; j < n4[j3]; ++j)
|
|
rts[j] = qrts[j][j3];
|
|
if (debug < 1)
|
|
printf("chris chose cubic %d %g %g\n\n",j3,v3[j3],worst3[j3]);
|
|
++nchr[30+n4[j3]];
|
|
++nchr[35+j3];
|
|
return(n4[j3]);
|
|
}; /* chris */
|
|
|