Skip to content

Commit

Permalink
Further improved bound and skipped bisection for isolated roots withi…
Browse files Browse the repository at this point in the history
…n interval size

Combined the Lagrange and Cauchy bounds to return the smallest. realRootIsolation no longer bisects intervals when the root is isolated and the interval is below the required size (e.g. when two roots are very close together, the interval widths with be smaller, but the interval widths do not need to be smaller for the other roots). Final list is also sorted in case bisection stops in an interval after the first.

Test also corrected (the new isolating intervals only look 'less nice' because the previous initial bounds happened to be divisible by 8).
  • Loading branch information
cel34-bath committed Feb 15, 2024
1 parent c044b8a commit 15d28bb
Showing 1 changed file with 17 additions and 18 deletions.
35 changes: 17 additions & 18 deletions M2/Macaulay2/packages/RealRoots.m2
Original file line number Diff line number Diff line change
Expand Up @@ -398,26 +398,31 @@ realRootIsolation (RingElement,A) := List => (f,r)->(

--bound for real roots
C := (listForm ((f-leadTerm(f))/leadCoefficient(f)))/last; --make the polynomial monic, and obtain list of coefficients of non-lead monomials.
M := max(1,sum(C,abs)); --obtains Lagrange bound.
M := min(1+max(apply(C,abs)),max(1,sum(C,abs))); --obtains Cauchy or Lagrange bound.

L := {{-M,M}};
midp := 0;
v := new MutableHashTable from {M=>variations apply(l,g->signAt(g,M)),-M=>variations apply(l,g->signAt(g,-M))};

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 (
{}
)
Expand Down Expand Up @@ -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 ///
Expand Down Expand Up @@ -1073,10 +1078,4 @@ TEST ///
-- assert(rationalUnivariateRepresentationresentation(I) == {Z^2 - (3/2)*Z - 9, x + y});
-- ///

end






end

0 comments on commit 15d28bb

Please sign in to comment.