#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 %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 */