diff --git a/toolbox/src/PolynomialRoots-1-Quadratic.cc b/toolbox/src/PolynomialRoots-1-Quadratic.cc index 49b4d50a..e4f8647c 100644 --- a/toolbox/src/PolynomialRoots-1-Quadratic.cc +++ b/toolbox/src/PolynomialRoots-1-Quadratic.cc @@ -25,7 +25,7 @@ namespace PolynomialRoots { - using std::abs; + using std::abs; static valueType const machepsi = std::numeric_limits::epsilon(); indexType @@ -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; diff --git a/toolbox/src/PolynomialRoots-2-Cubic.cc b/toolbox/src/PolynomialRoots-2-Cubic.cc index b62ac169..3e607ba1 100644 --- a/toolbox/src/PolynomialRoots-2-Cubic.cc +++ b/toolbox/src/PolynomialRoots-2-Cubic.cc @@ -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; @@ -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; } @@ -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 } diff --git a/toolbox/src/PolynomialRoots-3-Quartic.cc b/toolbox/src/PolynomialRoots-3-Quartic.cc index 49688cb4..55a1d2ad 100644 --- a/toolbox/src/PolynomialRoots-3-Quartic.cc +++ b/toolbox/src/PolynomialRoots-3-Quartic.cc @@ -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 ) { @@ -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; @@ -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; @@ -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; @@ -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; @@ -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) } } @@ -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 ); @@ -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) @@ -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; /* @@ -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 @@ -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; diff --git a/toolbox/src/PolynomialRoots-Jenkins-Traub.cc b/toolbox/src/PolynomialRoots-Jenkins-Traub.cc index d7a7b89d..ba03aeef 100644 --- a/toolbox/src/PolynomialRoots-Jenkins-Traub.cc +++ b/toolbox/src/PolynomialRoots-Jenkins-Traub.cc @@ -123,11 +123,11 @@ namespace PolynomialRoots { // Synthetic division of K by the quadratic 1, u, v QuadraticSyntheticDivision( N, u, v, K, qk, c, d ); - if ( std::abs(c) <= epsilon100 * std::abs(K[N-1]) && - std::abs(d) <= epsilon100 * std::abs(K[N-2]) ) return dumFlag; + if ( abs(c) <= epsilon100 * abs(K[N-1]) && + abs(d) <= epsilon100 * abs(K[N-2]) ) return dumFlag; h = v*b; - if ( std::abs(d) >= std::abs(c) ) { + if ( abs(d) >= abs(c) ) { dumFlag = 2; // TYPE = 2 indicates that all formulas are divided by d e = a / d; f = c / d; @@ -169,7 +169,7 @@ namespace PolynomialRoots { for ( indexType i = 2; i < N; ++i ) K[i] = qk[i-2]; } else { valueType temp = tFlag == 1 ? b : a; - if ( std::abs(a1) > epsilon10 * std::abs(temp) ) { + if ( abs(a1) > epsilon10 * abs(temp) ) { // Use scaled form of the recurrence a7 /= a1; a3 /= a1; @@ -265,11 +265,11 @@ namespace PolynomialRoots { qp[0] = pv; for ( indexType i = 1; i < NN; ++i ) qp[i] = pv = pv * s + p[i]; - valueType mp = std::abs(pv); + valueType mp = abs(pv); // Compute a rigorous bound on the error in evaluating p - valueType ms = std::abs(s); - valueType ee = 0.5 * std::abs(qp[0]); - for ( indexType i = 1; i < NN; ++i ) ee = ee * ms + std::abs(qp[i]); + valueType ms = abs(s); + valueType ee = 0.5 * abs(qp[0]); + for ( indexType i = 1; i < NN; ++i ) ee = ee * ms + abs(qp[i]); // Iteration has converged sufficiently if the polynomial // value is less than 20 times this bound if ( mp <= 20.0 * epsilon * (2*ee - mp) ) @@ -278,7 +278,7 @@ namespace PolynomialRoots { // Stop iteration after 10 steps if ( ++j > 10 ) break; if ( j >= 2 ) { - if ( (abs(t) <= 0.001 * std::abs(s-t)) && (mp > omp) ) { + if ( (abs(t) <= 0.001 * abs(s-t)) && (mp > omp) ) { // A cluster of zeros near the real axis has been encountered; // Return with iFlag set to initiate a quadratic iteration iFlag = 1; @@ -293,7 +293,7 @@ namespace PolynomialRoots { for ( indexType i = 1; i < N; ++i ) qk[i] = kv = kv * s + K[i]; - if ( std::abs(kv) > std::abs(K[N-1]) * epsilon10 ) { + if ( abs(kv) > abs(K[N-1]) * epsilon10 ) { // Use the scaled form of the recurrence if the value of K at s is non - zero valueType tt = -(pv / kv); K[0] = qp[0]; @@ -304,7 +304,7 @@ namespace PolynomialRoots { } kv = K[0]; for ( indexType i = 1; i < N; ++i ) kv = kv * s + K[i]; - t = std::abs(kv) > std::abs(K[N-1])*epsilon10 ? -(pv/kv) : 0; + t = abs(kv) > abs(K[N-1])*epsilon10 ? -(pv/kv) : 0; s += t; } } @@ -356,16 +356,16 @@ namespace PolynomialRoots { // Return if roots of the quadratic are real and not close // to multiple or nearly equal and of opposite sign. - if ( std::abs(abs(szr) - std::abs(lzr)) > 0.01 * std::abs(lzr) ) break; + if ( abs(abs(szr) - abs(lzr)) > 0.01 * abs(lzr) ) break; // Evaluate polynomial by quadratic synthetic division QuadraticSyntheticDivision( NN, u, v, p, qp, a, b ); - valueType mp = std::abs(a-(szr*b)) + std::abs(szi*b); + valueType mp = abs(a-(szr*b)) + abs(szi*b); // Compute a rigorous bound on the rounding error in evaluating p - valueType zm = std::sqrt(abs(v)); - valueType ee = 2 * std::abs(qp[0]); + valueType zm = sqrt(abs(v)); + valueType ee = 2 * abs(qp[0]); valueType t = -szr*b; - for ( indexType i = 1; i < N; ++i ) ee = ee * zm + std::abs(qp[i]); - ee = ee * zm + std::abs(a+t); + for ( indexType i = 1; i < N; ++i ) ee = ee * zm + abs(qp[i]); + ee = ee * zm + abs(a+t); ee = (9*ee+2*abs(t)-7*(abs(a+t)+zm*abs(b)))*epsilon; // Iteration has converged sufficiently if the polynomial // value is less than 20 times this bound @@ -376,7 +376,7 @@ namespace PolynomialRoots { if ( (relstp <= 0.01) && (mp >= omp) && !triedFlag ) { // A cluster appears to be stalling the convergence. // Five fixed shift steps are taken with a u, v close to the cluster. - relstp = (relstp < epsilon) ? std::sqrt(epsilon) : std::sqrt(relstp); + relstp = (relstp < epsilon) ? sqrt(epsilon) : sqrt(relstp); u -= u * relstp; v += v * relstp; QuadraticSyntheticDivision(NN, u, v, p, qp, a, b); @@ -397,7 +397,7 @@ namespace PolynomialRoots { // If vi is zero, the iteration is not converging if ( !isZero(vi) ) { - relstp = std::abs((vi-v)/vi); + relstp = abs((vi-v)/vi); u = ui; v = vi; } @@ -467,8 +467,8 @@ namespace PolynomialRoots { valueType tv = 1.0; if ( (j != 0) && (tFlag != 3) ) { // Compute relative measures of convergence of s and v sequences - tv = (vv != 0.0) ? std::abs((vv - ovv) / vv) : tv; - ts = (ss != 0.0) ? std::abs((ss - oss) / ss) : ts; + tv = (vv != 0.0) ? abs((vv - ovv) / vv) : tv; + ts = (ss != 0.0) ? abs((ss - oss) / ss) : ts; // If decreasing, multiply the two most recent convergence measures valueType tvv = (tv < otv) ? tv * otv : 1; valueType tss = (ts < ots) ? ts * ots : 1; @@ -569,7 +569,7 @@ namespace PolynomialRoots { /* // Scale if there are large or very small coefficients // Computes a scale factor to multiply the coefficients of the polynomial. - // The scaling is done to avoid overflow and to avoid undetected underflow + // The scaling is done to avoid overflow and to avoid undetected underflow // interfering with the convergence criterion. // The factor is a power of the base. */ @@ -585,7 +585,7 @@ namespace PolynomialRoots { .. The frexp() functions break the floating-point number value into .. a normalized fraction and an integral power of 2. .. They store the integer in the int object pointed to by exp. - .. The functions return a number x such that x has a magnitude in + .. The functions return a number x such that x has a magnitude in .. the interval [1/2, 1) or 0, and value = x*(2**exp). */ int max_exponent = std::numeric_limits::min(); @@ -614,7 +614,7 @@ namespace PolynomialRoots { valueType pt[N+1]; #endif - for ( indexType i = 0; i < N; ++i ) pt[i] = std::abs(p[i]); + for ( indexType i = 0; i < N; ++i ) pt[i] = abs(p[i]); pt[N] = -abs(p[N]); // Compute upper estimate of bound @@ -634,7 +634,7 @@ namespace PolynomialRoots { evalPoly( x, pt, N, f, df ); dx = f / df; x -= dx; - } while ( std::abs(dx) > std::abs(x)*0.005 ); + } while ( abs(dx) > abs(x)*0.005 ); return x; } @@ -661,7 +661,7 @@ namespace PolynomialRoots { } else if ( Degree == 3 ) { Cubic solve( p[0], p[1], p[2], p[3] ); switch ( solve.numRoots() ) { - case 3: solve.getRoot2( zeror[2], zeroi[2] ); + case 3: solve.getRoot1( zeror[2], zeroi[2] ); case 2: solve.getRoot1( zeror[1], zeroi[1] ); case 1: solve.getRoot0( zeror[0], zeroi[0] ); } @@ -697,7 +697,7 @@ namespace PolynomialRoots { #endif int N = Degree; - valueType xx = std::sqrt(0.5); //= 0.70710678 + valueType xx = sqrt(0.5); //= 0.70710678 valueType yy = -xx; // Remove zeros at the origin, if any @@ -735,7 +735,7 @@ namespace PolynomialRoots { K[j] = t * K[j-1] + p[j]; } K[0] = p[0]; - zerok = std::abs(K[NM1]) <= std::abs(bb) * epsilon10; + zerok = abs(K[NM1]) <= abs(bb) * epsilon10; } } diff --git a/toolbox/src/PolynomialRoots-Utils.cc b/toolbox/src/PolynomialRoots-Utils.cc index 2b02f8b9..4af6b3f5 100644 --- a/toolbox/src/PolynomialRoots-Utils.cc +++ b/toolbox/src/PolynomialRoots-Utils.cc @@ -4,7 +4,7 @@ | | | , __ , __ | | /|/ \ /|/ \ | - | | __/ _ ,_ | __/ _ ,_ | + | | __/ _ ,_ | __/ _ ,_ | | | \|/ / | | | | \|/ / | | | | | |(__/|__/ |_/ \_/|/|(__/|__/ |_/ \_/|/ | | /| /| | @@ -63,7 +63,7 @@ namespace PolynomialRoots { return res; } } - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // // op[0] * x^n + .... + op[n-1]*x + op[n] @@ -89,9 +89,9 @@ namespace PolynomialRoots { return res; } } - + //============================================================================ - + /* .. scale roots by pair.second, after scaling the polinomial has the form .. @@ -100,7 +100,7 @@ namespace PolynomialRoots { .. pair.first is the index such that p[pair.first] = +-1 .. */ - + static pair scalePolynomial( @@ -114,7 +114,7 @@ namespace PolynomialRoots { indexType i = n; while ( --i >= 0 ) { ps[i] = p[i]/an; - valueType scale1 = std::pow( std::abs(ps[i]), 1.0/(n-i) ); + valueType scale1 = std::pow( abs(ps[i]), 1.0/(n-i) ); if ( scale1 > scale ) { scale = scale1; i_max = i; } } // scale coeffs @@ -124,13 +124,13 @@ namespace PolynomialRoots { ps[n] = 1; return res; } - + //============================================================================ - + /* .. divide a(x) by (x-r) so that a(x) = (x-r) * b(x) */ - + static void deflatePolynomial( @@ -144,10 +144,10 @@ namespace PolynomialRoots { // Practical problems arising in the solution of polynomial equations. // J. Inst. Math. Appl. 8 (1971), 16–35. indexType i_cross = 0; - valueType v_cross = std::abs(a[0]); + valueType v_cross = abs(a[0]); valueType r1 = r; for ( indexType i = 1; i < n; ++i, r1 *= r ) { - valueType v_cross1 = std::abs(a[i]*r1); + valueType v_cross1 = abs(a[i]*r1); if ( v_cross1 > v_cross ) { v_cross = v_cross1; i_cross = i; } } @@ -206,7 +206,7 @@ namespace PolynomialRoots { return res; } #endif - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // // x^3 + A x^2 + B x + C @@ -223,11 +223,11 @@ namespace PolynomialRoots { indexType & i_case, valueType & scale ) { - - valueType a = std::abs(A); - valueType b = std::sqrt(abs(B)); - valueType c = std::cbrt(abs(C)); - + + valueType a = abs(A); + valueType b = sqrt(abs(B)); + valueType c = cbrt(abs(C)); + if ( a < b ) { if ( b < c ) i_case = 0; // a < b < c --> c MAX else i_case = 1; // a < b and c <= b --> b MAX @@ -235,7 +235,7 @@ namespace PolynomialRoots { if ( a < c ) i_case = 0; // b <= a < c --> c MAX else i_case = 2; // b <= a and c <= a --> a MAX } - + switch ( i_case ) { case 0: scale = c; @@ -257,7 +257,7 @@ namespace PolynomialRoots { break; } } - + // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // // a3*x^3 + a2*x^2 + a1*x + a0 = (x-r)*(a3*x^2+b1*x+b0) @@ -274,12 +274,12 @@ 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 ) i_cross = 3; switch ( i_cross ) { case 0: b1 = a2+a3*r; b0 = a1+r*b1; break;