Skip to content

Commit

Permalink
added missing files in toolbox
Browse files Browse the repository at this point in the history
  • Loading branch information
@ebertolazzi
ebertolazzi committed Feb 19, 2021
1 parent 51f9580 commit 5d2d189
Show file tree
Hide file tree
Showing 5 changed files with 90 additions and 96 deletions.
6 changes: 3 additions & 3 deletions toolbox/src/PolynomialRoots-1-Quadratic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@

namespace PolynomialRoots {

using std::abs;
using std::abs;
static valueType const machepsi = std::numeric_limits<valueType>::epsilon();

indexType
Expand Down Expand Up @@ -102,8 +102,8 @@ namespace PolynomialRoots {
}
} else { // Compute discriminant avoiding overflow.
valueType hb = B/2; // b now b/2
valueType abs_b = std::abs(hb);
valueType abs_c = std::abs(C);
valueType abs_b = abs(hb);
valueType abs_c = abs(C);
valueType e, d;
if ( abs_b < abs_c ) {
e = C < 0 ? -A : A;
Expand Down
14 changes: 7 additions & 7 deletions toolbox/src/PolynomialRoots-2-Cubic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,11 +242,11 @@ namespace PolynomialRoots {
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
converged = std::abs(dp) <= std::abs(x) * machepsi; // Newton convergence indicator
}
if ( bisection ) {
t = u - s; // initial bisection interval
while ( std::abs(t) > (1+std::abs(x)) * machepsi ) { // bisection iterates
while ( abs(t) > abs(x) * machepsi ) { // bisection iterates
++iter;
p = evalMonicCubic( x, a, b, c );
if ( p < 0 ) s = x;
Expand Down Expand Up @@ -366,17 +366,17 @@ namespace PolynomialRoots {
case 5:
r0 = a[1]-third;
r1 = a[0]+one27th;
trpx = std::abs(r0) <= machepsi && std::abs(r1) <= machepsi; // check for triple root
trpx = abs(r0) <= machepsi && abs(r1) <= machepsi; // check for triple root
if ( trpx ) { r0 = r1 = r2 = third * scale; nrts = 3; return; }
use_shifted = std::abs(r0) <= 0.01 && std::abs(r1) <= 0.01;
use_shifted = abs(r0) <= 0.01 && abs(r1) <= 0.01;
r2 = guess5(a);
break;
case 6:
r0 = a[1]-third;
r1 = a[0]-one27th;
trpx = std::abs(r0) <= machepsi && std::abs(r1) <= machepsi; // check for triple root
trpx = abs(r0) <= machepsi && abs(r1) <= machepsi; // check for triple root
if ( trpx ) { r0 = r1 = r2 = -third * scale; nrts = 3; return; }
use_shifted = std::abs(r0) <= 0.01 && std::abs(r1) <= 0.01;
use_shifted = abs(r0) <= 0.01 && abs(r1) <= 0.01;
r2 = guess6(a);
break;
}
Expand All @@ -400,7 +400,7 @@ namespace PolynomialRoots {
// y^3 + (a[1]-1/3)* y + (a[0]-a[1]/3+2/27), x = y+1/3
r2 += third; // shift guess
r1 -= a[1]/3-one27th;
//if ( std::abs(r3) < machepsi ) r3 = 0;
//if ( abs(r3) < machepsi ) r3 = 0;
iter = NewtonBisection( 0, r0, a[0]-a[1]/3+two27th, r2 );
r2 -= third; // unshift solution
}
Expand Down
64 changes: 29 additions & 35 deletions toolbox/src/PolynomialRoots-3-Quartic.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,10 +77,10 @@ namespace PolynomialRoots {
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)));
valueType a = abs(A);
valueType b = sqrt(abs(B));
valueType c = cbrt(abs(C));
valueType d = sqrt(sqrt(abs(D)));

if ( a < b ) {
if ( b < c ) {
Expand Down Expand Up @@ -201,14 +201,14 @@ namespace PolynomialRoots {
) {
indexType i_cross = 0;
valueType r2 = r*r;
valueType v_cross = std::abs(a0);
valueType v_cross1 = std::abs(a1*r);
valueType v_cross = abs(a0);
valueType v_cross1 = abs(a1*r);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 1; }
v_cross1 = std::abs(a2*r2);
v_cross1 = abs(a2*r2);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 2; }
v_cross1 = std::abs(a3*r*r2);
v_cross1 = abs(a3*r*r2);
if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = 3; }
v_cross1 = std::abs(a4*r2*r2);
v_cross1 = 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;
Expand Down Expand Up @@ -267,11 +267,11 @@ namespace PolynomialRoots {
//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
converged = std::abs(dp) <= std::abs(x) * machepsi; // Newton convergence indicator
}
if ( bisection || !converged ) {
t = u - s; // initial bisection interval
while ( std::abs(t) > (1+std::abs(x)) * machepsi ) { // bisection iterates
while ( abs(t) > abs(x) * machepsi ) { // bisection iterates
++iter;
p = evalMonicQuartic( x, a, b, c, d );
if ( p < 0 ) s = x;
Expand Down Expand Up @@ -319,11 +319,11 @@ namespace PolynomialRoots {
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
converged = abs(dp) <= abs(x) * machepsi; // Newton convergence indicator
}
if ( bisection ) {
t = u - s; // initial bisection interval
while ( std::abs(t) > std::abs(x) * machepsi ) { // bisection iterates
while ( abs(t) > abs(x) * machepsi ) { // bisection iterates
++iter;
p = evalHexic( x, q3, q2, q1, q0 );
if ( p < 0 ) s = x;
Expand Down Expand Up @@ -404,8 +404,8 @@ namespace PolynomialRoots {
x /= 2; // re
y /= 2; // im
valueType z = hypot(x,y);
y = std::sqrt(z - x);
x = std::sqrt(z + x);
y = sqrt(z - x);
x = sqrt(z + x);
r0 = -x;
r1 = y;
r2 = x;
Expand All @@ -414,16 +414,16 @@ namespace PolynomialRoots {
// real roots of quadratic are ordered x <= y
if ( x >= 0 ) { // y >= 0
nreal = 4;
x = std::sqrt(x); y = std::sqrt(y);
x = sqrt(x); y = 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);
x = sqrt(-x); y = 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);
x = sqrt(-x); y = sqrt(-y);
r0 = 0; r1 = x; r2 = 0; r3 = y; // 2 x (real,imaginary)
}
}
Expand Down Expand Up @@ -495,7 +495,6 @@ namespace PolynomialRoots {
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 );
Expand All @@ -522,13 +521,9 @@ namespace PolynomialRoots {
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;
}
if ( isZero(Qs) ) r3 = s;
else if ( isZero(Qu) ) r3 = u;
else nreal = 0;
}
} else {
// one single real root (only 1 minimum)
Expand All @@ -553,8 +548,7 @@ namespace PolynomialRoots {
.. oscillation brackets.
*/
if ( nreal > 0 ) {
if ( must_refine_r3 )
iter += zeroQuarticByNewtonBisection( q3, q2, q1, q0, r3 );
iter += zeroQuarticByNewtonBisection( q3, q2, q1, q0, r3 );
r3 *= scale;

/*
Expand Down Expand Up @@ -627,11 +621,11 @@ namespace PolynomialRoots {
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
c = sqrt(y); // 1st imaginary component -> c
d = 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
c = sqrt(A0 / tt - c); // 1st imaginary component -> c
d = sqrt(z); // 2nd imaginary component -> d
}

} else { // no bisection -> real components equal
Expand All @@ -645,12 +639,12 @@ namespace PolynomialRoots {
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
z = z > 0 ? 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
c = y < 0 ? 0 : sqrt(y); // ensure root of biquadratic > 0
d = x < 0 ? 0 : sqrt(x); // large magnitude imaginary component
}

ncplx = 4;
Expand Down
Loading

0 comments on commit 5d2d189

Please sign in to comment.