Program Listing for File PolynomialRoots-3-Quartic.cc¶
↰ Return to documentation for file (PolynomialRoots-3-Quartic.cc)
/*--------------------------------------------------------------------------*\
| |
| Copyright (C) 2017 |
| |
| , __ , __ |
| /|/ \ /|/ \ |
| | __/ _ ,_ | __/ _ ,_ |
| | \|/ / | | | | \|/ / | | | |
| |(__/|__/ |_/ \_/|/|(__/|__/ |_/ \_/|/ |
| /| /| |
| \| \| |
| |
| Enrico Bertolazzi |
| Dipartimento di Ingegneria Industriale |
| Universita` degli Studi di Trento |
| email: enrico.bertolazzi@unitn.it |
| |
\*--------------------------------------------------------------------------*/
#ifdef __GNUC__
#pragma GCC diagnostic ignored "-Wunused-function"
#endif
#ifdef __clang__
#pragma clang diagnostic ignored "-Wglobal-constructors"
#pragma clang diagnostic ignored "-Wvla-extension"
#pragma clang diagnostic ignored "-Wvla"
#pragma clang diagnostic ignored "-Wunused-function"
#pragma clang diagnostic ignored "-Wdocumentation-unknown-command"
#endif
#include "PolynomialRoots.hh"
#include <cmath>
#include <iostream>
#include <algorithm>
#include <limits>
#ifndef DOXYGEN_SHOULD_SKIP_THIS
#define MAX_ITER_SAFETY 50
namespace PolynomialRoots {
using std::abs;
static valueType const machepsi = std::numeric_limits<valueType>::epsilon();
indexType
Quartic::getRealRoots( valueType r[] ) const {
indexType nr = 0;
if ( !cplx0() ) r[nr++] = r0;
if ( !cplx1() ) r[nr++] = r1;
if ( !cplx2() ) r[nr++] = r2;
if ( !cplx3() ) r[nr++] = r3;
return nr;
}
indexType
Quartic::getPositiveRoots( valueType r[] ) const {
indexType nr = 0;
if ( !cplx0() && r0 > 0 ) r[nr++] = r0;
if ( !cplx1() && r1 > 0 ) r[nr++] = r1;
if ( !cplx2() && r2 > 0 ) r[nr++] = r2;
if ( !cplx3() && r3 > 0 ) r[nr++] = r3;
return nr;
}
indexType
Quartic::getNegativeRoots( valueType r[] ) const {
indexType nr = 0;
if ( !cplx0() && r0 < 0 ) r[nr++] = r0;
if ( !cplx1() && r1 < 0 ) r[nr++] = r1;
if ( !cplx2() && r2 < 0 ) r[nr++] = r2;
if ( !cplx3() && r3 < 0 ) r[nr++] = r3;
return nr;
}
indexType
Quartic::getRootsInRange( valueType a, valueType b, valueType r[] ) const {
indexType nr = 0;
if ( !cplx0() && r0 >= a && r0 <= b ) r[nr++] = r0;
if ( !cplx1() && r1 >= a && r1 <= b ) r[nr++] = r1;
if ( !cplx2() && r2 >= a && r2 <= b ) r[nr++] = r2;
if ( !cplx3() && r3 >= a && r3 <= b ) r[nr++] = r3;
return nr;
}
indexType
Quartic::getRootsInOpenRange( valueType a, valueType b, valueType r[] ) const {
indexType nr = 0;
if ( !cplx0() && r0 > a && r0 < b ) r[nr++] = r0;
if ( !cplx1() && r1 > a && r1 < b ) r[nr++] = r1;
if ( !cplx2() && r2 > a && r2 < b ) r[nr++] = r2;
if ( !cplx3() && r3 > a && r3 < b ) r[nr++] = r3;
return nr;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
//
// x^4 + A x^3 + B x^2 + C x + D
static
inline
void
scaleQuarticMonicPolynomial(
valueType A,
valueType B,
valueType C,
valueType D,
valueType & AS,
valueType & BS,
valueType & CS,
valueType & DS,
indexType & i_case,
valueType & scale
) {
valueType a = std::abs(A);
valueType b = std::sqrt(abs(B));
valueType c = std::cbrt(abs(C));
valueType d = std::sqrt(sqrt(abs(D)));
if ( a < b ) {
if ( b < c ) {
if ( c < d ) i_case = 0; // a < b < c < d --> d MAX
else i_case = 1; // a < b < c and d <= c --> c MAX
} else {
if ( b < d ) i_case = 0; // a < b and c <= b < d --> d MAX
else i_case = 2; // a < b and c <= b and d <= b --> b MAX
}
} else {
if ( a < c ) {
if ( c < d ) i_case = 0; // b <= a < c < d --> d MAX
else i_case = 1; // b <= a < c and d <= c --> c MAX
} else {
if ( a < d ) i_case = 0; // b <= a and c <= a and a < d --> d MAX
else i_case = 3; // b <= a and c <= a and d <= a --> a MAX
}
}
switch ( i_case ) {
case 0:
scale = d;
AS = A/d;
BS = (B/d)/d;
CS = ((C/d)/d)/d;
DS = D > 0 ? 1 : -1;
break;
case 1:
scale = c;
AS = A/c;
BS = (B/c)/c;
CS = C > 0 ? 1 : -1;
DS = (((D/c)/c)/c)/c;
break;
case 2:
scale = b;
AS = A/b;
BS = B > 0 ? 1 : -1;
CS = ((C/b)/b)/b;
DS = (((D/b)/b)/b)/b;
break;
case 3:
scale = a;
AS = A > 0 ? 1 : -1;
BS = (B/a)/a;
CS = ((C/a)/a)/a;
DS = (((D/a)/a)/a)/a;
break;
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
static
inline
valueType
evalHexic(
valueType x,
valueType q3,
valueType q2,
valueType q1,
valueType q0
) {
valueType t1 = x + q3;
valueType t2 = x + t1;
valueType t3 = x + t2;
t1 = t1 * x + q2;
t2 = t2 * x + t1;
t1 = t1 * x + q1;
valueType Q = t1 * x + q0; // Q(x)
valueType dQ = t2 * x + t1; // Q'(x)
valueType ddQ = t3 * x + t2; // Q''(x) / 2
valueType dddQ = x + t3; // Q'''(x) / 6
return Q * dddQ * dddQ - dQ * ddQ * dddQ + dQ * dQ; // H(x), usually < 0
}
static
inline
void
evalHexic(
valueType x,
valueType q3,
valueType q2,
valueType q1,
valueType q0,
valueType & p,
valueType & dp
) {
valueType t1 = x + q3;
valueType t2 = x + t1;
valueType t3 = x + t2;
t1 = t1 * x + q2;
t2 = t2 * x + t1;
t1 = t1 * x + q1;
valueType Q = t1 * x + q0; // Q(x)
valueType dQ = t2 * x + t1; // Q'(x)
valueType ddQ = t3 * x + t2; // Q''(x) / 2
valueType dddQ = x + t3; // Q'''(x) / 6
p = Q * dddQ * dddQ - dQ * ddQ * dddQ + dQ * dQ; // H(x), usually < 0
dp = 2 * dddQ * (4 * Q - dQ * dddQ - ddQ * ddQ); // H'(x)
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
//
// a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = (x-r)*(a4*x^3+b2*x^2+b1*x+b0)
static
void
deflateQuarticPolynomial(
valueType a4,
valueType a3,
valueType a2,
valueType a1,
valueType a0,
valueType r,
valueType & b2,
valueType & b1,
valueType & b0
) {
indexType i_cross = 0;
valueType r2 = r*r;
valueType v_cross = std::abs(a0);
valueType v_cross1 = std::abs(a1*r);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 1; }
v_cross1 = std::abs(a2*r2);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 2; }
v_cross1 = std::abs(a3*r*r2);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 3; }
v_cross1 = std::abs(a4*r2*r2);
if ( v_cross1 > v_cross ) i_cross = 4;
switch ( i_cross ) {
case 0: b2 = a3+a4*r; b1 = a2+r*b2; b0 = a1+r*b1; break;
case 1: b2 = a3+a4*r; b1 = a2+r*b2; b0 = -a0/r; break;
case 2: b2 = a3+a4*r; b0 = -a0/r; b1 = (b0-a1)/r; break;
case 3:
case 4: b0 = -a0/r; b1 = (b0-a1)/r; b2 = (b1-a2)/r; break;
}
}
/*
|| _ _ _ ____ _ _ _
|| | \ | | _____ _| |_ ___ _ __ | __ )(_)___ ___ ___| |_(_) ___ _ __
|| | \| |/ _ \ \ /\ / / __/ _ \| '_ \| _ \| / __|/ _ \/ __| __| |/ _ \| '_ \
|| | |\ | __/\ V V /| || (_) | | | | |_) | \__ \ __/ (__| |_| | (_) | | | |
|| |_| \_|\___| \_/\_/ \__\___/|_| |_|____/|_|___/\___|\___|\__|_|\___/|_| |_|
*/
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// Translate to C from Polynomial234RootSolvers
static
indexType
zeroQuarticByNewtonBisection(
valueType a,
valueType b,
valueType c,
valueType d,
valueType & x
) {
valueType p, dp;
evalMonicQuartic( x, a, b, c, d, p, dp );
valueType t = p; // save p(x) for sign comparison
x -= p/dp; // 1st improved root
indexType iter = 1;
indexType oscillate = 0;
indexType nconverged = 0;
bool bisection = false;
bool converged = false;
valueType s(0), u(0); // to mute warning
while ( ! ( nconverged > 1 || bisection ) && iter < MAX_ITER_SAFETY ) {
++iter;
valueType ddp;
evalMonicQuartic( x, a, b, c, d, p, dp, ddp );
if ( p*t < 0 ) { // does Newton start oscillating ?
if ( p < 0 ) {
++oscillate; // increment oscillation counter
s = x; // save lower bisection bound
} else {
u = x; // save upper bisection bound
}
t = p; // save current p(x)
}
//dp = (p*dp)/(dp*dp-p*ddp); // double zero correction
dp = (p*dp)/(dp*dp-0.5*p*ddp); // Halley correction
//dp = p/dp; // Newton correction
x -= dp; // new Newton root
bisection = oscillate > 2; // activate bisection
converged = std::abs(dp) <= (1+std::abs(x)) * machepsi; // Newton convergence indicator
if ( converged ) ++nconverged; else nconverged = 0;
}
if ( bisection || !converged ) {
t = u - s; // initial bisection interval
while ( std::abs(t) > (1+std::abs(x)) * machepsi && iter < MAX_ITER_SAFETY ) { // bisection iterates
++iter;
p = evalMonicQuartic( x, a, b, c, d );
if ( p < 0 ) s = x;
else u = x; // keep bracket on root
t = (u-s)/2; // new bisection interval
x = s + t; // new bisection root
}
}
return iter;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// Translate to C from Polynomial234RootSolvers
static
indexType
zeroHexicByNewtonBisection(
valueType q3,
valueType q2,
valueType q1,
valueType q0,
valueType & x
) {
valueType p, dp;
evalHexic( x, q3, q2, q1, q0, p, dp );
valueType t = p; // save p(x) for sign comparison
x -= p/dp; // 1st improved root
indexType iter = 1;
indexType oscillate = 0;
bool bisection = false;
bool converged = false;
valueType s(0), u(0); // to mute warning
while ( ! (converged||bisection) ) {
++iter;
evalHexic( x, q3, q2, q1, q0, p, dp );
if ( p*t < 0 ) { // does Newton start oscillating ?
if ( p < 0 ) {
++oscillate; // increment oscillation counter
s = x; // save lower bisection bound
} else {
u = x; // save upper bisection bound
}
t = p; // save current p(x)
}
dp = p/dp; // Newton correction
x -= dp; // new Newton root
bisection = oscillate > 2; // activate bisection
converged = std::abs(dp) <= std::abs(x) * machepsi; // Newton convergence indicator
}
if ( bisection ) {
t = u - s; // initial bisection interval
while ( std::abs(t) > std::abs(x) * machepsi ) { // bisection iterates
++iter;
p = evalHexic( x, q3, q2, q1, q0 );
if ( p < 0 ) s = x;
else u = x; // keep bracket on root
t = (u-s)/2; // new bisection interval
x = s + t; // new bisection root
}
}
return iter;
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
/*
|| _ __ _ _
|| __| | ___ / _| | __ _| |_ ___
|| / _` |/ _ \ |_| |/ _` | __/ _ \
|| | (_| | __/ _| | (_| | || __/
|| \__,_|\___|_| |_|\__,_|\__\___|
*/
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
//
// a3*x^3 + a2*x^2 + a1*x + a0 = (x-r)*(a3*x^2+b1*x+b0)
/*\
* Calculate the zeros of the quartic a*z^4 + b*z^3 + c*z^2 + d*x + e.
\*/
void
Quartic::findRoots() {
valueType const & A = ABCDE[0];
valueType const & B = ABCDE[1];
valueType const & C = ABCDE[2];
valueType const & D = ABCDE[3];
valueType const & E = ABCDE[4];
iter = nreal = ncplx = 0;
// special cases
if ( isZero(A) ) {
Cubic csolve( B, C, D, E );
nreal = csolve.numRoots();
switch ( nreal ) {
case 3: r2 = csolve.real_root2();
case 2: r1 = csolve.real_root1();
case 1: r0 = csolve.real_root0();
}
if ( csolve.complexRoots() ) { ncplx = 2; nreal -= 2; }
return;
}
if ( isZero(E) ) {
Cubic csolve( A, B, C, D );
nreal = csolve.numRoots();
r0 = r1 = r2 = r3 = 0;
switch ( nreal ) {
case 3: r2 = csolve.real_root2();
case 2: r1 = csolve.real_root1();
case 1: r0 = csolve.real_root0();
}
++nreal;
if ( csolve.complexRoots() ) { ncplx = 2; nreal -= 2; }
if ( nreal == 4 ) { // caso regolare, 4 radici reali maintain order
if ( r3 < r2 ) std::swap(r2,r3);
if ( r2 < r1 ) std::swap(r1,r2);
if ( r1 < r0 ) std::swap(r0,r1);
}
return;
}
if ( isZero(B) && isZero(D) ) { // biquadratic case
// A x^4 + C x^2 + E
Quadratic qsolve( A, C, E ) ;
valueType x = qsolve.real_root0() ;
valueType y = qsolve.real_root1() ;
if ( qsolve.complexRoots() ) {
// complex conjugate pair biquadratic roots x +/- iy.
ncplx = 4;
x /= 2; // re
y /= 2; // im
valueType z = hypot(x,y);
y = std::sqrt(z - x);
x = std::sqrt(z + x);
r0 = -x;
r1 = y;
r2 = x;
r3 = y;
} else {
// real roots of quadratic are ordered x <= y
if ( x >= 0 ) { // y >= 0
nreal = 4;
x = std::sqrt(x); y = std::sqrt(y);
r0 = -y; r1 = -x; r2 = x; r3 = y;
} else if ( y >= 0 ) { // x < 0 && y >= 0
nreal = ncplx = 2;
x = std::sqrt(-x); y = std::sqrt(y);
r0 = 0; r1 = x; // (real,imaginary)
r2 = -y; r3 = y;
} else { // x < 0 && y < 0
ncplx = 4;
x = std::sqrt(-x); y = std::sqrt(-y);
r0 = 0; r1 = x; r2 = 0; r3 = y; // 2 x (real,imaginary)
}
}
return;
}
/*
.. The general case. Rescale quartic polynomial, such that largest absolute
.. coefficient is (exactly!) equal to 1. Honor the presence of a special
.. quartic case that might have been obtained during the rescaling process
.. (due to underflow in the coefficients).
*/
valueType const A3 = B/A;
valueType const A2 = C/A;
valueType const A1 = D/A;
valueType const A0 = E/A;
valueType q0, q1, q2, q3, scale;
indexType i_case;
scaleQuarticMonicPolynomial( A3, A2, A1, A0, q3, q2, q1, q0, i_case, scale );
// Q(x) = q0 + q1 * x + q2 * x^2 + q3 * x^3 + x^4
// Q'(x)/4 = q1/4 + (q2/2) * x + (3/4)*q3 * x^2 + x^3
/*
.. The general quartic case. Search for stationary points. Set the first
.. derivative polynomial (cubic) equal to zero and find its roots.
.. Check, if any minimum point of Q(x) is below zero, in which case we
.. must have real roots for Q(x). Hunt down only the real root, which
.. will potentially converge fastest during Newton iterates. The remaining
.. roots will be determined by deflation Q(x) -> cubic.
..
.. The best roots for the Newton iterations are the two on the opposite
.. ends, i.e. those closest to the +2 and -2. Which of these two roots
.. to take, depends on the location of the Q(x) minima x = s and x = u,
.. with s > u. There are three cases:
..
.. 1) both Q(s) and Q(u) < 0
.. ----------------------
..
.. The best root is the one that corresponds to the lowest of
.. these minima. If Q(s) is lowest -> start Newton from +2
.. downwards (or zero, if s < 0 and a0 > 0). If Q(u) is lowest
.. -> start Newton from -2 upwards (or zero, if u > 0 and a0 > 0).
..
.. 2) only Q(s) < 0
.. -------------
..
.. With both sides +2 and -2 possible as a Newton starting point,
.. we have to avoid the area in the Q(x) graph, where inflection
.. points are present. Solving Q''(x) = 0, leads to solutions
.. x = -a3/4 +/- discriminant, i.e. they are centered around -a3/4.
.. Since both inflection points must be either on the r.h.s or l.h.s.
.. from x = s, a simple test where s is in relation to -a3/4 allows
.. us to avoid the inflection point area.
..
.. 3) only Q(u) < 0
.. -------------
..
.. Same of what has been said under 2) but with x = u.
*/
valueType c2 = 3 * (q3/4);
valueType c1 = q2/2;
valueType c0 = q1/4;
Cubic qsolve( 1, c2, c1, c0 );
valueType u = qsolve.real_root0(); // root according to paper
valueType t = qsolve.real_root1();
valueType s = qsolve.real_root2();
bool must_refine_r3 = true;
if ( !qsolve.complexRoots() ) {
valueType Qs = evalMonicQuartic( s, q3, q2, q1, q0 );
valueType Qu = evalMonicQuartic( u, q3, q2, q1, q0 );
bool q0pos = q0 > 0;
nreal = 1;
if ( Qs < 0 && Qu < 0 ) {
if ( Qs < Qu ) {
r3 = 2;
if ( q0pos && s < 0 ) r3 = 0;
} else {
r3 = -2;
if ( q0pos && u > 0 ) r3 = 0;
}
} else if ( Qs < 0 ) {
if ( 4*s < -q3 ) {
r3 = -2;
if ( q0pos && s > 0 ) r3 = 0;
} else {
r3 = 2;
if ( q0pos && s < 0 ) r3 = 0;
}
} else if ( Qu < 0 ) {
if ( 4*u < -q3 ) {
r3 = -2;
if ( q0pos && u > 0 ) r3 = 0;
} else {
r3 = 2;
if ( q0pos && u < 0 ) r3 = 0;
}
} else {
// check for astrological combination when s or u are root of the quartic
/*
DO NOT WORK
valueType epsi = 10 * ( 1 + std::abs(q3) +
std::abs(q2) +
std::abs(q1) +
std::abs(q0) ) * machepsi;
if ( Qs <= epsi ) r3 = s;
else if ( Qu <= epsi ) r3 = u;
else nreal = 0;
*/
if ( isZero(Qs) ) {
must_refine_r3 = false; r3 = s;
} else if ( isZero(Qu) ) {
must_refine_r3 = false; r3 = u;
} else {
nreal = 0;
}
}
} else {
// one single real root (only 1 minimum)
valueType Qs = evalMonicQuartic( s, q3, q2, q1, q0 );
if ( Qs <= 0 ) {
valueType tmp = q0 >= 0 ? 0 : 2;
r3 = u > 0 ? -tmp : -2;
nreal = 1;
}
}
/*
.. Do all necessary Newton iterations. In case we have more than 2 oscillations,
.. exit the Newton iterations and switch to bisection. Note, that from the
.. definition of the Newton starting point, we always have Q(x) > 0 and Q'(x)
.. starts (-ve/+ve) for the (-2/+2) starting points and (increase/decrease) smoothly
.. and staying (< 0 / > 0). In practice, for extremely shallow Q(x) curves near the
.. root, the Newton procedure can overshoot slightly due to rounding errors when
.. approaching the root. The result are tiny oscillations around the root. If such
.. a situation happens, the Newton iterations are abandoned after 3 oscillations
.. and further location of the root is done using bisection starting with the
.. oscillation brackets.
*/
if ( nreal > 0 ) {
if ( must_refine_r3 )
iter += zeroQuarticByNewtonBisection( q3, q2, q1, q0, r3 );
r3 *= scale;
/*
.. Find remaining roots -> reduce to cubic. The reduction to a cubic polynomial
.. is done using composite deflation to minimize rounding errors. Also, while
.. the composite deflation analysis is done on the reduced quartic, the actual
.. deflation is being performed on the original quartic again to avoid enhanced
.. propagation of root errors.
*/
deflateQuarticPolynomial( A, B, C, D, E, r3, q2, q1, q0 );
Cubic csolve( A, q2, q1, q0 );
r0 = csolve.real_root0();
r1 = csolve.real_root1();
r2 = csolve.real_root2();
if ( csolve.complexRoots() ) {
nreal = ncplx = 2;
if ( r2 > r3 ) std::swap( r2, r3 );
} else {
nreal = 4;
if ( r2 > r3 ) std::swap( r2, r3 );
if ( r1 > r2 ) std::swap( r1, r2 );
if ( r0 > r1 ) std::swap( r0, r1 );
}
} else {
/*
.. If no real roots have been found by now, only complex roots are
.. possible. Find real parts of roots first, followed by imaginary
.. components.
*/
s = q3/2;
t = s * s - q2;
u = s * t + q1; // value of Q'(-q3/4) at stationary point -q3/4
bool notZero = (abs(u) >= machepsi); // H(-a3/4) is considered > 0 at stationary point
bool minimum;
if ( !isZero(q3) ) {
s = q1 / q3;
minimum = q0 > s * s; // H''(-q3/4) > 0 -> minimum
} else {
minimum = 4 * q0 > q2 * q2; // H''(-q3/4) > 0 -> minimum
}
bool iterate = notZero || (!notZero && minimum);
valueType a, b, c, d;
if ( iterate ) {
valueType x = q3 >= 0 ? 2 : -2; // initial root -> target = smaller mag root
iter += zeroHexicByNewtonBisection( q3, q2, q1, q0, x );
a = x*scale; // 1st real component -> a
b = -A3/2 - a; // 2nd real component -> b
x = 4 * a + A3; // Q'''(a)
valueType y = x + 2*A3;
y = y * a + 2*A2; // Q'(a)
y = y * a + A1;
y /= x;
if ( y < 0 ) y = 0; // ensure Q'(a) / Q'''(a) >= 0
x = 4 * b + A3; // Q'''(b)
valueType z = x + 2*A3;
z = z * b + 2*A2; // Q'(b)
z = z * b + A1;
z /= x;
if ( z < 0 ) z = 0; // ensure Q'(b) / Q'''(b) >= 0
c = a * a; // store a^2 for later
d = b * b; // store b^2 for later
valueType ss = c + y; // magnitude^2 of (a + iy) root
valueType tt = d + z; // magnitude^2 of (b + iz) root
if ( ss > tt ) { // minimize imaginary error
c = std::sqrt(y); // 1st imaginary component -> c
d = std::sqrt(A0 / ss - d); // 2nd imaginary component -> d
} else {
c = std::sqrt(A0 / tt - c); // 1st imaginary component -> c
d = std::sqrt(z); // 2nd imaginary component -> d
}
} else { // no bisection -> real components equal
a = -A3/4; // 1st real component -> a
b = a; // 2nd real component -> b = a
valueType x = a + A3;
x = x * a + A2; // Q(a)
x = x * a + A1;
x = x * a + A0;
valueType y = A2/2 - 3*(a*a); // Q''(a) / 2
valueType z = y * y - x;
z = z > 0 ? std::sqrt(z) : 0; // force discriminant to be >= 0
// square root of discriminant
y = y > 0 ? y + z : y - z; // larger magnitude root
x /= y; // smaller magnitude root
c = y < 0 ? 0 : std::sqrt(y); // ensure root of biquadratic > 0
d = x < 0 ? 0 : std::sqrt(x); // large magnitude imaginary component
}
ncplx = 4;
if (a > b) { r0 = a; r1 = c; r2 = b; r3 = d; }
else if (a < b) { r0 = b; r1 = d; r2 = a; r3 = c; }
else { r0 = a; r1 = c; r2 = a; r3 = d; }
}
}
void
Quartic::info( std::ostream & s ) const {
valueType const & A = ABCDE[0];
valueType const & B = ABCDE[1];
valueType const & C = ABCDE[2];
valueType const & D = ABCDE[3];
valueType const & E = ABCDE[4];
s << "\npoly a=" << A << " b=" << B << " c=" << C << " d=" << D << " e=" << E
<< "\nn. complex = " << ncplx
<< "\nn. real = " << nreal;
if ( ncplx > 0 ) {
s << "\nx0 = (" << r0 << "," << r1 << ')'
<< "\nx1 = (" << r0 << "," << -r1 << ')';
} else {
if ( nreal > 0 ) s << "\nx0 = " << r0;
if ( nreal > 1 ) s << "\nx1 = " << r1;
}
if ( ncplx > 2 ) {
s << "\nx2 = (" << r2 << "," << r3 << ')'
<< "\nx3 = (" << r2 << "," << -r3 << ')';
} else {
if ( nreal > 2 || (ncplx > 0 && nreal > 0) ) s << "\nx2 = " << r2;
if ( nreal > 3 || (ncplx > 0 && nreal > 1) ) s << "\nx3 = " << r3;
}
s << '\n';
}
bool
Quartic::check( std::ostream & s ) const {
valueType const & A = ABCDE[0];
valueType const & B = ABCDE[1];
valueType const & C = ABCDE[2];
valueType const & D = ABCDE[3];
valueType const & E = ABCDE[4];
bool ok = true;
valueType epsi = 1000 * ( ( std::abs(A) +
std::abs(B) +
std::abs(C) +
std::abs(D) +
std::abs(E) )*machepsi) ;
if ( ncplx > 0 ) {
valueType z0 = std::abs(eval( root0() ));
valueType z1 = std::abs(eval( root1() ));
s << "|p(r0)| = " << z0 << "\n|p(r1)| = " << z1 << '\n';
ok = ok && std::abs(z0) < epsi && std::abs(z1) < epsi;
} else {
if ( nreal > 0 ) {
valueType z0 = eval( real_root0() );
s << "p(r0) = " << z0 << '\n';
ok = ok && std::abs(z0) < epsi;
}
if ( nreal > 1 ) {
valueType z1 = eval( real_root1() );
s << "p(r1) = " << z1 << '\n';
ok = ok && std::abs(z1) < epsi;
}
}
if ( ncplx > 2 ) {
valueType z2 = std::abs(eval( root2() ));
valueType z3 = std::abs(eval( root3() ));
s << "|p(r2)| = " << z2 << "\n|p(r3)| = " << z3 << '\n';
ok = ok && std::abs(z2) < epsi && std::abs(z3) < epsi;
} else {
if ( nreal > 2 || (ncplx > 0 && nreal > 0) ) {
valueType z2 = eval( real_root2() );
s << "p(r2) = " << z2 << '\n';
ok = ok && std::abs(z2) < epsi;
}
if ( nreal > 3 || (ncplx > 0 && nreal > 1) ) {
valueType z3 = eval( real_root3() );
s << "p(r3) = " << z3 << '\n';
ok = ok && std::abs(z3) < epsi;
}
}
return ok;
}
}
#endif
// EOF: PolynomialRoots-3-Quartic.cc
