polexpr (9dec18)

git-svn-id: svn://tug.org/texlive/trunk@49365 c570f23f-e606-0410-a88d-b1316a301751

1 files changed, 888 insertions, 280 deletions

diff --git a/Master/texmf-dist/tex/latex/polexpr/polexpr.sty b/Master/texmf-dist/tex/latex/polexpr/polexpr.sty
index f163fe10f29..601446cf671 100644
--- a/Master/texmf-dist/tex/latex/polexpr/polexpr.sty
+++ b/Master/texmf-dist/tex/latex/polexpr/polexpr.sty
@@ -1,7 +1,7 @@
 % author: Jean-François Burnol
 % License: LPPL 1.3c (author-maintained)
 \ProvidesPackage{polexpr}%
-  [2018/11/20 v0.6 Polynomial expressions with rational coefficients (JFB)]%
+  [2018/12/08 v0.7 Polynomial expressions with rational coefficients (JFB)]%
 \RequirePackage{xintexpr}[2018/06/17]% xint 1.3c for \ifxintglobaldefs boolean
 \edef\POL@restorecatcodes
     {\catcode`\noexpand\_ \the\catcode`\_ %
@@ -28,11 +28,28 @@
 }%
 
 %% AUXILIARIES
+\catcode`! 3
+%% added at 0.7
+\newcommand\polexprsetup[1]{\POL@setup_parsekeys #1,=!,\xint_bye}%
+\def\POL@setup_parsekeys #1=#2#3,{%
+    \ifx!#2\expandafter\xint_bye\fi
+    \csname POL@setup_setkey_\xint_zapspaces #1 \xint_gobble_i\endcsname
+    \xint_firstoftwo
+    {\PackageWarning{polexpr}{The \detokenize{#1} key is unknown! ignoring}}%
+    {\xintZapLastSpaces{#2#3}}%
+    \POL@setup_parsekeys
+}%
+\catcode`! 11 
+\def\POL@setup_setkey_norr     #1#2{\edef\POL@norr}%
+\def\POL@setup_setkey_sqfnorr  #1#2{\edef\POL@sqfnorr}%
+\polexprsetup{norr=_norr, sqfnorr=_sqf_norr}
+
 \newcount\POL@count
 \newif\ifPOL@pol
 \newif\ifxintveryverbose
 \newif\ifpoltypesetall
-\newif\ifPOL@sturm@declareunnormalized
+\newif\ifPOL@tosturm@makefirstprimitive
+\POL@tosturm@makefirstprimitivetrue
 \newif\ifPOL@isolz@nextwillneedrefine
 \newif\ifpoltoexprall
 %% the main exchange structure (stored in macros \POLuserpol@<name>)
@@ -681,9 +698,16 @@
 \def\POL@makeprim@macro#1%
    {\xintREZ{\xintNum{\xintDiv{#1}{\POL@makeprim@icontent}}}}%
 \newcommand\PolMakePrimitive[1]{%
+    % This does not need a full user declared polynomial on input, only
+    % a \POLuserpol@name macro, but on output it is fully declared
     \edef\POL@makeprim@icontent{\PolIContent{#1}}%
     \PolMapCoeffs\POL@makeprim@macro{#1}%
 }%
+\def\POL@makeprimitive#1{%
+    % Avoids declaring the polynomial, internal usage in \PolToSturm
+    \edef\POL@makeprim@icontent{\PolIContent{#1}}%
+    \POL@mapcoeffs\POL@makeprim@macro{#1}%
+}%
 
 
 %% Sturm Algorithm (polexpr 0.4)
@@ -696,15 +720,16 @@
 %%       holding the coefficients in memory
 %% 0.6 fixes the case of a constant polynomial P which caused division
 %%     by zero error from P'.
-\newcommand\PolToSturm{\@ifstar
-   {\POL@sturm@declareunnormalizedtrue\POL@ToSturm}%
-   {\POL@sturm@declareunnormalizedfalse\POL@ToSturm}%
-}%
+\newcommand\PolToSturm{\@ifstar{\PolToSturm@@}{\PolToSturm@}}%
 \def\POL@aux@toint#1{\xintREZ{\xintNum{#1}}}% for polynomials with int. coeffs!
-\def\POL@ToSturm#1#2{%
+%% Attention that some macros rely upon this one setting \POL@sturmname
+%% and \POL@sturm@N as it does
+\def\PolToSturm@#1#2{%
   \edef\POL@sturmname{#2}%
   % 0.6 uses 2 underscores (one before index, one after) to keep in memory
   % the unnormalized chain
+  % This supposes #1 to be a genuine polynomial, not only a name with
+  % a \POLuserpol@#1 macro
   \POL@let{\POL@sturmname _0_}{#1}%
   \ifnum\PolDegree{#1}=\z@
      \def\POL@sturm@N{0}%
@@ -713,7 +738,8 @@
      % if constant is negative. I also don't worry if polynomial is zero.
      \@namedef{POLuserpol@\POL@sturmname _0}{0.\empty{1/1[0]}}%
   \else
-     \POL@ToSturm@DoSturm
+     \ifPOL@tosturm@makefirstprimitive\POL@makeprimitive{\POL@sturmname _0_}\fi
+     \POL@tosturm@dosturm
   \fi
   \expandafter
   \let\csname PolSturmChainLength_\POL@sturmname\endcsname\POL@sturm@N
@@ -724,22 +750,21 @@
   \unless\ifnum\POL@sturm@N=\POL@count
     \advance\POL@count\@ne
   \repeat
+}%
+\def\PolToSturm@@#1#2{\PolToSturm@{#1}{#2}\POL@tosturm@declareunnormalized}%
+\def\POL@tosturm@declareunnormalized{%
   % optionally declare also the unnormalized ones
   \POL@count\z@
-  \ifPOL@sturm@declareunnormalized
-    \POL@count\z@
-    \xintloop
-        \POL@newpol{\POL@sturmname _\the\POL@count _}%
-    \unless\ifnum\POL@sturm@N=\POL@count
-        \advance\POL@count\@ne
-    \repeat
-  \fi
+  \xintloop
+      \POL@newpol{\POL@sturmname _\the\POL@count _}%
+  \unless\ifnum\POL@sturm@N=\POL@count
+      \advance\POL@count\@ne
+  \repeat
 }%
-\def\POL@ToSturm@DoSturm{%
-  \PolMakePrimitive{\POL@sturmname _0_}%
+\def\POL@tosturm@dosturm{%
   \POL@Diff@@one{\POL@sturmname _0_}{\POL@sturmname _1_}%
   % re-utiliser \POL@varcoeffs directement?
-  \PolMakePrimitive{\POL@sturmname _1_}%
+  \POL@makeprimitive{\POL@sturmname _1_}% does not do \POL@newpol
   \POL@count\@ne
   \xintloop
     \POL@divide{\POL@sturmname _\the\numexpr\POL@count-\@ne\relax _}%
@@ -750,6 +775,7 @@
     \expandafter\let
     \csname POLuserpol@\POL@sturmname _\the\POL@count _\endcsname\POL@R
     \edef\POL@makeprim@icontent{-\POL@icontent\POL@polR}%
+    % this avoids the \POL@newpol from \PolMapCoeffs
     \POL@mapcoeffs\POL@makeprim@macro{\POL@sturmname _\the\POL@count _}%
   \repeat
   \edef\POL@sturm@N{\the\POL@count}%
@@ -763,13 +789,13 @@
       \expandafter
        \let\csname POLuserpol@\POL@sturmname _\the\POL@count\endcsname\POL@Q
       % quotient actually belongs to Z[X] and is primitive
-      \POL@mapcoeffs{\POL@aux@toint}{\POL@sturmname _\the\POL@count}%    
+      \POL@mapcoeffs\POL@aux@toint{\POL@sturmname _\the\POL@count}%    
     \ifnum\POL@count>\z@
     \repeat
     \@namedef{POLuserpol@\POL@sturmname _\POL@sturm@N}{0.\empty{1/1[0]}}%
   \else % they are already normalized
-      \advance\POL@count\@ne % attention to include last one also
-      \xintloop
+    \advance\POL@count\@ne % attention to include last one also
+    \xintloop
       \advance\POL@count\m@ne
       \expandafter\let
        \csname POLuserpol@\POL@sturmname _\the\POL@count\expandafter\endcsname
@@ -777,7 +803,7 @@
     \ifnum\POL@count>\z@
     \repeat
   \fi
-  % Back to \POL@ToSturm
+  % Back to \PolToSturm@, \POL@count holds 0
 }%
 \newcommand\PolSturmChainLength[1]
     {\romannumeral`^^@\csname PolSturmChainLength_#1\endcsname}%
@@ -791,19 +817,19 @@
 }%
 \def\POL@sturmchain@getSV@at#1{% ATTENTION USES \POL@count
   \def\POL@sturmchain@SV{0}%
-  \edef\POL@sturmchain@sign{\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{#1}}}%
+  \edef\POL@sturmchain@sign{\xintiiSgn{\POL@eval{\POL@sturmname _0}{#1}}}%
   \let\POL@isolz@lastsign\POL@sturmchain@sign
   \POL@count \z@
   \ifnum\POL@isolz@lastsign=\z@
     \edef\POL@isolz@lastsign
-        {\xintiiSgn{\Pol@Eval{\POL@sturmname _1}{#1}}}%
+        {\xintiiSgn{\POL@eval{\POL@sturmname _1}{#1}}}%
     \POL@count \@ne
   \fi
   \xintloop
   \unless\ifnum\POL@sturmlength=\POL@count
   \advance\POL@count \@ne
     \edef\POL@isolz@newsign
-      {\xintiiSgn{\Pol@Eval{\POL@sturmname _\the\POL@count}{#1}}}%
+      {\xintiiSgn{\POL@eval{\POL@sturmname _\the\POL@count}{#1}}}%
     \ifnum\POL@isolz@newsign=\numexpr-\POL@isolz@lastsign\relax
       \edef\POL@sturmchain@SV{\the\numexpr\POL@sturmchain@SV+\@ne}%
       \let\POL@isolz@lastsign=\POL@isolz@newsign
@@ -827,11 +853,431 @@
 }%
 
 
+% 0.6 added starred variant to count multiplicities
+% 0.7 added double starred variant to locate all rational roots
 \newcommand\PolSturmIsolateZeros{\@ifstar
       {\PolSturmIsolateZerosAndGetMultiplicities}%
       {\PolSturmIsolateZeros@}%
 }%
-\newcommand\PolSturmIsolateZerosAndGetMultiplicities[2][\empty]{%
+\newcommand\PolSturmIsolateZerosAndGetMultiplicities{\@ifstar
+      {\PolSturmIsolateZerosGetMultiplicitiesAndRationalRoots}%
+      {\PolSturmIsolateZerosAndGetMultiplicities@}%
+}%
+% on aurait besoin de ça dans xint, mais il aurait un \xintRaw{#1} alors
+\def\POL@xintfrac@getNDE #1%
+   {\expandafter\POL@xintfrac@getNDE@i\romannumeral`^^@#1}%
+\def\POL@xintfrac@getNDE@i #1/#2[#3]#4#5#6{\def#4{#1}\def#5{#2}\def#6{#3}}%
+\newcommand\PolSturmIsolateZerosGetMultiplicitiesAndRationalRoots[2][\empty]{%
+  \PolSturmIsolateZerosAndFindRationalRoots[#1]{#2}%
+  \ifnum\POL@isolz@NbOfRoots>\z@
+     % get multiplicities of irrational (real) roots, if any
+     \ifnum\POL@findrat@nbofirrroots>\z@
+        \POL@findrat@getirrmult
+     \fi
+     \POL@isolzmult@defvar@M
+  \fi
+}%
+% added at 0.7
+\newcommand\PolSturmIsolateZerosAndFindRationalRoots[2][\empty]{%
+  % #1 optional E such that roots are searched in -10^E < x < 10^E
+  % both -10^E and +10^E must not be roots!
+  % #2 name of Sturm chain (already pre-computed)
+  \edef\POL@sturmname{#2}%
+  \edef\POL@sturm@N{\@nameuse{PolSturmChainLength_\POL@sturmname}}%
+  % isolate the roots (detects case of constant polynomial)
+  \PolSturmIsolateZeros@{\POL@sturmname}%
+  \ifnum\POL@isolz@NbOfRoots=\z@
+     % no real roots, define empty arrays nevertheless
+     \begingroup\globaldefs\@ne
+     \expandafter\xintAssignArray\expandafter\to\csname POL_ZeroMult\POL@sturmname\endcsname
+     \expandafter\xintAssignArray\expandafter\to\csname POL_RRIndex\POL@sturmname\endcsname
+     \endgroup
+  \else
+     % all we currently know is that multiplicities are at least one
+     \begingroup\globaldefs\@ne
+     \expandafter\POL@initarray\csname POL_ZeroMult\POL@sturmname\endcsname{1}%
+     \endgroup
+     % on ne va pas utiliser de Horner, mais des divisions par X - x, et ces
+     % choses vont évoluer, ainsi que le coefficient dominant entier
+     % (pour \POL@divide entre autres if faut des noms de user pol)
+     \expandafter\let
+         \csname POLuserpol@\POL@sturmname\POL@sqfnorr\expandafter\endcsname
+         \csname POLuserpol@\POL@sturmname _0\endcsname
+     \expandafter\let
+         \csname POLuserpol@\POL@sturmname\POL@norr\expandafter\endcsname
+         \csname POLuserpol@\POL@sturmname _0_\endcsname
+     % attention formé avec\xintREZ d'où le \xintAbs pas \xintiiAbs
+     % D and its exponent E will get updated along the way
+     \edef\POL@findrat@D{\xintAbs{\PolLeadingCoeff{\POL@sturmname _0}}}%
+     \POL@xintfrac@getNDE\POL@findrat@D\POL@findrat@Dint\POL@_\POL@findrat@Dexp
+     \xintiiifOne{\POL@findrat@Dint}
+        {\let\POL@findrat@E\POL@findrat@Dexp} % aussi ok pour 1[0]
+        {\edef\POL@findrat@E{\the\numexpr\xintLen{\POL@findrat@Dint}%
+                                                 +\POL@findrat@Dexp}}%
+     \POL@initarray\POL@IfMultIsKnown\xint_secondoftwo
+     \let\POL@findrat@nbofirrroots\POL@isolz@NbOfRoots
+     % find all rational roots, and their multiplicities,
+     % factor them out in passing from original (Sturm root) polynomial
+     \ifnum\POL@findrat@E<7
+         \PolEnsureIntervalLength{\POL@sturmname}{1}{-\POL@findrat@E}%
+         \def\POL@findrat@index{1}%
+         \POL@findrat@loop@secondpass@direct
+     \else
+     % we do a first pass scanning for "small" roots p/q (i.e. q < 1000)
+         \def\POL@findrat@index{1}%
+         \POL@findrat@loop@firstpass
+     % and now we do the final pass finding them all
+         \def\POL@findrat@index{1}%
+         \PolEnsureIntervalLength{\POL@sturmname}{1}{-\POL@findrat@E}%
+         \POL@findrat@loop@secondpass
+     \fi
+     % declare the new polynomials
+     \POL@newpol{\POL@sturmname\POL@sqfnorr}%  without multiplicities
+     \POL@newpol{\POL@sturmname\POL@norr}% with multiplicities
+     % declare the array holding the interval indices for the rational roots
+     \expandafter\POL@findrat@doRRarray\csname POL_RRIndex\POL@sturmname\endcsname
+  \fi
+}%
+\def\POL@findrat@doRRarray#1{%
+  % il faudrait un \xintAssignArray* qui fasse même expansion que \xintFor*
+  \edef\POL@temp{%
+    \xintiloop[1+1]
+    \romannumeral0\csname POL_ZeroIsKnown\POL@sturmname\xintiloopindex\endcsname
+    \xintbracediloopindex % I should have named it \xintiloopbracedindex...
+    {}%
+    \ifnum\xintiloopindex<\POL@isolz@NbOfRoots\space
+    \repeat }%
+  \begingroup\globaldefs1
+  % attention de ne surtout pas faire un \expandafter ici, car en cas d'un
+  % seul item, \xintAssignArray l'unbraces...
+  \xintAssignArray\POL@temp\to#1%
+  \endgroup
+}%
+\def\POL@findrat@loop@firstpass{%
+  \PolSturmIfZeroExactlyKnown{\POL@sturmname}{\POL@findrat@index}%
+     \POL@findrat@loop@decimal% get its multiplicity
+     \POL@findrat@loop@a
+  \edef\POL@findrat@index{\the\numexpr\POL@findrat@index+\@ne}%
+  \ifnum\POL@findrat@index>\POL@isolz@NbOfRoots
+  \else
+    \expandafter\POL@findrat@loop@firstpass
+  \fi
+}%
+\def\POL@findrat@loop@secondpass{%
+  \PolSturmIfZeroExactlyKnown{\POL@sturmname}{\POL@findrat@index}%
+     {}% nothing more to be done, already stored
+     \POL@findrat@loop@b
+  \edef\POL@findrat@index{\the\numexpr\POL@findrat@index+\@ne}%
+  \ifnum\POL@findrat@index>\POL@isolz@NbOfRoots
+  \else
+     \PolEnsureIntervalLength
+         {\POL@sturmname}{\POL@findrat@index}{-\POL@findrat@E}% dynamic
+     \expandafter\POL@findrat@loop@secondpass
+  \fi
+}%
+\def\POL@findrat@loop@secondpass@direct{%
+  \PolSturmIfZeroExactlyKnown{\POL@sturmname}{\POL@findrat@index}%
+     \POL@findrat@loop@decimal
+     \POL@findrat@loop@b
+  \edef\POL@findrat@index{\the\numexpr\POL@findrat@index+\@ne}%
+  \ifnum\POL@findrat@index>\POL@isolz@NbOfRoots
+  \else
+     \PolEnsureIntervalLength
+         {\POL@sturmname}{\POL@findrat@index}{-\POL@findrat@E}% dynamic
+     \expandafter\POL@findrat@loop@secondpass@direct
+  \fi
+}%
+\def\POL@findrat@loop@decimal{% we have an already found decimal root
+  % we do not go via @storeit, as it is already stored
+  % j'ai beaucoup hésité néanmoins, car je pourrais faire \xintIrr ici,
+  % mais attention aussi à l'interaction avec le \PolDecToString. Les racines
+  % trouvées directement (qui peuvent être des nombres décimaux) sont elles
+  % stockées comme fraction irréductibles (modulo action additionnelle de 
+  % \PolDecToString). 
+  \POL@xintfrac@getNDE
+    {\xintIrr{\POL@xintexprGetVar{\POL@sturmname L_\POL@findrat@index}}[0]}%
+    \POL@findrat@xN\POL@findrat@xD\POl@_
+  % we can't move this to updatequotients because other branch will
+  % need to do the division first anyhow
+  \edef\POLuserpol@_findrat@oneterm{1.\noexpand\empty
+        {\xintiiOpp\POL@findrat@xN/1[0]}{\POL@findrat@xD/1[0]}}%
+  \POL@divide{\POL@sturmname\POL@sqfnorr}{_findrat@oneterm}% the one without mult.
+  \POL@findrat@loop@updatequotients
+  \POL@findrat@loop@getmultiplicity
+}%
+% lacking from xint 1.3c, but \xintSgn has overhead, so we define ii version
+\def\xintiiifNeg{\romannumeral0\xintiiifneg }%
+\def\xintiiifneg #1%
+{%
+    \ifcase \xintiiSgn{#1}
+               \expandafter\xint_stop_atsecondoftwo
+            \or\expandafter\xint_stop_atsecondoftwo
+          \else\expandafter\xint_stop_atfirstoftwo
+    \fi
+}%
+\def\POL@findrat@getE #1/1[#2]{#2}% /1 as it should be there.
+% so an error will arise if not but cf \POL@refine@getE where I did not put it
+\def\POL@findrat@loop@a{%
+  % we do a first pass to identify roots with denominators < 1000
+  \PolEnsureIntervalLength{\POL@sturmname}{\POL@findrat@index}{-6}%
+  % attention that the width may have been already smaller than 10^{-6}
+  % also attention that one of the bound may be zero
+  \POL@get@IsoLeft@rawin
+  \POL@get@IsoRight@rawin
+  \edef\POL@findrat@localW
+      {\the\numexpr-\expandafter\POL@findrat@getE
+  % do I really need the \xintREZ?
+       \romannumeral0\xintrez
+         {\xintSub{\POL@IsoRight@rawin}{\POL@IsoLeft@rawin}}%
+      }% at least 6, maybe larger
+  \expandafter\POL@get@Int@aux
+              \POL@IsoLeft@rawin\POL@IsoLeft@Int{-\POL@findrat@localW}%
+  \expandafter\POL@get@Int@aux
+              \POL@IsoRight@rawin\POL@IsoRight@Int{-\POL@findrat@localW}%
+  % in case of odd, some waste here
+  \edef\POL@findrat@halflocalW{\the\numexpr(\POL@findrat@localW+1)/2-1}%
+  % Legendre Theorem will be used now but we separate a branch where
+  % everything can be done with \numexpr
+  \ifnum\POL@findrat@localW>10
+    % not implemented yet by lazyness!
+    % this root will be handled in second pass only
+  \else
+    \POL@findrat@gcdloop
+  \fi
+}%
+\def\POL@findrat@gcdloop{%
+  % we must be careful with sign (attention one of the bounds may be zero)
+  \let\POL@findrat@IfNeg\xint_secondoftwo
+  \xintiiifSgn\POL@IsoLeft@Int
+    \POL@findrat@gcdloop@n
+    \POL@findrat@gcdloop@zero
+    \POL@findrat@gcdloop@p
+}%
+\def\POL@findrat@gcdloop@n{%
+  \let\POL@findrat@IfNeg\xint_firstoftwo
+  \let\POL@temp\POL@IsoRight@Int
+  \edef\POL@IsoRight@Int{\xintiiOpp{\POL@IsoLeft@Int}}%
+  \edef\POL@IsoLeft@Int{\xintiiOpp{\POL@temp}}%
+  \xintiiifSgn\POL@IsoLeft@Int
+     \POL@error % impossible branch
+     \POL@findrat@gcdloop@zero
+     \POL@findrat@gcdloop@p
+}%
+\def\POL@findrat@gcdloop@zero{%
+  % the continued fraction would be the one of 1/2, so only 1/2 to test...
+  \edef\POL@findrat@x
+      {1/2\romannumeral\xintreplicate{\POL@findrat@localW}{0}[0]}%
+  \POL@findrat@gcdloop@testit
+}%
+\def\POL@findrat@gcdloop@p{%
+  \edef\POL@findrat@gcdloop@Ap{\xintDec{\xintDouble\POL@IsoRight@Int}}%
+  \edef\POL@findrat@gcdloop@A
+      {2\romannumeral\xintreplicate\POL@findrat@localW{0}}%
+  \xintAssign
+    \xintiiDivision\POL@findrat@gcdloop@Ap\POL@findrat@gcdloop@A
+  \to\POL@findrat@gcdloop@B\POL@findrat@gcdloop@An
+  % on fait de la tambouille pour n'utiliser que \numexpr par la suite
+  % le reste @An est < 2.10^10 au pire donc ok pour \numexpr
+  % we will drop integral part in our updating P
+  \let\POL@findrat@gcdloop@Binitial\POL@findrat@gcdloop@B
+  \def\POL@findrat@gcdloop@B{0}%  do as if B1 = 0
+  \def\POL@findrat@gcdloop@Pp{1}% P0
+  \def\POL@findrat@gcdloop@P{0}%  P1
+  \def\POL@findrat@gcdloop@Qp{0}% Q0
+  \def\POL@findrat@gcdloop@Q{1}%  Q1
+  % A2=An can not be zero, as Ap (=A0) is odd and A (=A1=200...000) is even
+  % first Binitial + P1/Q1 ( = Binitial) can not be root
+  \let\POL@findrat@gcdloop@Ap\POL@findrat@gcdloop@A % A1
+  \let\POL@findrat@gcdloop@A\POL@findrat@gcdloop@An % A2
+  \def\next{\POL@findrat@gcdloop@update}%
+  \def\POL@findrat@gcdloop@done{0}%
+  \POL@findrat@gcdloop@body
+}%
+\def\POL@findrat@gcdloop@body{%
+    % annoying that \numexpr has no divmod... use counts? but groups annoying
+    \edef\POL@findrat@gcdloop@B
+        {\the\numexpr(\POL@findrat@gcdloop@Ap+\POL@findrat@gcdloop@A/2)/%
+                      \POL@findrat@gcdloop@A - \@ne}%
+    \edef\POL@findrat@gcdloop@An
+        {\the\numexpr\POL@findrat@gcdloop@Ap-%
+                     \POL@findrat@gcdloop@B*\POL@findrat@gcdloop@A}%
+    \edef\POL@findrat@gcdloop@Pn
+        {\the\numexpr\POL@findrat@gcdloop@Pp+%
+                     \POL@findrat@gcdloop@B*\POL@findrat@gcdloop@P}%
+    \edef\POL@findrat@gcdloop@Qn
+        {\the\numexpr\POL@findrat@gcdloop@Qp+%
+                     \POL@findrat@gcdloop@B*\POL@findrat@gcdloop@Q}%
+    \ifnum\expandafter\xintLength\expandafter{\POL@findrat@gcdloop@Qn}%
+          >\POL@findrat@halflocalW\space
+       \let\next\empty % no solution was found
+    \else
+    % with these conditions on denom, only candidates are by Legendre
+    % theorem among the convergents as computed here
+       \ifnum\POL@findrat@gcdloop@Qn>\POL@findrat@gcdloop@An\space
+         % means that P/Q is in interval and is thus a candidate
+         % it is automatically irreducible
+         \edef\POL@findrat@x{\xintiiAdd
+           {\xintiiMul{\POL@findrat@gcdloop@Qn}{\POL@findrat@gcdloop@Binitial}}%
+           {\POL@findrat@gcdloop@Pn}/\POL@findrat@gcdloop@Qn[0]}%
+         \POL@findrat@gcdloop@testit
+         \if1\POL@findrat@gcdloop@done
+            \let\next\empty % a solution was found
+         \fi
+       \fi
+    \fi
+    \next
+}%
+\def\POL@findrat@gcdloop@update{%
+    \ifnum\POL@findrat@gcdloop@An>\z@
+        \let\POL@findrat@gcdloop@Ap\POL@findrat@gcdloop@A
+        \let\POL@findrat@gcdloop@A\POL@findrat@gcdloop@An
+        \let\POL@findrat@gcdloop@Pp\POL@findrat@gcdloop@P
+        \let\POL@findrat@gcdloop@P\POL@findrat@gcdloop@Pn
+        \let\POL@findrat@gcdloop@Qp\POL@findrat@gcdloop@Q
+        \let\POL@findrat@gcdloop@Q\POL@findrat@gcdloop@Qn
+    \expandafter\POL@findrat@gcdloop@body
+    \fi
+}%
+\def\POL@findrat@gcdloop@testit{%
+  % zero should never occur here
+  \POL@findrat@IfNeg{\edef\POL@findrat@x{-\POL@findrat@x}}{}%
+  \POL@xintfrac@getNDE\POL@findrat@x\POL@findrat@xN\POL@findrat@xD\POL@_
+  \edef\POLuserpol@_findrat@oneterm{1.\noexpand\empty
+      {\xintiiOpp{\POL@findrat@xN}/1[0]}{\POL@findrat@xD/1[0]}}%
+  \POL@divide{\POL@sturmname\POL@sqfnorr}{_findrat@oneterm}% the one without mult.
+  \expandafter\POL@split\POL@R;\POL@degR\POL@polR
+  \ifnum\POL@degR=\m@ne % found a root
+    \POL@findrat@loop@storeit
+    \POL@findrat@loop@updatequotients
+    \POL@findrat@loop@getmultiplicity % will continue updating the mult. one
+    \def\POL@findrat@gcdloop@done{1}%
+  \else
+  \fi
+}%            
+% This is second phase
+\def\POL@findrat@loop@b{%
+  \edef\POL@findrat@Lscaled{\xintMul{\POL@findrat@D}%
+       {\POL@xintexprGetVar{\POL@sturmname L_\POL@findrat@index}}}%
+  \edef\POL@findrat@Rscaled{\xintMul{\POL@findrat@D}%
+       {\POL@xintexprGetVar{\POL@sturmname R_\POL@findrat@index}}}%
+  \xintiiifNeg{\POL@findrat@Lscaled}% using ii version is an abuse
+    {% negative interval (right bound possibly zero!)
+     % truncate towards zero (i.e. to the right) the left bound
+     \edef\POL@findrat@Num{\xintNum{\POL@findrat@Lscaled}/1[0]}%
+     % interval boundaries are not root hence in case that was exact
+     % this will not be found as a root; check if in interval
+     \xintifLt\POL@findrat@Num\POL@findrat@Rscaled
+         \POL@findrat@loop@c
+         {}% iterate
+    }%
+    {% positive interval (left bound possibly zero!)
+     % truncate towards zero (i.e. to the left) the right bound
+     \edef\POL@findrat@Num{\xintNum{\POL@findrat@Rscaled}/1[0]}%
+     % check if in interval
+     \xintifGt\POL@findrat@Num\POL@findrat@Lscaled
+         \POL@findrat@loop@c
+         {}% iterate
+    }%
+}%
+\def\POL@findrat@loop@c{%
+  % safer to do the edef as \POL@findrat@x used later in storeit
+  \edef\POL@findrat@x{\xintIrr{\xintDiv\POL@findrat@Num\POL@findrat@D}[0]}%
+  \POL@xintfrac@getNDE\POL@findrat@x\POL@findrat@xN\POL@findrat@xD\POL@_
+  \edef\POLuserpol@_findrat@oneterm{1.\noexpand\empty
+      {\xintiiOpp{\POL@findrat@xN}/1[0]}{\POL@findrat@xD/1[0]}}%
+  \POL@divide{\POL@sturmname\POL@sqfnorr}{_findrat@oneterm}% the one without mult.
+  \expandafter\POL@split\POL@R;\POL@degR\POL@polR
+  \ifnum\POL@degR=\m@ne % found a root
+    \POL@findrat@loop@storeit
+    \POL@findrat@loop@updatequotients
+    \POL@findrat@loop@getmultiplicity % will continue updating the mult. one
+  \fi
+  % iterate
+}%
+\def\POL@findrat@loop@storeit{%
+  % update storage, I can not use storeleftandright here (due to rawout etc...)
+  \expandafter
+  \xdef\csname POL_ZeroInt\POL@sturmname L\POL@findrat@index\endcsname
+                                     {\PolDecToString{\POL@findrat@x}}%
+  \global\expandafter
+  \let\csname POL_ZeroInt\POL@sturmname R\POL@findrat@index\expandafter\endcsname
+       \csname POL_ZeroInt\POL@sturmname L\POL@findrat@index\endcsname
+  \global\expandafter
+  \let\csname POL_ZeroIsKnown\POL@sturmname\POL@findrat@index\endcsname
+        \xint_stop_atfirstoftwo
+  \begingroup\xintglobaldefstrue
+  \xintdefvar
+    \POL@sturmname L_\POL@findrat@index,%
+    \POL@sturmname R_\POL@findrat@index,%
+    \POL@sturmname Z_\POL@findrat@index _isknown
+    := qfrac(\POL@findrat@x),qfrac(\POL@findrat@x),1;%
+  \endgroup
+}%
+\def\POL@findrat@loop@updatequotients{%
+  % attention last division must have been one testing vanishing of\POL@sqfnorr
+  \expandafter\let\csname POLuserpol@\POL@sturmname\POL@sqfnorr\endcsname\POL@Q
+  % quotient belongs to Z[X] and is primitive
+  \POL@mapcoeffs\POL@aux@toint{\POL@sturmname\POL@sqfnorr}%
+  % update the one with multiplicities
+  \POL@divide{\POL@sturmname\POL@norr}{_findrat@oneterm}%
+  \expandafter\let\csname POLuserpol@\POL@sturmname\POL@norr\endcsname\POL@Q
+  \POL@mapcoeffs\POL@aux@toint{\POL@sturmname\POL@norr}
+  % updating of \POL@findrat@D at end of execution of getmultiplicity
+}%
+\def\POL@findrat@loop@getmultiplicity{%
+  % the one without multiplicity must not be divided again!
+  % check if we have remaining multiplicity
+  \POL@divide{\POL@sturmname\POL@norr}{_findrat@oneterm}%
+  \expandafter\POL@split\POL@R;\POL@degR\POL@polR
+  \ifnum\POL@degR=\m@ne % yes
+    \expandafter\let\csname POLuserpol@\POL@sturmname\POL@norr\endcsname\POL@Q
+    \POL@mapcoeffs\POL@aux@toint{\POL@sturmname\POL@norr}%
+    \expandafter
+    \xdef
+    \csname POL_ZeroMult\POL@sturmname\POL@findrat@index\endcsname
+      {\the\numexpr
+       \csname POL_ZeroMult\POL@sturmname\POL@findrat@index\endcsname+\@ne}%
+    \expandafter\POL@findrat@loop@getmultiplicity
+  \else 
+  % done with multiplicity for this rational root, update stuff
+    \edef\POL@findrat@nbofirrroots
+         {\the\numexpr\POL@findrat@nbofirrroots-\@ne}%
+    \@namedef{POL@IfMultIsKnown\POL@findrat@index}{\xint_firstoftwo}%
+    \edef\POL@findrat@D{\xintAbs{\PolLeadingCoeff{\POL@sturmname\POL@sqfnorr}}}%
+    \POL@xintfrac@getNDE\POL@findrat@D\POL@findrat@Dint\POL@_\POL@findrat@Dexp
+    \xintiiifOne{\POL@findrat@Dint}
+        {\let\POL@findrat@E\POL@findrat@Dexp} % aussi ok pour 1[0]
+        {\edef\POL@findrat@E{\the\numexpr\xintLen{\POL@findrat@Dint}%
+                                                 +\POL@findrat@Dexp}}%
+  \fi
+}%
+\def\POL@findrat@getirrmult{%
+  % first get the GCD of remaining pol with its derivative
+  \POL@divide{\POL@sturmname\POL@norr}{\POL@sturmname\POL@sqfnorr}%
+  \expandafter\let
+    % attention au _ (cf. grosse astuce pour \POL@isolzmult@loop)
+    \csname POLuserpol@@_1\POL@sturmname _\endcsname\POL@Q
+  \ifnum\PolDegree{@_1\POL@sturmname _}>\z@
+    % il reste des multiplicités (mais peut-être pour des racines complexes)
+    % (ou pour des racines en-dehors de l'intervalle optionnel)
+    % attention recyclage ici de \POL@isolzmult@loop qui dépend de
+    % la  grosse astuce avec \@gobble
+    \POL@makeprimitive{@_1\POL@sturmname _}%
+    \let\POL@originalsturmname\POL@sturmname
+    % trick to get isolzmult@loop to define @@lastGCD to @_1sturmname_
+    % because it will do \POL@sturmname _\POL@sturm@N _
+    \edef\POL@sturmname{@_1\POL@sturmname}%
+    \let\POL@sturm@N\@gobble% !
+    \let\POL@isolz@NbOfRoots@with_unknown_mult\POL@findrat@nbofirrroots
+    \POL@tosturm@makefirstprimitivefalse
+    \POL@isolzmult@loop
+    \POL@tosturm@makefirstprimitivetrue
+    \let\POL@sturmname\POL@originalsturmname
+  \fi
+}%
+
+
+\newcommand\PolSturmIsolateZerosAndGetMultiplicities@[2][\empty]{%
   % #1 optional E such that roots are searched in -10^E < x < 10^E
   % both -10^E and +10^E must not be roots!
   % #2 name of Sturm chain (already pre-computed)
@@ -845,9 +1291,6 @@
      \expandafter\xintAssignArray\expandafter\to\csname POL_ZeroMult\POL@sturmname\endcsname
      \endgroup
   \else
-     % store Sturm chain name for usage in the main loop
-     \let\POL@originalsturmname\POL@sturmname
-     \edef\POL@isolzmult@indices{\xintSeq{1}{\POL@isolz@NbOfRoots}}%
      % all we currently know is that multiplicities are at least one
      \begingroup\globaldefs\@ne
      \expandafter\POL@initarray\csname POL_ZeroMult\POL@sturmname\endcsname{1}%
@@ -857,61 +1300,92 @@
      \ifnum\PolDegree{\POL@sturmname _\POL@sturm@N _}>\z@
          % scratch array of flags to signal known multiplicities
          \POL@initarray\POL@IfMultIsKnown\xint_secondoftwo
+         % this count has utility for the case there are other roots
+         % either complex or outside interval (in case of optional argument)
          \let\POL@isolz@NbOfRoots@with_unknown_mult\POL@isolz@NbOfRoots
-         \expandafter\expandafter\expandafter\POL@isolzmult@loop
+         % store Sturm chain name, it is needed and altered in isolzmult@loop
+         \let\POL@originalsturmname\POL@sturmname
+         \POL@tosturm@makefirstprimitivefalse
+         \POL@isolzmult@loop
+         \POL@tosturm@makefirstprimitivetrue
+         \let\POL@sturmname\POL@originalsturmname
      \fi
+     \POL@isolzmult@defvar@M
   \fi
 }%
+\def\POL@isolzmult@defvar@M{%
+  % Attention that is used not only in ...GetMultiplicities@ but also
+  % in FindRationalRoots
+  \begingroup\xintglobaldefstrue
+  % added at 0.7
+  \let\x\POL@isolz@NbOfRoots
+    \xintloop
+      \xintdefvar \POL@sturmname M_\x
+                  := \csname POL_ZeroMult\POL@sturmname\x\endcsname ;%
+    \edef\x{\the\numexpr\x-\@ne}%
+    \ifnum\x>\z@
+    \repeat
+  \endgroup
+}%
 \def\POL@isolzmult@loop{%
-   % we are here only if last iteration gave a new PGCD still of degree > 0
-   % As 0.6 \PolToSturm keeps memory of unnormalized Sturm chain, we use the
-   % PGCD from last iteration and generate a new Sturm chain.
-   % ATTENTION: first argument of \PolToSturm MUST NOT CONTAIN \POL@sturmname
-   \let\POL@@sturmname\POL@sturmname
-   % ATTENTION: we could use an underscore prefix to the name, but attention
-   % to tacit multiplication if used in an expression; however \PolEvalAt
-   % does not use expression parsing as \PolEvalAtExpr so this would be
-   % relatively safe. We must also not overwrite privately used names
-   % by polexpr or xint... Using prefix @_1 appears safe. They will accumulate.
-   % As the loop may break at any moment, depending on original P, not only
-   % on current polynomial which is examined to see if it has zeros, it does
-   % not seem to make sense to think about interface to keep memory of all
-   % the defined polynomials.
-   % \POL@sturm@N supposedly the one from last iteration
-   \PolToSturm{\POL@@sturmname _\POL@sturm@N _}{@_1\POL@@sturmname}%
+   % we are here only if last iteration gave a new GCD still of degree > 0
+   % \POL@sturm@N is the one from last iteration
+   % Attention to not use \POL@sturmname directly in first arg. of \PolToSturm
+   % Attention that we need for the case of known roots also to have the last
+   % GCD (with its multiplicities) known as a genuine polynomial
+   % - because of usage of \POL@eval in @isknown branch
+   % - because \PolToSturm@ does a \POL@let which would be anomalous
+   %   if the extended structure is not existing
+   \edef\POL@isolzmult@lastGCD{\POL@sturmname _\POL@sturm@N _}%
+   \edef\POL@isolzmult@newsturmname{@_1\POL@sturmname}%
+   \POL@newpol{\POL@isolzmult@lastGCD}%
+   \PolToSturm@{\POL@isolzmult@lastGCD}{\POL@isolzmult@newsturmname}%
    % now both \POL@sturmname and \POL@sturm@N have changed
-   % if GCD is now a constant, we will not come back here
-   \edef\POL@sturmfinaldeg{\PolDegree{\POL@sturmname _\POL@sturm@N _}}%
-   \xintFor* ##1 in {\POL@isolzmult@indices}\do
-     {%
-      \csname POL@IfMultIsKnown##1\endcsname
-          {}% nothing to do
-          {\def\POL@isolzmult@index{##1}%
-           \POL@SturmIfZeroExactlyKnown{\POL@originalsturmname}{##1}%
-              \POL@isolzmult@loop@zero_isknown
-              \POL@isolzmult@loop@zero_isnotknown
-           \POL@isolzmult@loop@sharedbody
-          }%
-     }%
-   \ifnum\POL@sturmfinaldeg>\z@
+   \edef\POL@isolzmult@newGCDdegree{\PolDegree{\POL@sturmname _\POL@sturm@N _}}%
+   \let\POL@isolzmult@index\POL@isolz@NbOfRoots
+   \xintloop
+     % ATTENTION that this executes macros which also modifies \POL@sturmname!
+     % (but not \POL@sturm@N)
+     \POL@isolzmult@doone
+     \edef\POL@isolzmult@index{\the\numexpr\POL@isolzmult@index-\@ne}%
+   \if1\ifnum\POL@isolz@NbOfRoots@with_unknown_mult=\z@ 0\fi
+       \ifnum\POL@isolzmult@index=\z@ 0\fi 1%
+   \repeat
+   \let\POL@sturmname\POL@isolzmult@newsturmname
+   \if1\ifnum\POL@isolz@NbOfRoots@with_unknown_mult=\z@ 0\fi
+       % (if new GCD is constant, time to abort)
+       \ifnum\POL@isolzmult@newGCDdegree=\z@ 0\fi 1%
      \expandafter\POL@isolzmult@loop
    \fi
 }%
-\def\POL@isolzmult@loop@zero_isknown{%
+\def\POL@isolzmult@doone{%
+   \csname POL@IfMultIsKnown\POL@isolzmult@index\endcsname
+     {}% nothing to do
+     {\POL@SturmIfZeroExactlyKnown{\POL@originalsturmname}%
+                                  {\POL@isolzmult@index}%
+           \POL@isolzmult@loop@isknown
+           \POL@isolzmult@loop@isnotknown
+      \POL@isolzmult@loop@sharedbody
+     }%
+}%
+\def\POL@isolzmult@loop@isknown{%
    \xintifZero
-     {\Pol@Eval{\POL@sturmname _0_}%
+   % attention that \POL@eval requires a declared polynomial
+     {\POL@eval{\POL@isolzmult@lastGCD}%
                {\POL@xintexprGetVar{\POL@originalsturmname L_\POL@isolzmult@index}}}%
    {\let\POL@isolzmult@haszero\@ne}%
    {\let\POL@isolzmult@haszero\z@}%
 }%
-\def\POL@isolzmult@loop@zero_isnotknown{%
+\def\POL@isolzmult@loop@isnotknown{%
    \edef\POL@isolzmult@loop@A 
         {\POL@xintexprGetVar{\POL@originalsturmname L_\POL@isolzmult@index}}
    \edef\POL@isolzmult@loop@B
-        {\POL@xintexprGetVar{\POL@originalsturmname R_\POL@isolzmult@index}}
+        {\POL@xintexprGetVar{\POL@originalsturmname
+            R_\POL@isolzmult@index}}
+   % attention that \PolSetToNbOfZerosWithin sets \POL@sturmname to 2nd argument
    \PolSetToNbOfZerosWithin
       \POL@isolzmult@haszero  % nb of zeros A < x <= B, here 0 or 1
-      \POL@sturmname
+      \POL@isolzmult@newsturmname
       \POL@isolzmult@loop@A
       \POL@isolzmult@loop@B
 }%
@@ -928,10 +1402,6 @@
      \@namedef{POL@IfMultIsKnown\POL@isolzmult@index}{\xint_firstoftwo}%
      \edef\POL@isolz@NbOfRoots@with_unknown_mult
        {\the\numexpr\POL@isolz@NbOfRoots@with_unknown_mult-\@ne}%
-     \ifnum\POL@isolz@NbOfRoots@with_unknown_mult=\z@
-        \def\POL@sturmfinaldeg{0}% flag to force termination
-        \expandafter\expandafter\expandafter\xintBreakFor
-     \fi
    \fi
 }%
 
@@ -1039,7 +1509,7 @@
 }%
 % utility macro for a priori bound on root decimal exponent, via Float Rounding
 \def\POL@isolz@updateE #1e#2;%
-{\unless\ifnum#2<\POL@isolz@E\space\edef\POL@isolz@E{\the\numexpr#2+\@ne}\fi}%
+  {\unless\ifnum#2<\POL@isolz@E\space\edef\POL@isolz@E{\the\numexpr#2+\@ne}\fi}%
 \def\POL@isolz@getaprioribound{%
   \PolAssign{\POL@sturmname _0}\toarray\POL@arrayA
   \edef\POL@isolz@leading{\POL@arrayA{\POL@arrayA{0}}}%
@@ -1067,12 +1537,14 @@
 \def\POL@IsoRight@raw{\POL@IsoRight@Int/1[\POL@isolz@E]}%
 \def\POL@IsoLeft@raw {\POL@IsoLeft@Int/1[\POL@isolz@E]}%
 \def\POL@IsoRight@rawout{%
-    \ifnum\POL@IsoRightSign=\z@\expandafter\xintREZ\fi\POL@IsoRight@raw}%
+    \ifnum\POL@IsoRightSign=\z@\expandafter\xintREZ\fi\POL@IsoRight@raw
+}%
 \def\POL@IsoLeft@rawout{%
     \ifnum\POL@IsoRightSign=\z@
        \expandafter\xint_firstoftwo\else\expandafter\xint_secondoftwo
     \fi{\xintREZ\POL@IsoRight@raw}%
-       {\POL@IsoLeft@Int/1[\POL@isolz@E]}}%
+       {\POL@IsoLeft@Int/1[\POL@isolz@E]}%
+}%
 \def\POL@isolz@main {%
 % NOTE 2018/02/16. THIS WILL PRESUMABLY BE RE-ORGANIZED IN FUTURE TO DO
 % FIRST POSITIVE ROOTS THEN NEGATIVE ROOTS VIA CHANGE OF VARIABLE TO OPPOSITE.
@@ -1090,116 +1562,137 @@
     \edef\POL@IsoRightSV{\the\numexpr\POL@IsoRightSV+\@ne}%
 % subtlety here if original polynomial had multiplicities, but ok. I checked!
     \edef\POL@IsoRightSign % evaluated twice, but that's not so bad
-      {\xintiiOpp{\xintiiSgn{\Pol@Eval{\POL@sturmname _1}{0/1[0]}}}}%
+      {\xintiiOpp{\xintiiSgn{\POL@eval{\POL@sturmname _1}{0/1[0]}}}}%
   \fi
   \def\POL@IsoLeft@Int{-1}% -10^E isn't a root!
   \let\POL@IsoLeftSV  \POL@isolz@minusinf@SV
   \let\POL@IsoLeftSign\POL@isolz@minusinf@sign
+  % \POL@IsoRight@SV was modified if zero is a root
   \edef\POL@isolz@NbOfNegRoots{\the\numexpr\POL@IsoLeftSV-\POL@IsoRightSV}%
   \gdef\POL@isolz@IntervalIndex{0}%
-  \begingroup
-  \let\POL@IsoAtZeroSV\POL@IsoRightSV    % locally shifted if root at zero
-  \let\POL@IsoAtZeroSign\POL@IsoRightSign
+  \let\POL@isolz@@E\POL@isolz@E
   \ifnum\POL@isolz@NbOfNegRoots>\z@
-      \def\POL@IsoRight@Int{-1}%
-      \xintloop
-        \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}%
-        \POL@sturmchain@getSV@at\POL@IsoRight@raw
-        \let\POL@IsoRightSV  \POL@sturmchain@SV
-        \let\POL@IsoRightSign\POL@sturmchain@sign
-      % would an \ifx test be quicker? (to be checked)
-      \ifnum\POL@IsoRightSV=\POL@IsoLeftSV\space
-        % no roots in-between, sign and SV kept
-      \repeat
-      \def\POL@IsoLeft@Int{-10}%
-      \let\POL@@IsoRightSign\POL@IsoRightSign % zero possible
-      \let\POL@@IsoRightSV\POL@IsoRightSV
-      \xintloop
-        \edef\POL@IsoRight@Int{\the\numexpr\POL@IsoLeft@Int+\@ne}%
-% we could arguably do a more efficient dichotomy here
-        \POL@sturmchain@getSV@at\POL@IsoRight@raw
-        \let\POL@IsoRightSV  \POL@sturmchain@SV
-        \let\POL@IsoRightSign\POL@sturmchain@sign
-        \POL@isolz@check
-      \ifnum\POL@isolz@IntervalIndex=\POL@isolz@NbOfNegRoots\space
-        \expandafter\xintbreakloop
-      \fi
-        \let\POL@IsoLeft@Int\POL@IsoRight@Int
-        \let\POL@IsoLeftSign\POL@IsoRightSign
-        \let\POL@IsoLeftSV\POL@IsoRightSV
-      \ifnum\POL@IsoRight@Int < -\tw@
-      \repeat
-      \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfNegRoots\space
-        \def\POL@IsoRight@Int{-1}%
-        \let\POL@IsoRightSign\POL@@IsoRightSign
-        \let\POL@IsoRightSV\POL@@IsoRightSV
-        \POL@isolz@check
-        \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfNegRoots\space
-          \def\POL@IsoLeft@Int{-1}%
-          \let\POL@IsoLeftSign\POL@IsoRightSign
-          \let\POL@IsoLeftSV\POL@IsoRightSV
-          \def\POL@IsoRight@Int{0}%
-          \let\POL@IsoRightSV\POL@IsoAtZeroSV    % altered if 0 was a root
-          \let\POL@IsoRightSign\POL@IsoAtZeroSign% id.
-% this will recurse to locate roots with smaller decimal exponents
-          \POL@isolz@check % attention that this should not re-evaluate at 0
-       \fi
-     \fi
+% refactored at 0.7 to fix cases leading to an intervals with zero as end-point
+    \POL@isolz@findroots@neg
   \fi
-  \endgroup
+  \let\POL@isolz@E\POL@isolz@@E
   \def\POL@IsoLeft@Int{0}%
-  \let\POL@IsoLeftSV  \POL@IsoAtZeroSV
-  \let\POL@IsoLeftSign\POL@IsoAtZeroSign
+  \let\POL@IsoLeftSV  \POL@IsoAtZeroSV  % véritable SV en zéro
+  \let\POL@IsoLeftSign\POL@IsoAtZeroSign% véritable signe en zéro
   \ifnum\POL@IsoLeftSign=\z@
     \xdef\POL@isolz@IntervalIndex{\the\numexpr\POL@isolz@IntervalIndex+\@ne}%
-    \global\POL@isolz@nextwillneedrefinetrue
-  \else
-    \global\POL@isolz@nextwillneedrefinefalse
   \fi
   \let\POL@@IsoRightSV  \POL@isolz@plusinf@SV
   \let\POL@@IsoRightSign\POL@isolz@plusinf@sign % 10^E not a root!
   \edef\POL@isolz@NbOfPosRoots
     {\the\numexpr\POL@IsoLeftSV-\POL@@IsoRightSV}% attention @@
   \ifnum\POL@isolz@NbOfPosRoots>\z@
-      \def\POL@IsoRight@Int{1}%
-      \xintloop
-        \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}%
-        \POL@sturmchain@getSV@at\POL@IsoRight@raw
-        \let\POL@IsoRightSV  \POL@sturmchain@SV
-        \let\POL@IsoRightSign\POL@sturmchain@sign
-      \ifnum\POL@IsoRightSV=\POL@@IsoRightSV\space
-        \let\POL@@IsoRightSign\POL@IsoRightSign % root here possible!
-      \repeat
-      \unless\ifnum\POL@IsoRightSV=\POL@IsoLeftSV\space
-         \POL@isolz@check % will recurse inside groups if needed
-      \fi
-      \def\POL@IsoLeft@Int{1}%
-      \let\POL@IsoLeftSV\POL@IsoRightSV
-      \let\POL@IsoLeftSign\POL@IsoRightSign
-      \xintloop
-% we could arguably do a more efficient dichotomy here
-        \edef\POL@IsoRight@Int{\the\numexpr\POL@IsoLeft@Int+\@ne}%
-        \POL@sturmchain@getSV@at\POL@IsoRight@raw
-        \let\POL@IsoRightSV  \POL@sturmchain@SV
-        \let\POL@IsoRightSign\POL@sturmchain@sign
-        \POL@isolz@check
-        \let\POL@IsoLeft@Int\POL@IsoRight@Int
-        \let\POL@IsoLeftSign\POL@IsoRightSign
-        \let\POL@IsoLeftSV\POL@IsoRightSV
-      \ifnum\POL@isolz@IntervalIndex=\POL@isolz@NbOfRoots\space
-          \expandafter\xintbreakloop
-      \fi
-      \ifnum\POL@IsoLeft@Int < \xint_c_ix
-      \repeat
-      \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfRoots\space
-      % get now the last, rightmost, root (or roots)
-        \def\POL@IsoRight@Int{10}%
-        \let\POL@IsoRightSign\POL@@IsoRightSign
-        \let\POL@IsoRightSV\POL@@IsoRightSV
-        \POL@isolz@check
-      \fi
+    % always do that to avoid zero as end-point whether it is a root or not
+    \global\POL@isolz@nextwillneedrefinetrue
+    \POL@isolz@findroots@pos
+  \fi
+}%
+\def\POL@isolz@findroots@neg{%
+  \def\POL@IsoRight@Int{-1}%
+  \POL@isolz@findnextzeroboundeddecade@neg
+  \def\POL@IsoLeft@Int{-10}%
+  \let\POL@@IsoRightSign\POL@IsoRightSign % a zero there is possible
+  \let\POL@@IsoRightSV  \POL@IsoRightSV
+  % this will do possibly recursive \POL@isolz@check's
+  \POL@isolz@explorenexteightsubdecades@neg
+  \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfNegRoots\space
+    % above did not explore -2, -1 for this optimization (SV known at Right)
+    \def\POL@IsoRight@Int{-1}%
+    \let\POL@IsoRightSign\POL@@IsoRightSign
+    \let\POL@IsoRightSV  \POL@@IsoRightSV
+    \POL@isolz@check
+    \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfNegRoots\space
+      \def\POL@IsoLeft@Int{-1}%
+      \let\POL@IsoLeftSign\POL@@IsoRightSign
+      \let\POL@IsoLeftSV  \POL@@IsoRightSV
+      % I don't like being inside TeX conditionals
+      \expandafter\expandafter\expandafter\POL@isolz@findroots@neg
+    \fi
+  \fi
+}%
+\def\POL@isolz@findnextzeroboundeddecade@neg{%
+  \xintloop
+    \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}%
+    \POL@sturmchain@getSV@at\POL@IsoRight@raw
+    \let\POL@IsoRightSV  \POL@sturmchain@SV
+    \let\POL@IsoRightSign\POL@sturmchain@sign
+    % would an \ifx test be quicker? (to be checked)
+  \ifnum\POL@IsoRightSV=\POL@IsoLeftSV\space
+  % no roots in-between, iterate
+  \repeat
+}%
+\def\POL@isolz@explorenexteightsubdecades@neg{%
+  \xintloop
+    \edef\POL@IsoRight@Int{\the\numexpr\POL@IsoLeft@Int+\@ne}%
+     % we could arguably do a more efficient dichotomy here
+    \POL@sturmchain@getSV@at\POL@IsoRight@raw
+    \let\POL@IsoRightSV  \POL@sturmchain@SV
+    \let\POL@IsoRightSign\POL@sturmchain@sign
+    \POL@isolz@check % may recurse if multiple roots are to be found
+    \ifnum\POL@isolz@IntervalIndex=\POL@isolz@NbOfNegRoots\space
+      \expandafter\xintbreakloop
+    \fi
+    \let\POL@IsoLeft@Int\POL@IsoRight@Int
+    \let\POL@IsoLeftSign\POL@IsoRightSign
+    \let\POL@IsoLeftSV\POL@IsoRightSV
+  \ifnum\POL@IsoRight@Int < -\tw@
+  \repeat
+}%
+\def\POL@isolz@findroots@pos{%
+  % remark (2018/12/08), this needs some refactoring, I hardly understand
+  % the logic and it hides most into the recursion done by \POL@isolz@check
+  % It would probably make more sense to proceed like done for the negative
+  % but here finding the largest roots first.
+  \def\POL@IsoRight@Int{1}%
+  \POL@isolz@findnextzeroboundeddecade@pos
+  \unless\ifnum\POL@IsoRightSV=\POL@IsoLeftSV\space
+    % this actually explores the whole of some interval (0, 10^{e-1}]
+    % in a context where some roots are known to be in (10^{e-1}, 10^{e}]
+    % and none are larger
+    \POL@isolz@check % will recurse inside groups if needed with modified E
+  \fi
+  % we know get the roots in the last 9 decades from 10^{e-1} to 10^{e}
+  % we should arguably do a more efficient dichotomy here
+  \def\POL@IsoLeft@Int{1}%
+  \let\POL@IsoLeftSV\POL@IsoRightSV
+  \let\POL@IsoLeftSign\POL@IsoRightSign
+  \xintloop
+    \edef\POL@IsoRight@Int{\the\numexpr\POL@IsoLeft@Int+\@ne}%
+    \POL@sturmchain@getSV@at\POL@IsoRight@raw
+    \let\POL@IsoRightSV  \POL@sturmchain@SV
+    \let\POL@IsoRightSign\POL@sturmchain@sign
+    \POL@isolz@check % recurses in needed
+    \let\POL@IsoLeft@Int\POL@IsoRight@Int
+    \let\POL@IsoLeftSign\POL@IsoRightSign
+    \let\POL@IsoLeftSV\POL@IsoRightSV
+    \ifnum\POL@isolz@IntervalIndex=\POL@isolz@NbOfRoots\space
+      \expandafter\xintbreakloop
+    \fi
+    \ifnum\POL@IsoLeft@Int < \xint_c_ix
+  \repeat
+  \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfRoots\space
+  % get now the last, rightmost, root (or roots)
+    \def\POL@IsoRight@Int{10}%
+    \let\POL@IsoRightSign\POL@@IsoRightSign
+    \let\POL@IsoRightSV\POL@@IsoRightSV
+    \POL@isolz@check
   \fi
 }%
+\def\POL@isolz@findnextzeroboundeddecade@pos{%
+  \xintloop
+    \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}%
+    \POL@sturmchain@getSV@at\POL@IsoRight@raw
+    \let\POL@IsoRightSV  \POL@sturmchain@SV
+    \let\POL@IsoRightSign\POL@sturmchain@sign
+  \ifnum\POL@IsoRightSV=\POL@@IsoRightSV\space
+    \let\POL@@IsoRightSign\POL@IsoRightSign % root here possible!
+  \repeat
+}%
 \def\POL@isolz@check{% \POL@IsoRightSign must be ready for use here
 % \ifxintverbose
 %   \xintMessage{polexpr}{Info}%
@@ -1220,7 +1713,7 @@
      \ifPOL@isolz@nextwillneedrefine
        \expandafter\expandafter\expandafter\POL@isolz@refine
      \else
-       % \POL@IsoRightSign is zero iff root now exactly know
+       % \POL@IsoRightSign is zero iff root now exactly known
        \POL@refine@storeleftandright
        \ifnum\POL@IsoRightSign=\z@
         \global\POL@isolz@nextwillneedrefinetrue
@@ -1275,7 +1768,7 @@
     \edef\POL@IsoLeft@Int {\xintDSL{\POL@IsoLeft@Int}}%
     \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}%
     \edef\POL@IsoRightSign
-        {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+        {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
   \ifnum\POL@IsoRightSign=\POL@@IsoRightSign\space
   \repeat
   % now second root has been separated from the one at left end point
@@ -1294,7 +1787,7 @@
   \else
     \edef\POL@IsoRight@Int{\xintDec{\POL@@IsoRight@Int}}%
     \edef\POL@IsoRightSign
-      {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+      {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
     \ifnum\POL@IsoRightSign=\POL@@IsoRightSign\space
       \POL@refine@doonce % we need to locate in interval (1, 9) in local scale
     \else
@@ -1319,17 +1812,17 @@
   \let\POL@@IsoRightSign\POL@IsoRightSign
   \edef\POL@IsoRight@Int{\xintiiAdd{4}{\POL@IsoLeft@Int}}% 5
   \edef\POL@IsoRightSign
-      {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+      {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
   \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space
     \let\POL@IsoLeft@Int\POL@IsoRight@Int % 5
     \edef\POL@IsoRight@Int{\xintiiAdd{2}{\POL@IsoLeft@Int}}%
     \edef\POL@IsoRightSign
-        {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+        {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
     \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space
       \let\POL@IsoLeft@Int\POL@IsoRight@Int % 7
       \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}%
       \edef\POL@IsoRightSign
-          {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+          {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
       \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space
         \let\POL@IsoLeft@Int\POL@IsoRight@Int % 8
         \let\POL@IsoRight@Int\POL@@IsoRight@Int % 9
@@ -1343,7 +1836,7 @@
       \let\POL@@IsoRight@Int\POL@IsoRight@Int % 7
       \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% 6
       \edef\POL@IsoRightSign
-          {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+          {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
       \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space
           \let\POL@IsoLeft@Int\POL@IsoRight@Int   % 6
           \let\POL@IsoRight@Int\POL@@IsoRight@Int % 7
@@ -1358,12 +1851,12 @@
     \let\POL@@IsoRight@Int\POL@IsoRight@Int % 5
     \edef\POL@IsoRight@Int{\xintiiAdd{2}{\POL@IsoLeft@Int}}%
     \edef\POL@IsoRightSign
-        {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+        {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
     \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space
       \let\POL@IsoLeft@Int\POL@IsoRight@Int % 3
       \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% 4
       \edef\POL@IsoRightSign
-          {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+          {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
       \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space
           \let\POL@IsoLeft@Int\POL@IsoRight@Int   % 4
           \let\POL@IsoRight@Int\POL@@IsoRight@Int % 5
@@ -1377,7 +1870,7 @@
       \let\POL@@IsoRight@Int\POL@IsoRight@Int % 3
       \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% 2
       \edef\POL@IsoRightSign
-          {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+          {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
       \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space
           \let\POL@IsoLeft@Int\POL@IsoRight@Int   % 2
           \let\POL@IsoRight@Int\POL@@IsoRight@Int % 3
@@ -1395,35 +1888,43 @@
     \xdef\csname POL_ZeroInt\POL@sturmname
                  R\POL@isolz@IntervalIndex\endcsname
          {\PolDecToString{\POL@IsoRight@rawout}}%
-    \begingroup\xintglobaldefstrue
-    \xintdefvar\POL@sturmname
-               L_\POL@isolz@IntervalIndex:=qfrac(\POL@IsoLeft@rawout);%
-    \xintdefvar\POL@sturmname
-               R_\POL@isolz@IntervalIndex:=qfrac(\POL@IsoRight@rawout);%
-    \endgroup
-    % added at 0.6+
+    % added at 0.6
     \ifnum\POL@IsoRightSign=\z@
       \global
       \expandafter
       \let\csname POL_ZeroIsKnown\POL@sturmname\POL@isolz@IntervalIndex\endcsname
           \xint_stop_atfirstoftwo
     \fi
+    \begingroup\xintglobaldefstrue
+    \xintdefvar
+      \POL@sturmname L_\POL@isolz@IntervalIndex,%
+      \POL@sturmname R_\POL@isolz@IntervalIndex,%
+      % added at 0.7
+      \POL@sturmname Z_\POL@isolz@IntervalIndex _isknown
+        := qfrac(\POL@IsoLeft@rawout),%
+           qfrac(\POL@IsoRight@rawout),%
+           \ifnum\POL@IsoRightSign=\z@ 1\else 0\fi;%
+    \endgroup
 }%
 
 
 %% \PolRefineInterval
 \def\POL@xintexprGetVar#1{\expandafter\expandafter\expandafter
     \XINT_expr_unlock\csname XINT_expr_var_#1\endcsname}%
-\def\POL@set@IsoLeft@rawin{%
+% attention, also used by \POL@findrat@loop@a
+\def\POL@get@IsoLeft@rawin{%
     \edef\POL@IsoLeft@rawin
      {\POL@xintexprGetVar{\POL@sturmname L_\POL@isolz@IntervalIndex}}%
 }%
-\def\POL@set@IsoRight@rawin{%
+% attention, also used by \POL@findrat@loop@a
+\def\POL@get@IsoRight@rawin{%
     \edef\POL@IsoRight@rawin
      {\POL@xintexprGetVar{\POL@sturmname R_\POL@isolz@IntervalIndex}}%
 }%
-\def\POL@set@IsoLeft@Int #1/1[#2]{%
-    \edef\POL@IsoLeft@Int{\xintDSH{\POL@isolz@E-#2}{#1}}%
+% attention, also used by \POL@findrat@loop@a
+\def\POL@get@Int@aux #1/1[#2]#3#4{\edef#3{\xintDSH{#4-#2}{#1}}}%
+\def\POL@get@IsoLeft@Int{%
+    \expandafter\POL@get@Int@aux\POL@IsoLeft@rawin\POL@IsoLeft@Int\POL@isolz@E
 }%
 \newcommand\PolRefineInterval{\@ifstar\POL@srefine@start\POL@refine@start}%
 \newcommand\POL@refine@start[3][1]{%
@@ -1441,19 +1942,19 @@
      \POL@refine@main}%
 }%
 \def\POL@refine@sharedbody#1{%
-  \POL@set@IsoLeft@rawin
+  \POL@get@IsoLeft@rawin
   \edef\POL@IsoLeftSign
-      {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoLeft@rawin}}}%
+      {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoLeft@rawin}}}%
   \ifnum\POL@IsoLeftSign=\z@
   % do nothing if that interval was already a singleton
   \else
   % else both end-points are not roots and there is a single one in-between
-    \POL@set@IsoRight@rawin
+    \POL@get@IsoRight@rawin
     \edef\POL@IsoRightSign{\the\numexpr-\POL@IsoLeftSign}%
     \edef\POL@isolz@E{\expandafter\POL@refine@getE
          % je pense que le xintrez ici est superflu
          \romannumeral0\xintrez{\xintSub{\POL@IsoRight@rawin}{\POL@IsoLeft@rawin}}}%
-    \expandafter\POL@set@IsoLeft@Int\POL@IsoLeft@rawin
+    \POL@get@IsoLeft@Int
     \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}%
     #1%
     \POL@refine@storeleftandright % \POL@IsoRightSign not zero
@@ -1478,7 +1979,7 @@
   \let\POL@@IsoRightSign\POL@IsoRightSign
   \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}%
   \edef\POL@IsoRightSign
-      {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+      {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
   \ifnum\POL@IsoRightSign=\z@
    \let\POL@IsoLeft@Int\POL@IsoRight@Int % root at 1
    \def\POL@IsoLeftSign{0}%
@@ -1491,7 +1992,7 @@
     \let\POL@IsoLeft@Int\POL@IsoRight@Int
     \edef\POL@IsoRight@Int{\xintDec{\POL@@IsoRight@Int}}%
     \edef\POL@IsoRightSign
-      {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+      {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
     \ifnum\POL@IsoRightSign=\z@
      \let\POL@IsoLeft@Int\POL@IsoRight@Int % root at 9
      \def\POL@IsoLeftSign{0}%
@@ -1525,11 +2026,16 @@
 
 \newcommand\PolEnsureIntervalLengths[2]{% #1 = Sturm chain name,
    % localize roots in intervals of length at most 10^{#2}
-   \POL@count\z@
-   % \POL@count used by \POL@sturmchain@getSV@at but latter not used
    \edef\POL@sturmname{#1}%
    \edef\POL@ensure@targetE{\the\numexpr#2}%
-   \edef\POL@nbofroots{\csname POL_ZeroInt\POL@sturmname L\endcsname 0}%
+   \edef\POL@nbofroots{\csname POL_ZeroInt\POL@sturmname L0\endcsname}%
+   \ifnum\POL@nbofroots>\z@
+      \expandafter\POL@ensureintervallengths
+   \fi
+}%
+\def\POL@ensureintervallengths{%
+   \POL@count\z@
+   % \POL@count used by \POL@sturmchain@getSV@at but latter not used
    \xintloop
      \advance\POL@count\@ne
      \edef\POL@isolz@IntervalIndex{\the\POL@count}%
@@ -1543,20 +2049,27 @@
    \edef\POL@sturmname{#1}%
    \edef\POL@ensure@targetE{\the\numexpr#3}%
    \edef\POL@isolz@IntervalIndex{\the\numexpr#2}%
-   \POL@ensure@one
+% peut-être autoriser -1, -2, ... ?
+   \ifnum\POL@isolz@IntervalIndex>\z@
+% 0.7, add this safeguard but attention means this structure must be in place
+   \ifnum\csname POL_ZeroInt\POL@sturmname L0\endcsname>\z@
+% je ne fais pas les \expandafter mais je préfèrerai ne pas être à l'intérieur
+     \POL@ensure@one
+   \fi
+   \fi
 }%
 \def\POL@ensure@one{%
-    \POL@set@IsoLeft@rawin
-    \POL@set@IsoRight@rawin
+    \POL@get@IsoLeft@rawin
+    \POL@get@IsoRight@rawin
     \edef\POL@ensure@delta{\xintREZ{\xintSub{\POL@IsoRight@rawin}{\POL@IsoLeft@rawin}}}%
     \xintiiifZero{\POL@ensure@delta}
       {}
       {\edef\POL@isolz@E{\expandafter\POL@refine@getE\POL@ensure@delta}%
-       \expandafter\POL@set@IsoLeft@Int\POL@IsoLeft@rawin
+       \POL@get@IsoLeft@Int
        \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}%
        \ifnum\POL@isolz@E>\POL@ensure@targetE\space
          \edef\POL@IsoLeftSign
-          {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoLeft@raw}}}%
+          {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoLeft@raw}}}%
           % at start left and right are not roots, and values of opposite signs
          % \edef\POL@IsoRightSign{\the\numexpr-\POL@IsoLeftSign}%
          \xintloop
@@ -1576,7 +2089,7 @@
     \xintloop
       \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}%
       \edef\POL@IsoRightSign
-          {\xintiiSgn{\Pol@Eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
+          {\xintiiSgn{\POL@eval{\POL@sturmname _0}{\POL@IsoRight@raw}}}%
           % if we have found a zero at right boundary the \ifnum test will fail
           % and we exit the loop
           % else we exit the loop if sign at right boundary is opposite of
@@ -1595,54 +2108,119 @@
 
 
 \catcode`_ 8
-\newcommand\PolPrintIntervals[2][Z]{%
-   \POL@count \@nameuse{POL_ZeroInt#2L}{0}
-   \ifnum\POL@count=\z@
-%     No real roots.\par
+\newcommand\PolPrintIntervals
+           {\@ifstar{\PolPrintIntervals@@}{\PolPrintIntervals@}}%
+\newcommand\PolPrintIntervals@@{%
+ \begingroup
+   \def\POL@AfterPrintIntervals{\endgroup}%
+   \def\arraystretch{2}%
+   \let\PolPrintIntervalsPrintExactZero\POL@@PrintIntervalsPrintExactZero
+   \let\PolPrintIntervalsUnknownRoot\POL@@PrintIntervalsUnknownRoot
+   \let\PolPrintIntervalsKnownRoot\POL@@PrintIntervalsKnownRoot
+   \def\PolPrintIntervalsBeginEnv{\[\begin{array}{cl}}%\]
+   \def\PolPrintIntervalsEndEnv{\end{array}\]}%
+ \PolPrintIntervals@
+}%
+\newcommand\PolPrintIntervals@[2][Z]{\POL@PrintIntervals{#1}{#2}}%
+\newcommand\POL@PrintIntervals[2]{%
+   \def\PolPrintIntervalsTheSturmName{#2}%
+   \def\PolPrintIntervalsTheVar{#1}%
+   \ifnum\@nameuse{POL_ZeroInt#2L}{0}=\z@
+     \PolPrintIntervalsNoRealRoots
    \else
-%     There are \the\POL@count\space distinct real roots:\par
-     \[\count@\POL@count
-       \global\POL@count\@ne
-       \begin{array}{rcccl}
-        \xintloop
-        \POL@SturmIfZeroExactlyKnown{#2}\POL@count
-        {% exact root
-         &&
-         #1_{\the\POL@count}&=&
-         \POL@printintervals@prepare{#2R}%
-         \PolPrintIntervalsPrintExactZero
-        }%
-        {% interval with root in its strict interior
-         \POL@printintervals@prepare{#2L}%
-         \PolPrintIntervalsPrintLeftEndPoint&<&
-         #1_{\the\POL@count}&<&
-         \POL@printintervals@prepare{#2R}%
-         \PolPrintIntervalsPrintRightEndPoint
-        }%
-        \global\advance\POL@count\@ne
-        \unless\ifnum\POL@count>\count@
-        \\%
-        \repeat
-     \end{array}\]
+     \gdef\PolPrintIntervalsTheIndex{1}%
+     \POL@PrintIntervals@DoDefs
+     \begingroup\edef\POL@tmp{\endgroup
+       \unexpanded\expandafter{\PolPrintIntervalsBeginEnv}%
+       \unexpanded\expandafter{\POL@PrintIntervals@Loop}%
+       \unexpanded\expandafter{\PolPrintIntervalsEndEnv}%
+     }\POL@tmp
    \fi
+   \POL@AfterPrintIntervals
+}%
+\let\POL@AfterPrintIntervals\@empty
+\newcommand\PolPrintIntervalsNoRealRoots{}%
+\newcommand\PolPrintIntervalsBeginEnv{\[\begin{array}{rcccl}}%
+\newcommand\PolPrintIntervalsEndEnv{\end{array}\]}%
+\newcommand\PolPrintIntervalsKnownRoot{%
+  &&\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}%
+  &=&\PolPrintIntervalsPrintExactZero
+}%
+\newcommand\PolPrintIntervalsUnknownRoot{%
+  \PolPrintIntervalsPrintLeftEndPoint&<&%
+  \PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}&<&%
+  \PolPrintIntervalsPrintRightEndPoint
+}%
+\newcommand\PolPrintIntervalsPrintExactZero    {\PolPrintIntervalsTheLeftEndPoint}%
+\newcommand\PolPrintIntervalsPrintLeftEndPoint {\PolPrintIntervalsTheLeftEndPoint}%
+\newcommand\PolPrintIntervalsPrintRightEndPoint{\PolPrintIntervalsTheRightEndPoint}%
+\newcommand\PolPrintIntervalsPrintMultiplicity{(\mbox{mult. }\PolPrintIntervalsTheMultiplicity)}%
+%
+\newcommand\POL@@PrintIntervalsKnownRoot{%
+   \PolPrintIntervalsPrintMultiplicity&%
+   \PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}=%
+   \PolPrintIntervalsPrintExactZero
+}%
+\newcommand\POL@@PrintIntervalsPrintExactZero{%
+   \displaystyle
+   \xintSignedFrac{\PolPrintIntervalsTheLeftEndPoint}%
+}%
+\newcommand\POL@@PrintIntervalsUnknownRoot{%
+   \PolPrintIntervalsPrintMultiplicity&%
+   \xintifSgn{\PolPrintIntervalsTheLeftEndPoint}%
+     {\xintifSgn{\PolPrintIntervalsTheRightEndPoint}
+       {\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}=%
+        \PolPrintIntervalsPrintRightEndPoint\dots}%
+      {0>\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}>%
+        \PolPrintIntervalsPrintLeftEndPoint}%
+      {\PolErrorThisShouldNotHappenPleaseReportToAuthorA}}%
+    {\xintifSgn{\PolPrintIntervalsTheRightEndPoint}
+      {\PolErrorThisShouldNotHappenPleaseReportToAuthorB}%
+      {\PolErrorThisShouldNotHappenPleaseReportToAuthorC}%
+      {0<\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}<%
+        \PolPrintIntervalsPrintRightEndPoint}}%
+    {\xintifSgn{\PolPrintIntervalsTheRightEndPoint}
+      {\PolErrorThisShouldNotHappenPleaseReportToAuthorD}%
+      {\PolErrorThisShouldNotHappenPleaseReportToAuthorE}%
+      {\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}=%
+        \PolPrintIntervalsPrintLeftEndPoint\dots}}%
 }%
+%
 \catcode`_ 11
-\newcommand\PolPrintIntervalsPrintExactZero {\PolPrintIntervalsTheEndPoint}%
-\newcommand\PolPrintIntervalsPrintLeftEndPoint {\PolPrintIntervalsTheEndPoint}%
-\newcommand\PolPrintIntervalsPrintRightEndPoint{\PolPrintIntervalsTheEndPoint}%
-\def\POL@printintervals@prepare#1{%
-    \edef\PolPrintIntervalsTheIndex{\the\POL@count}%
-    \edef\PolPrintIntervalsTheEndPoint{\@nameuse{POL_ZeroInt#1}\POL@count}%
-    \xintiiifSgn{\POL@xintexprGetVar{#1_\PolPrintIntervalsTheIndex}}
-       {\let\PolIfEndPointIsPositive\xint_secondoftwo
-        \let\PolIfEndPointIsNegative\xint_firstoftwo
-        \let\PolIfEndPointIsZero\xint_secondoftwo}
-       {\let\PolIfEndPointIsPositive\xint_secondoftwo
-        \let\PolIfEndPointIsNegative\xint_secondoftwo
-        \let\PolIfEndPointIsZero\xint_firstoftwo}
-       {\let\PolIfEndPointIsPositive\xint_firstoftwo
-        \let\PolIfEndPointIsNegative\xint_secondoftwo
-        \let\PolIfEndPointIsZero\xint_secondoftwo}%
+\def\POL@PrintIntervals@Loop{%
+  \POL@SturmIfZeroExactlyKnown\PolPrintIntervalsTheSturmName
+                              \PolPrintIntervalsTheIndex
+    \PolPrintIntervalsKnownRoot
+    \PolPrintIntervalsUnknownRoot
+  \xdef\PolPrintIntervalsTheIndex{\the\numexpr\PolPrintIntervalsTheIndex+\@ne}%
+  \unless\ifnum\PolPrintIntervalsTheIndex>
+               \@nameuse{POL_ZeroInt\PolPrintIntervalsTheSturmName L0}
+     \POL@PrintIntervals@DoDefs
+     \xint_afterfi{\\\POL@PrintIntervals@Loop}%
+  \fi
+}%
+\def\POL@PrintIntervals@DoDefs{%
+  \xdef\PolPrintIntervalsTheLeftEndPoint{%
+     \csname POL_ZeroInt%
+             \PolPrintIntervalsTheSturmName L\PolPrintIntervalsTheIndex
+     \endcsname
+  }%
+  \xdef\PolPrintIntervalsTheRightEndPoint{%
+     \csname POL_ZeroInt%
+             \PolPrintIntervalsTheSturmName R\PolPrintIntervalsTheIndex
+     \endcsname
+  }%
+  \xdef\PolPrintIntervalsTheMultiplicity{%
+     \ifcsname POL_ZeroMult%
+             \PolPrintIntervalsTheSturmName\PolPrintIntervalsTheIndex
+     \endcsname
+     \csname POL_ZeroMult%
+             \PolPrintIntervalsTheSturmName\PolPrintIntervalsTheIndex
+     \endcsname
+     \else
+     ?% or use 0 ?
+     \fi
+  }%
 }%
 
 
@@ -1650,18 +2228,38 @@
     \romannumeral0\csname POL_ZeroIsKnown#1\endcsname{#2}%
 }%
 \newcommand\POL@SturmIfZeroExactlyKnown[2]{% #1 = sturmname, #2=index
-    \romannumeral0\csname POL_ZeroIsKnown#1\the\numexpr#2\relax\endcsname
+    \romannumeral0\csname POL_ZeroIsKnown#1\the\numexpr#2\endcsname
 }%
 \newcommand\PolSturmIsolatedZeroMultiplicity[2]{%
     \romannumeral`^^@\csname POL_ZeroMult#1\endcsname{#2}%
 }%
 \newcommand\PolSturmIsolatedZeroLeft[2]{%
-    \romannumeral`^^@\csname POL_ZeroInt#1L\endcsname{#2}}%
+    \romannumeral`^^@\csname POL_ZeroInt#1L\endcsname{#2}%
+}%
 \newcommand\PolSturmIsolatedZeroRight[2]{%
-    \romannumeral`^^@\csname POL_ZeroInt#1R\endcsname{#2}}%
+    \romannumeral`^^@\csname POL_ZeroInt#1R\endcsname{#2}%
+}%
 \newcommand\PolSturmNbOfIsolatedZeros[1]{%
     \romannumeral`^^@\csname POL_ZeroInt#1L0\endcsname
 }%
+\newcommand\PolSturmRationalRoot[2]{%
+    \romannumeral`^^@\csname POL_ZeroInt#1L%
+                     \csname POL_RRIndex#1\endcsname{#2}\endcsname
+}%
+\newcommand\PolSturmRationalRootIndex[2]{%
+    \romannumeral`^^@\csname POL_RRIndex#1\endcsname{#2}%
+}%
+\newcommand\PolSturmRationalRootMultiplicity[2]{%
+    \romannumeral`^^@\csname POL_ZeroMult#1%
+                     \csname POL_RRIndex#1\endcsname{#2}\endcsname
+}%
+\newcommand\PolSturmNbOfRationalRoots[1]{%
+    \romannumeral`^^@\csname POL_RRIndex#10\endcsname
+}%
+\newcommand\PolSturmNbOfRationalRootsWithMultiplicities[1]{%
+% means the \POL@norr must not have been changed in-between...
+    \the\numexpr\PolDegree{#1}-\PolDegree{#1\POL@norr}\relax
+}%
 
 
 \let\PolDecToString\xintDecToString
@@ -2023,18 +2621,18 @@
 
 
 %% EXPANDABLE MACROS
-\def\Pol@Eval@fork#1\At#2#3\krof{#2}%
-\newcommand\PolEval[3]{\romannumeral`^^@\Pol@Eval@fork
+\def\POL@eval@fork#1\At#2#3\krof{#2}%
+\newcommand\PolEval[3]{\romannumeral`^^@\POL@eval@fork
      #2\PolEvalAt
      \At\PolEvalAtExpr\krof {#1}{#3}%
 }%
 \newcommand\PolEvalAt[2]
     {\xintpraw{\csname XINT_expr_userfunc_#1\endcsname{#2}}}%
-\newcommand\Pol@Eval[2]
+\newcommand\POL@eval[2]
     {\csname XINT_expr_userfunc_#1\endcsname{#2}}%
 \newcommand\PolEvalAtExpr[2]{\xinttheexpr #1(#2)\relax}%
 %
-\newcommand\PolEvalReduced[3]{\romannumeral`^^@\Pol@Eval@fork
+\newcommand\PolEvalReduced[3]{\romannumeral`^^@\POL@eval@fork
      #2\PolEvalReducedAt
      \At\PolEvalReducedAtExpr\krof {#1}{#3}%
 }%
@@ -2047,7 +2645,7 @@
     {\xintIrr{\romannumeral`^^@\xintthebareeval#1(#2)\relax}[0]}%
 }%
 %
-\newcommand\PolFloatEval[3]{\romannumeral`^^@\Pol@Eval@fork
+\newcommand\PolFloatEval[3]{\romannumeral`^^@\POL@eval@fork
      #2\PolFloatEvalAt
      \At\PolFloatEvalAtExpr\krof {#1}{#3}%
 }%
@@ -2056,37 +2654,48 @@
 \newcommand\PolFloatEvalAtExpr[2]{\xintthefloatexpr #1(#2)\relax}%
 
 
-\newcommand\PolSturmMultiplicity[3]{\romannumeral`^^@\Pol@Eval@fork
-     #2\PolSturmMultiplicityAt
-     \At\PolSturmMultiplicityAtExpr\krof {#1}{#3}%
+\newcommand\PolSturmIntervalIndex[3]{\the\numexpr\POL@eval@fork
+     #2\PolSturmIntervalIndexAt
+     \At\PolSturmIntervalIndexAtExpr\krof {#1}{#3}%
 }%
-\newcommand\PolSturmMultiplicityAtExpr[2]
-    {\PolSturmMultiplicityAt{#1}{\xinttheexpr#2\relax}}%
-\newcommand\PolSturmMultiplicityAt[2]
-    {\expandafter\POL@sturm@mult@at\romannumeral`^^@#2!{#1}}%
-\def\POL@sturm@mult@at#1!#2%
+\newcommand\PolSturmIntervalIndexAtExpr[2]
+    {\PolSturmIntervalIndexAt{#1}{\xinttheexpr#2\relax}}%
+\newcommand\PolSturmIntervalIndexAt[2]
+    {\expandafter\POL@sturm@index@at\romannumeral`^^@#2!{#1}\xint_bye\relax}%
+\def\POL@sturm@index@at#1!#2%
 {%
-    \xintifZero{\Pol@Eval{#2_0}{#1}}%
-    {\POL@sturm@mult@at@iloop 1!{#2}{#1}}% we have a zero
-    0% not a zero
+    \expandafter\POL@sturm@index@at@iloop 
+    \romannumeral`^^@\PolSturmNbOfIsolatedZeros{#2}!{#2}{#1}%
 }%
-\def\POL@sturm@mult@at@iloop #1!#2#3%
+% implementation is sub-optimal as it should use some kind of binary tree
+% search rather than comparing to the intervals from right to left as here
+\def\POL@sturm@index@at@iloop #1!%
+{%
+    \ifnum #1=\z@ 0\expandafter\xint_bye\fi
+    \POL@sturm@index@at@iloop@a #1!%
+}%
+\def\POL@sturm@index@at@iloop@a #1!#2#3%
 {% #1 = index, #2 = sturmname, #3 value
-    \PolSturmIfZeroExactlyKnown{#2}{#1}%
-    {\xintifEq{\POL@xintexprGetVar{#2L_#1}}{#3}%
-     {\PolSturmIsolatedZeroMultiplicity{#2}{#1}}%
-% catcode of ! is 11 in polexpr.sty
-     {\expandafter\POL@sturm@mult@at@iloop\the\numexpr#1+\@ne !{#2}{#3}}%
+    \PolSturmIfZeroExactlyKnown{#2}{#1}
+    {\xintifCmp{#3}{\POL@xintexprGetVar{#2L_#1}}%
+       {}%
+       {#1\xint_bye}%
+       {0\xint_bye}%
     }%
-    {\xintifLt{#3}{\POL@xintexprGetVar{#2R_#1}}%
-     {\PolSturmIsolatedZeroMultiplicity{#2}{#1}}%
-     {\expandafter\POL@sturm@mult@at@iloop\the\numexpr#1+\@ne !{#2}{#3}}%
+    {\xintifGt{#3}{\POL@xintexprGetVar{#2L_#1}}%
+       {\xintifLt{#3}{\POL@xintexprGetVar{#2R_#1}}%
+        {#1\xint_bye}%
+        {0\xint_bye}%
+       }%
+       {}%
     }%
+    % catcode of ! is 11 in polexpr.sty
+    \expandafter\POL@sturm@index@at@iloop\the\numexpr#1-\@ne !{#2}{#3}%
 }%
 
 
-\def\Pol@LessThanOrEqualTo@fork#1\LessThanOrEqualTo#2#3\krof{#2}%
-\newcommand\PolSturmNbOfRootsOf[3]{\romannumeral`^^@\Pol@LessThanOrEqualTo@fork
+\def\POL@leq@fork#1\LessThanOrEqualTo#2#3\krof{#2}%
+\newcommand\PolSturmNbOfRootsOf[3]{\romannumeral`^^@\POL@leq@fork
      #2\PolNbOfRootsLessThanOrEqualTo
      \LessThanOrEqualTo\PolNbOfRootsLessThanOrEqualToExpr\krof {#1}{#3}%
 }%
@@ -2106,7 +2715,7 @@
 \def\POL@nbofrootsleq@prep#1!#2%
 {%
     \expandafter\POL@nbofrootsleq@iloop\expandafter 1\expandafter !%
-    \romannumeral0\xintsgn{\Pol@Eval{#2_0}{#1}}!%
+    \romannumeral0\xintsgn{\POL@eval{#2_0}{#1}}!%
     #1!{#2}%
 }%
 \def\POL@nbofrootsleq@iloop#1!#2!#3!#4%
@@ -2120,7 +2729,7 @@
       % the test \xintifLt will be negative
       {\xintifLt{#3}{\POL@xintexprGetVar{#4R_#1}}%
          {\POL@nbofrootsleq@return
-            #1\ifnum#2=\xintSgn{\Pol@Eval{#4_0}{\POL@xintexprGetVar{#4L_#1}}}
+            #1\ifnum#2=\xintSgn{\POL@eval{#4_0}{\POL@xintexprGetVar{#4L_#1}}}
               -\@ne\fi !%
          }%
          {\ifnum#1=\PolSturmNbOfIsolatedZeros{#4}
@@ -2135,9 +2744,8 @@
     \the\numexpr\@ne+#1!#2!#3!#4{#1}%
 
 
-\def\Pol@LessThanOrEqualTo@fork#1\LessThanOrEqualTo#2#3\krof{#2}%
 \newcommand\PolSturmNbWithMultOfRootsOf[3]
-{\the\numexpr0\Pol@LessThanOrEqualTo@fork
+{\the\numexpr0\POL@leq@fork
      #2\PolNbWithMultOfRootsLessThanOrEqualTo
      \LessThanOrEqualTo\PolNbWithMultOfRootsLessThanOrEqualToExpr\krof {#1}{#3}%
 }%
@@ -2158,7 +2766,7 @@
 \def\POL@nbwmofrootsleq@prep#1!#2%
 {%
     \expandafter\POL@nbwmofrootsleq@iloop\expandafter 1\expandafter !%
-    \romannumeral0\xintsgn{\Pol@Eval{#2_0}{#1}}!%
+    \romannumeral0\xintsgn{\POL@eval{#2_0}{#1}}!%
     #1!{#2}%
 }%
 \def\POL@nbwmofrootsleq@iloop#1!#2!#3!#4%
@@ -2174,7 +2782,7 @@
       {\xintifLt{#3}{\POL@xintexprGetVar{#4R_#1}}%
          {\POL@nbwmofrootsleq@return
             \unless
-            \ifnum#2=\xintSgn{\Pol@Eval{#4_0}{\POL@xintexprGetVar{#4L_#1}}}
+            \ifnum#2=\xintSgn{\POL@eval{#4_0}{\POL@xintexprGetVar{#4L_#1}}}
             +\PolSturmIsolatedZeroMultiplicity{#4}{#1}\fi !%
          }%
          {+\PolSturmIsolatedZeroMultiplicity{#4}{#1}%