diff --git a/M2/Macaulay2/packages/RealRoots.m2 b/M2/Macaulay2/packages/RealRoots.m2 index a5d7e140be9..383ef531193 100644 --- a/M2/Macaulay2/packages/RealRoots.m2 +++ b/M2/Macaulay2/packages/RealRoots.m2 @@ -391,14 +391,14 @@ realRootIsolation (RingElement,A) := List => (f,r)->( if not isUnivariatePolynomial(f) then error "Error: Expected univariate polynomial"; R := ring f; - f = sub(f/gcd(f,diff(variable f,f)),R); + f = sub(f/gcd(f,diff(variable f,f)),R); --makes polynomial squarefree if (SturmCount(f)>0) then ( l := SturmSequence(f); --bound for real roots - C := (listForm f)/last; - M := (sum(C,abs))/(leadCoefficient f); + C := (listForm ((f-leadTerm(f))/leadCoefficient(f)))/last; --make the polynomial monic, and obtain list of coefficients of non-lead monomials. + M := min(1+max(apply(C,abs)),max(1,sum(C,abs))); --obtains Cauchy or Lagrange bound. L := {{-M,M}}; midp := 0; @@ -406,18 +406,23 @@ realRootIsolation (RingElement,A) := List => (f,r)->( while (max apply(L,I-> I#1-I#0) > r) or (max apply(L,I-> v#(I#0)-v#(I#1)) > 1) do ( for I in L do ( - midp = (sum I)/2; - v#midp = variations apply(l,g->signAt(g,midp)); - L = drop(L,1); - if (v#(I#0)-v#midp > 0) then ( - L = append(L,{I#0,midp}) - ); - if (v#midp-v#(I#1) > 0) then ( - L = append(L,{midp,I#1}) + if ((v#(I#0)-v#(I#1) == 1) and (I#1-I#0 <= r)) then ( + L = take(L,{1,#L})|{L#0}; -- skip bisection if root is identified and bound is within interval size. + ) + else ( + midp = (sum I)/2; + v#midp = variations apply(l,g->signAt(g,midp)); + L = drop(L,1); + if (v#(I#0)-v#midp > 0) then ( + L = append(L,{I#0,midp}) + ); + if (v#midp-v#(I#1) > 0) then ( + L = append(L,{midp,I#1}) + ); ) ) ); - L + sort L --list is ordered from smallest to largest interval, in case the while look stops after a larger interval ) else ( {} ) @@ -1033,7 +1038,7 @@ TEST /// TEST /// R = QQ[t]; f = (t-1)^2*(t+3)*(t+5)*(t-6); - assert(realRootIsolation(f,1/2) == {{-161/32, -299/64}, {-207/64, -23/8}, {23/32, 69/64}, {23/4, 391/64}}); + assert(realRootIsolation(f,1/2) == {{-1365/256,-637/128},{-819/256,-91/32},{91/128,273/256},{91/16,1547/256}}); /// TEST /// @@ -1073,10 +1078,4 @@ TEST /// -- assert(rationalUnivariateRepresentationresentation(I) == {Z^2 - (3/2)*Z - 9, x + y}); -- /// -end - - - - - - +end \ No newline at end of file