diff --git a/PolynomialRoots.cpp b/PolynomialRoots.cpp index 38a01f0..6b1939b 100644 --- a/PolynomialRoots.cpp +++ b/PolynomialRoots.cpp @@ -201,14 +201,14 @@ QuadraticPolynomialRoots( DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) vdRoot.push_back( 0) ; vdRoot.push_back( - vdPoly[1] / vdPoly[2]) ; if ( vdRoot[0] < vdRoot[1]) - swap( vdRoot[0], vdRoot[1]) ; + swap( vdRoot[0], vdRoot[1]) ; return 2 ; } - // caso generale + // caso generale ( x^2 + 2 p x + q = 0) double p = vdPoly[1] / ( 2.0 * vdPoly[2]) ; double q = vdPoly[0] / vdPoly[2] ; double Delta = p * p - q ; - if ( abs( Delta) < DBL_EPSILON) { + if ( abs( Delta) < DBL_EPSILON * DBL_EPSILON) { vdRoot.clear() ; vdRoot.reserve( 1) ; vdRoot.push_back( -p) ; @@ -219,7 +219,11 @@ QuadraticPolynomialRoots( DBLVECTOR& vdPoly, DBLVECTOR& vdRoot, int* pnIter) } // Delta positivo double dSqrtD = sqrt( Delta) ; - vdRoot.push_back( - p + dSqrtD) ; - vdRoot.push_back( - p - dSqrtD) ; + double r1 = -( p + copysign( dSqrtD, p)) ; + double r2 = q / r1 ; + if ( r1 > r2) + swap( r1, r2) ; + vdRoot.push_back( r1) ; + vdRoot.push_back( r2) ; return 2 ; }