Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

correction to RealRoots:-realRootIsolation #3100

Merged
merged 2 commits into from
Feb 15, 2024
Merged

correction to RealRoots:-realRootIsolation #3100

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c044b8a
correction to RealRoots:-realRootIsolation
cel34-bath Feb 2, 2024
15d28bb
Further improved bound and skipped bisection for isolated roots withi…
cel34-bath Feb 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 19 additions & 20 deletions M2/Macaulay2/packages/RealRoots.m2
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -391,33 +391,38 @@ 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;
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
Loading