polexpr (10jan22)

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

5 files changed, 1004 insertions, 112 deletions

diff --git a/Master/texmf-dist/tex/generic/polexpr/polexpr-examples.tex b/Master/texmf-dist/tex/generic/polexpr/polexpr-examples.tex
new file mode 100644
index 00000000000..fa2f330513b
--- /dev/null
+++ b/Master/texmf-dist/tex/generic/polexpr/polexpr-examples.tex
@@ -0,0 +1,881 @@
+% -*- sentence-end-double-space: t -*-
+\documentclass[a4paper,svgnames,dvipsnames,dvipdfmx]{article}
+\usepackage{geometry}
+\usepackage{shortvrb}
+\usepackage{xcolor}
+\usepackage{graphicx}
+\usepackage{polexpr}
+\usepackage{hyperref}
+\usepackage{bookmark}
+% \usepackage{amsmath}
+\usepackage{framed}
+\usepackage{newtxtext,newtxmath}
+\title{\pkg{polexpr} root localization examples}
+\author{Jean-François Burnol}
+\date{To access the reference documentation:\\
+% faut-il vraiment le ./ ?
+      \texttt{texdoc} \href{run:./polexpr.html}{polexpr.html}}
+%\usepackage{parskip}
+\MakeShortVerb{\|}
+
+% Note 25 juin 2021 et 8 janvier 2022 pour polexpr-examples.tex
+%
+% Ceci est extrait quasi verbatim de xint.dtx
+%
+% Le code est pour LaTeX, pas pour Plain TeX.
+% Prévu pour latex+dvipdfmx, et j'ai supprimé ajouts spécifiques pour
+% pdflatex/xetex/lualatex en ce qui concerne les couleurs et polices.
+%
+% Regarding latex+dvipdfmx, a document with many more usage
+% of everbatim* could run into color stack overflow problems
+% Refer to xint.dtx for how this problem is avoided there
+% via direct usage of \special rather than \color in the
+% \everbatimxprehook
+
+% verbatim macros and environments
+% ================================
+%
+% June 2013, then October 2014.
+% -----------------------------
+%
+\makeatletter
+\catcode`_ 11
+% some of my verbatim environments do not make the space active (\lverb e.g.). Then
+% \do@noligs must be modified, \char`#1 must be followed by a space token, else,
+% the `#1 expansion will swallow one space.
+\def\do@noligs #1{%
+  \catcode`#1\active
+  \begingroup
+     \lccode`~`#1\relax
+     \lowercase{%
+  \endgroup\def~{\leavevmode\kern\z@\char`#1 }}%
+}
+% \lowast
+% Pas forcément adapté à toutes les polices
+\def\lowast{\raisebox{-.25\height}{*}}
+\catcode`* 13
+\def\makestarlowast {\let*\lowast\catcode`\*\active}%
+\catcode`* 12
+
+
+%--- straight quotes, added (finally...) Nov 2, 2014
+%--- obsolete with use of newtxtt 1.05, late 2014
+\begingroup\makeatletter
+  \catcode`\'\active
+  \catcode`\`\active
+\@firstofone {\endgroup
+  \def\makequotesstraight{% assumes textcomp package
+% à propos textcomp est automatique avec pdflatex depuis Février 2020
+     \let`\textasciigrave
+     \let'\textquotesingle
+     \catcode39\active
+     \catcode96\active }%
+}
+
+% \verb
+% =====
+% Initially, June 2013, then Sep 9, 2014, and Oct 9-12 2014
+%
+% pour les short verb |...|
+
+\def\MicroFont{\ttfamily\makestarlowast\makequotesstraight}% default
+\def\verb
+{%
+  \relax \ifmmode\else\leavevmode\null\fi
+  \bgroup
+  \let\do\@makeother \dospecials
+  \@ifstar{\@sverb}% unused
+          {\MicroFont
+           \catcode 32 10 \endlinechar 32 % allows to fetch across line breaks
+           \frenchspacing
+           \@@jfverb}%
+}%
+% Note (Oct 12, 2014): in the improbable situation a newlinechar is
+% found in the ##1, \scantokens will convert this to an end of line in
+% its "write" phase, which will be then ignored in its "read" phase due
+% to \endlinechar-1. This also avoids possible creation of \par which
+% would defeat \@@jfverb@@. Thus it is good.
+\def\@@jfverb #1{%
+   \ifcat\noexpand#1\noexpand~\catcode`#1\active\fi
+% No problem with the EOL for the line where the short verb delimiter stands.
+   \def\next ##1#1{%
+            \@vobeyspaces\everyeof{\relax}\endlinechar\m@ne
+            \expandafter\@@jfverb_a\scantokens\expandafter{##1}}%
+% hack with \@empty to prevent brace stripping if catcodes have been
+% frozen earlier, like in footnotes.
+   \next \@empty
+}
+% We don't want a \discretionary at the very start.
+% But then an empty argument is forbidden!
+\def\@@jfverb_a #1{#1\@@jfverb_b }
+
+\def\@@jfverb_b #1{\ifx\relax #1%
+        \egroup
+      \else
+% \penalty\z@, or rather (Oct 11, 2014) but I then adjust the textwidth
+% precisely:
+      \discretionary{\copy\SoftWrapIcon}{}{}%
+        #1\expandafter\@@jfverb_b\fi
+}
+% \SoftWrapIcon box for line-breaking using discretionaries
+% =========================================================
+\DeclareFontFamily{U}{MdSymbolC}{}
+\DeclareFontShape {U}{MdSymbolC}{m}{n}{<-> MdSymbolC-Regular}{}
+\newbox\SoftWrapIcon
+% Emacs/AUCTeX uses very strange comment-like highlighting for \usefont{U}...
+\def\SetSoftWrapIcon{%
+    \global\setbox\SoftWrapIcon\hb@xt@\z@
+    {\hb@xt@\fontdimen2\font
+        {\hss{\color{verbsoftwrapiconcolor}%
+              \usefont{U}{MdSymbolC}{m}{n}\char"97}\hss}%
+     \hss}%
+   }
+\AtBeginDocument {{\ttfamily\SetSoftWrapIcon}}%
+
+\catcode`_ 8
+\makeatother
+
+% everbatim environment
+% =====================
+
+% October 13-14, 2014
+% Verbatim with an \everypar hook, mainly to have background color, followed by
+% execution of the contents (not limited by a group-scope)
+
+\makeatletter
+\catcode`_ 11
+
+\def\everbatimtop {\MacroFont\small}% default
+\let\everbatimbottom\empty
+\let\everbatimhook\empty
+
+\def\everbatim {\s@everbatim\@everbatim }
+\@namedef{everbatim*}{\s@everbatim\expandafter\@everbatimx\expandafter
+                      {\the\newlinechar}}
+\let\everbatimbgcolorcmd\empty
+\def\everbatimeverypar{\strut
+                   {\everbatimbgcolorcmd\vrule\@width\linewidth }%
+                   \kern-\linewidth
+                   \kern\everbatimindent }
+\def\everbatimindent {\z@}
+% voir plus loin atbegindocument
+
+\def\endeverbatim  {\if@newlist \leavevmode\fi\endtrivlist }
+
+\@namedef{endeverbatim*}{\endeverbatim\aftergroup\everbatimundoparskip}
+\def\everbatimundoparskip{\vbox{}\kern-\baselineskip\kern-\parskip}
+% Note 25 juin 2021
+% On ne peut pas emboîter un everbatim à l'intérieur d'un everbatim
+% ou un everbatim* à l'intérieur d'un everbatim*...
+\def\s@everbatim {%
+%     \ineverbtrue
+     \everbatimtop % put there size changes
+       \topsep    \z@skip
+       \partopsep \z@skip
+       \itemsep   \z@skip
+       \parsep    \z@skip
+       \parskip   \z@skip
+       \lineskip  \z@skip
+     \let\do\@makeother \dospecials
+     \let\do\do@noligs  \verbatim@nolig@list
+     \makestarlowast
+     \makequotesstraight
+     \everbatimhook
+     \trivlist
+       \@topsepadd \z@skip
+       \item\relax
+       \leftskip   \@totalleftmargin
+       \rightskip  \z@skip
+       \parindent  \z@
+       \parfillskip\@flushglue
+       \parskip    \z@skip
+       \@@par
+       \def\par{\leavevmode\null\@@par\pagebreak[1]}%
+       \everypar\expandafter{\the\everypar \unpenalty
+                \everbatimeverypar
+                \everypar \expandafter{\the\everypar\everbatimeverypar}%
+       }%
+       \obeylines \@vobeyspaces
+}
+
+\begingroup
+\lccode`X 13
+\catcode`X \active
+\lccode`Y `* % this is because of \makestarlowast.
+% I have to think whether this is useful: obviously if I were to provide
+% everbatim and everbatim*  in a package I wouldn't do that.
+\catcode`Y  \active
+\catcode`| 0 \catcode`[ 1 \catcode`] 2 \catcode`* 12
+\catcode`{ 12 \catcode`} 12 |catcode`\\ 12
+|lowercase[|endgroup% both freezes catcodes and converts X to active ^^M
+|def|@everbatim #1X#2\end{everbatim}%
+  [#2|end[everbatim]|everbatimbottom ]
+|def|@everbatimx #1#2X#3\end{everbatimY}]%
+  {#3\end{everbatim*}%
+     \everbatimbottom
+     \newlinechar 13
+% execution as LaTeX code of contents
+     \everbatimxprehook
+     \scantokens {#3}%
+     \newlinechar #1\relax
+     \everbatimxposthook
+}%
+
+% L'espace venant du endofline final mis par \scantokens sera inhibé si #3 se
+% termine par un % ou un \x, etc...
+
+\let\everbatimxfgcolorcmd\empty
+\def\everbatimxprehook {\colorlet{everbsavedcolor}{.}%
+                        \everbatimxfgcolorcmd
+                        \smallskip % pour polexpr-examples.tex
+                                   % à cause de problèmes avec les
+                                   % output en "display"
+                       }%
+                       
+\def\everbatimxposthook {\color{everbsavedcolor}}
+
+
+\catcode`_ 8
+\makeatother
+
+\newcommand\pkg[2][]{\if\relax\detokenize{#1}\relax
+                       \href{https://www.ctan.org/pkg/#2}{#2}%
+                     \else
+                       \href{https://www.ctan.org/pkg/#1}{#2}%
+                     \fi
+                    }
+
+% Colors for \verb and everbatim
+% \MacroFont and \MicroFont
+% font size in verbatim blocks
+
+% \verb
+%\colorlet{verbcolor}{DarkCyan}
+\colorlet{verbcolor}{black}
+\colorlet{verbsoftwrapiconcolor}{DarkBlue}
+\def\MicroFont{\ttfamily%\color{verbcolor}
+               \makestarlowast\makequotesstraight}%
+
+% everbatim/everbatim*
+\def\everbatimtop{\MacroFont\small}
+% the \small is not in \MacroFont in case of a document with macrocode (doc.sty)
+% and some customization is desired
+%\colorlet{everbatimfgcolor}{Olive}
+\colorlet{everbatimfgcolor}{DarkBlue}
+\def\MacroFont{\ttfamily\color{everbatimfgcolor}}
+
+%\colorlet{everbatimbgcolor}{WhiteSmoke}
+\colorlet{everbatimbgcolor}{Ivory}
+\def\everbatimbgcolorcmd{\color{everbatimbgcolor}}
+
+%\colorlet{everbatimxfgcolor}{MidnightBlue}
+\colorlet{everbatimxfgcolor}{OliveDrab}
+\def\everbatimxfgcolorcmd{\color{everbatimxfgcolor}}
+
+% Notice that \macrocode uses \macro@font which stores the \MacroFont meaning
+% in force at \begin{document}. But doc.sty's verbatim uses current \MacroFont
+% not the meaning at \begin{document}. Comprenne qui pourra...
+
+\begin{document}
+\maketitle
+
+
+The package provides a parser |\poldef| of algebraic polynomial
+expressions.  
+
+Once defined, a polynomial is usable by its name either as a numerical
+function in |\xintexpr/\xinteval|, or for additional polynomial
+definitions, or as argument to the package macros.
+
+%   The localization of
+% real roots to arbitrary precision as well as the determination of all
+% rational roots is implemented via such macros.
+
+% Since release |0.8|, polexpr extends the \pkg{xintexpr}
+% syntax to recognize
+% polynomials as a new variable type (and not only as functions).
+% Functionality which previously was implemented via macros such as the
+% computation of a greatest common divisor is now available directly in
+% |\xintexpr|, |\xinteval| or |\poldef| via infix or functional
+% syntax.
+
+This document illustrates root localization via usage of macros such as
+|\PolToSturm| and |\PolSturmIsolateZeros| which implement the
+\href{https://en.wikipedia.org/wiki/Sturm%27s_theorem}{Sturm theorem}:
+\begin{itemize}
+\item Root localization based on
+  \href{https://en.wikipedia.org/wiki/Sturm%27s_theorem}{Sturm theorem} was
+  added at release |0.4| (2018/02/16).
+\item Ability to find all rational roots was added at release |0.7.2|
+  (2018/12/09).
+\end{itemize}
+As of |0.8| (2021/03/29), \pkg{polexpr} is usable with Plain \TeX\ and not
+only with \LaTeX. The examples here use most of the time a syntax which works
+with both.
+
+Copying-pasting from |pdf| the example source may lose formatting.
+Formerly, they were included verbatim in the |html| documentation.  Here
+they are both rendered verbatim and got executed during the \LaTeX\ run
+which created this |pdf| file, with the output shown after the source code.
+
+% Perhaps future releases will implement other approaches, which are known
+% to be generically computationally more efficient, at least in high
+% degrees, than the \href{https://en.wikipedia.org/wiki/Sturm%27s_theorem}{Sturm theorem} based approach.  This is not
+% immediate priority though (perhaps support of multivariate polynomials
+% would be more important feature; or localization of complex roots).
+
+Regarding how polynomial coefficients are printed on the typeset page by
+|\PolTypeset|:
+\begin{itemize}\def\everbatimtop {\MacroFont}% sans le \small
+\item The default for |\PolTypesetOne| is to use |\xintTeXsignedFrac| with
+  \LaTeX, |\xintTeXsignedOver| with Plain.  See the \pkg{xintexpr}
+  documentation for a description of what these macros do.  A sensible
+  definition is: 
+\begin{everbatim}
+  \def\PolTypesetOne#1{\PolDecToString{\xintREZ{#1}}}%
+\end{everbatim}
+% le \smallskip est ennuyeux ici
+  It means to use decimal notation, with perhaps a trailing denominator if the
+  argument is a fraction, and will suppress trailing zeros after the decimal
+  mark.
+
+\item
+  As these are expandable macros, they are usable to redefine
+  |\PolToExprCmd| as well:
+\begin{everbatim}
+  \def\PolToExprCmd#1{\PolDecToString{\xintREZ{#1}}}%
+\end{everbatim}
+% le \smallskip est ennuyeux ici
+  This will customize the output of |\PolToExpr| (which a priori is
+  destined for writes to external files but may also be used on the typeset
+  page).
+\end{itemize}
+
+With |\xintverbosetrue| in the \TeX\ source extra information relative
+to the internal data manipulated by the macros will be written to the |.log|
+file.
+
+\begin{framed}
+  Package macros related to root localization create (user-level) new
+  polynomials, or numeric variables, via a naming scheme using the given
+  |<sturmname>| as prefix. It is thus advisable to keep this
+  |<sturmname>| name-space separate from the one used to name polynomial or
+  scalar variables.
+\end{framed}
+
+\begin{framed}
+  Regrettably all examples here use the condemnable |\PolToSturm{f}{f}|
+  practice which means that internally defined polynomials will use as prefix
+  the original polynomial name.  This merge of namespaces may cause
+  overwriting previously defined data and may lead to hard-to-debug problems.
+\end{framed}
+
+
+\section{A first example}
+
+
+In this example the polynomial is square-free.
+\begin{everbatim*}
+  \poldef f(x) := x^7 - x^6 - 2x + 1;
+
+  \PolToSturm{f}{f}
+  \PolSturmIsolateZeros{f}
+  The \PolTypeset{f} polynomial has \PolSturmNbOfIsolatedZeros{f} distinct real
+  roots which are located in the following intervals:
+  \PolPrintIntervals{f}
+\end{everbatim*}
+
+\begin{everbatim*}
+  Here is the second root with ten more decimal digits:
+  \PolRefineInterval[10]{f}{2}
+  $$\PolSturmIsolatedZeroLeft{f}{2}<Z_2<\PolSturmIsolatedZeroRight{f}{2}$$
+\end{everbatim*}
+
+\begin{everbatim*}
+  And here is the first root with twenty digits after decimal mark:
+  \PolEnsureIntervalLength{f}{1}{-20}
+  $$\PolSturmIsolatedZeroLeft{f}{1}<Z_1<\PolSturmIsolatedZeroRight{f}{1}$$
+\end{everbatim*}
+
+\begin{everbatim*}
+  The first element of the Sturm chain has degree $\PolDegree{f_0}$. As
+  this is the original degreee $\PolDegree{f}$ we know that $f$ is square free.
+  Its derivative is up to a constant \PolTypeset{f_1} (in this example
+  it is identical with it).
+  \PolToSturm{f_1}{f_1}\PolSturmIsolateZeros{f_1}%
+  The derivative has \PolSturmNbOfIsolatedZeros{f_1} distinct real
+  roots:
+  \PolPrintIntervals[W]{f_1}%
+  \PolEnsureIntervalLengths{f_1}{-10}%
+  Here they are with ten digits after decimal mark:
+  \PolPrintIntervals[W]{f_1}
+\end{everbatim*}
+
+\begin{everbatim*}
+  \PolDiff{f_1}{f''}
+  \PolToSturm{f''}{f''}
+  \PolSturmIsolateZeros{f''}
+  The second derivative is \PolTypeset{f''}.
+  It has \PolSturmNbOfIsolatedZeros{f''} distinct real
+  roots:
+  \PolPrintIntervals[X]{f''}%
+  Here is the positive one with 20 digits after decimal mark:
+  \PolEnsureIntervalLength{f''}{2}{-20}%
+  $$X_2 = \PolSturmIsolatedZeroLeft{f''}{2}\dots$$
+\end{everbatim*}
+%The more mathematically advanced among our dear readers will be able
+%  to give the exact value for $X_2$!
+
+\section{A degree four polynomial with nearby roots}
+
+
+Notice that this example is a bit outdated as |0.7| release has
+added |\PolSturmIsolateZeros**{<sturmname>}| which would find exactly
+the roots. The steps here retain their interest when one is interested
+in finding isolating intervals for example to prepare some demonstration
+of dichotomy method.
+
+
+\begin{everbatim*}
+  \PolDef{Q}{(x-1.050001)(x-1.105001)(x-1.110501)(x-1.111051)}
+  \PolTypeset{Q}
+  \PolToSturm{Q}{Q} % it is allowed to use same prefix for Sturm chain
+  \PolSturmIsolateZeros{Q}
+  \PolPrintIntervals{Q}
+\end{everbatim*}
+
+\begin{everbatim*}
+  \PolRefineInterval*{Q}{1}
+  \PolRefineInterval*{Q}{2}
+  \PolRefineInterval*{Q}{3}
+  \PolRefineInterval*{Q}{4}
+  \PolPrintIntervals{Q}
+\end{everbatim*}
+
+\begin{everbatim*}
+  \PolEnsureIntervalLengths{Q}{-6}
+  \PolPrintIntervals{Q}
+  % finds here all roots exactly
+\end{everbatim*}
+
+\section{The degree nine polynomial with 0.99, 0.999, 0.9999 as triple roots}
+
+Define a user command (\pkg[xint]{xinttools} is loaded automatically by
+\pkg{polexpr}):
+
+\begin{everbatim*}
+  \def\showmultiplicities#1{% #1 = "sturmname"
+  \xintFor* ##1 in {\xintSeq{1}{\PolSturmNbOfIsolatedZeros{#1}}}\do{%
+      The multiplicity is \PolSturmIsolatedZeroMultiplicity{#1}{##1}
+      \PolSturmIfZeroExactlyKnown{#1}{##1}%
+      {at the root $x=\PolSturmIsolatedZeroLeft{#1}{##1}$}
+      {for the root such that
+      $\PolSturmIsolatedZeroLeft{#1}{##1}<x<\PolSturmIsolatedZeroRight{#1}{##1}$}
+      \par
+  }}%
+\end{everbatim*}
+
+\begin{everbatim*}
+  \PolDef{f}{(x-0.99)^3(x-0.999)^3(x-0.9999)^3}
+  \def\PolTypesetOne#1{\PolDecToString{\xintREZ{#1}}}
+  \PolTypeset{f}\par
+\end{everbatim*}
+
+\begin{everbatim*}
+  \PolToSturm{f}{f}% it is allowed to use "polname" as "sturmname" too
+  \PolSturmIsolateZerosAndGetMultiplicities{f}% use the "sturmname" here
+  % or \PolSturmIsolateZeros*{f} which is exactly the same, but shorter..
+
+  \showmultiplicities{f}
+\end{everbatim*}
+% In this example, the output will look like this (but using math mode)::
+
+%   x^9 - 8.9667x^8 + 35.73400293x^7 - 83.070418400109x^6 + 124.143648875193123x^5
+%   - 123.683070924326075877x^4 + 82.149260397553075617891x^3 
+%   - 35.07602992699900159127007x^2 + 8.7364078733314648368671733x
+%   - 0.967100824643585986488103299
+
+%   The multiplicity is 3 at the root x = 0.99
+%   The multiplicity is 3 at the root x = 0.999
+%   The multiplicity is 3 at the root x = 0.9999
+
+% On first pass, these rational roots were found (due to their relative
+% magnitudes, using |\PolSturmIsolateZeros**| was not needed here). But
+% multiplicity computation works also with (decimal) roots not yet
+% identified or with non-decimal or irrational roots.
+
+It is fun to modify only a tiny bit the polynomial and see if polexpr
+survives:
+
+\begin{everbatim*}
+  \PolDef{g}{f(x)+1e-27}
+  \PolTypeset{g}\par
+  \PolToSturm{g}{g}
+  \PolSturmIsolateZeros*{g}
+
+  \showmultiplicities{g}
+\end{everbatim*}
+% This produces::
+
+%   x^9 - 8.9667x^8 + 35.73400293x^7 - 83.070418400109x^6 + 124.143648875193123x^5
+%   - 123.683070924326075877x^4 + 82.149260397553075617891x^3 
+%   - 35.07602992699900159127007x^2 + 8.7364078733314648368671733x
+%   - 0.967100824643585986488103298
+
+%   The multiplicity is 1 for the root such that 0.98 < x < 0.99
+%   The multiplicity is 1 for the root such that 0.9991 < x < 0.9992
+%   The multiplicity is 1 for the root such that 0.9997 < x < 0.9998
+
+This means that the multiplicity-3 roots each became a real and a pair of
+complex ones. Let's see them better:
+
+\begin{everbatim*}
+  \PolEnsureIntervalLengths{g}{-10}
+
+  \showmultiplicities{g}
+\end{everbatim*}
+% which produces::
+
+%   The multiplicity is 1 for the root such that 0.9899888032 < x < 0.9899888033
+%   The multiplicity is 1 for the root such that 0.9991447980 < x < 0.9991447981
+%   The multiplicity is 1 for the root such that 0.9997663986 < x < 0.9997663987
+
+\section{A degree five polynomial with three rational roots}
+
+\begin{everbatim*}
+  \poldef Q(x) :=  1581755751184441 x^5
+                 -14907697165025339 x^4
+                 +48415668972339336 x^3
+                 -63952057791306264 x^2
+                 +46833913221154895 x
+                 -49044360626280925;
+
+  \PolToSturm{Q}{Q}
+    \def\PolTypesetCmdPrefix#1{\allowbreak\xintiiifSgn{#1}{}{+}{+}}%
+    $Q_0(x) = \PolTypeset{Q_0}$
+  \PolSturmIsolateZeros**{Q}
+  \PolPrintIntervals{Q}
+
+  $Q_{norr}(x) = \PolTypeset{Q_norr}$
+\end{everbatim*}
+
+Here, all real roots are rational. Let's get their decimal expansion too:
+\begin{everbatim*}
+\begingroup
+  % print decimal expansion of the found roots
+  \def\PolPrintIntervalsPrintExactZero
+              {\xintTrunc{20}{\PolPrintIntervalsTheLeftEndPoint}\dots}
+  \PolPrintIntervals{Q}
+\endgroup % we localized the modified \PolPrintIntervalsPrintExactZero
+\end{everbatim*}
+
+
+  % Z_1 = 3.14159265358107777120...
+  % Z_2 = 3.14159265358979340254...
+  % Z_3 = 3.14159292035398230088...
+
+\section{A Mignotte type polynomial}
+
+
+\begin{everbatim*}
+  \PolDef{P}{x^10 - (10x-1)^2}%
+  \PolTypeset{P}              % prints it in expanded form
+  \PolToSturm{P}{P}           % we can use same prefix for Sturm chain
+  \PolSturmIsolateZeros{P}    % finds 4 real roots
+  This polynomial has \PolSturmNbOfIsolatedZeros{P} distinct real roots:
+  \PolPrintIntervals{P}%
+\end{everbatim*}
+Let us refine the second and third intervals to separate the corresponding
+roots:
+\begin{everbatim*}
+  \PolRefineInterval*{P}{2}% will refine to 0.0999990 < Z_2 < 0.0999991
+  \PolRefineInterval*{P}{3}% will refine to 0.100001 < Z_3 < 0.100002
+  \PolPrintIntervals{P}%
+\end{everbatim*}
+Let us now get to know all roots with 10 digits after decimal mark:
+\begin{everbatim*}
+  \PolEnsureIntervalLengths{P}{-10}%
+  \PolPrintIntervals{P}% now all roots are known 10 decimal digits after mark
+\end{everbatim*}
+Finally, we display 20 digits of the second root:
+\begin{everbatim*}
+  \PolEnsureIntervalLength{P}{2}{-20}% makes Z_2 known with 20 digits after mark
+  $$\PolSturmIsolatedZeroLeft{P}{2}<Z_2<\PolSturmIsolatedZeroRight{P}{2}$$
+\end{everbatim*}
+% The last line produces::
+
+%   0.09999900004999650028 < Z_2 < 0.09999900004999650029
+
+\section{The Wilkinson polynomial}
+
+
+See \url{https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial}.
+
+
+\begin{everbatim*}
+  %\xintverbosetrue % for the curious...
+
+  \poldef f(x) := mul((x - i), i = 1..20);
+
+  \def\PolTypesetCmdPrefix#1{\allowbreak\xintiiifSgn{#1}{}{+}{+}}%
+  \def\PolTypesetOne#1{\xintDecToString{#1}}%
+
+  \noindent\PolTypeset{f}
+
+  \PolToSturm{f}{f}
+  \PolSturmIsolateZeros{f}
+  \PolPrintIntervals{f}
+  % \vfill\eject
+
+  % This page is commented out because it takes about 30s on a 2GHz CPU
+  % \poldef g(x) := f(x) - 2**{-23} x**19;
+
+  % \PolToSturm{g}{g}
+  % \noindent\PolTypeset{g_0}% integer coefficient primitive polynomial
+
+  % \PolSturmIsolateZeros{g}
+  % \PolEnsureIntervalLengths{g}{-10}
+
+  % \PolPrintIntervals*{g}
+\end{everbatim*}
+
+The first polynomial
+%
+%   f(x) = x**20
+%   - 210 x**19
+%   + 20615 x**18
+%   - 1256850 x**17
+%   + 53327946 x**16
+%   - 1672280820 x**15
+%   + 40171771630 x**14
+%   - 756111184500 x**13
+%   + 11310276995381 x**12
+%   - 135585182899530 x**11
+%   + 1307535010540395 x**10
+%   - 10142299865511450 x**9
+%   + 63030812099294896 x**8
+%   - 311333643161390640 x**7
+%   + 1206647803780373360 x**6
+%   - 3599979517947607200 x**5
+%   + 8037811822645051776 x**4
+%   - 12870931245150988800 x**3
+%   + 13803759753640704000 x**2
+%   - 8752948036761600000 x
+%   + 2432902008176640000
+%
+is handled fast enough, but the modified one |f(x) -
+2**-23 x**19| takes about 20x longer.
+
+Its Sturm chain polynomials
+have integer coefficients with up to 321 digits, whereas (surprisingly
+perhaps) those of the Sturm chain polynomials derived from |f| never
+have more than 21 digits ...
+
+Once the Sturm chain is computed and the zeros isolated, obtaining their
+decimal digits is relatively faster. Here are the ten real roots of
+|f(x) - 2**-23 x**19| which would be computed by the commented-out code above:
+\begin{everbatim}
+  Z_1 = 0.9999999999...
+  Z_2 = 2.0000000000...
+  Z_3 = 2.9999999999...
+  Z_4 = 4.0000000002...
+  Z_5 = 4.9999999275...
+  Z_6 = 6.0000069439...
+  Z_7 = 6.9996972339...
+  Z_8 = 8.0072676034...
+  Z_9 = 8.9172502485...
+  Z_10 = 20.8469081014...
+\end{everbatim}
+
+\section{The second Wilkinson polynomial}
+
+
+\begin{everbatim*}
+  \poldef f(x) := mul(x - 2^-i, i = 1..20);
+
+  %\PolTypeset{f}
+
+  \PolToSturm{f}{f}
+  \PolSturmIsolateZeros**{f}
+  \PolPrintIntervals{f}
+\end{everbatim*}
+
+This takes more time than the polynomial with 1, 2, .., 20 as roots but
+less than the latter modified by the |2**-23| tiny change to one of its
+coefficient.
+
+% Here is the output (with release 0.7.2)::
+
+%   Z_1  = 0.00000095367431640625
+%   Z_2  = 0.0000019073486328125
+%   Z_3  = 0.000003814697265625
+%   Z_4  = 0.00000762939453125
+%   Z_5  = 0.0000152587890625
+%   Z_6  = 0.000030517578125
+%   Z_7  = 0.00006103515625
+%   Z_8  = 0.0001220703125
+%   Z_9  = 1/4096
+%   Z_10 = 1/2048
+%   Z_11 = 1/1024
+%   Z_12 = 1/512
+%   Z_13 = 1/256
+%   Z_14 = 1/128
+%   Z_15 = 0.015625
+%   Z_16 = 0.03125
+%   Z_17 = 0.0625
+%   Z_18 = 0.125
+%   Z_19 = 0.25
+%   Z_20 = 0.5
+
+There is some incoherence in output format which has its source in the
+fact that some roots are found in branches which can only find decimal
+roots, whereas some are found in branches which could find general
+fractions and they use |\xintIrr| before storage of the found root.
+This may evolve in future.
+
+
+\section{The degree 41 polynomial with -2, -1.9, -1.8, ..., 0, 0.1, ..., 1.9, 2 as roots}
+
+
+
+\begin{everbatim*}
+  \PolDef{P}{mul((x-i*1e-1), i=-20..20)}% i/10 is same but less efficient
+\end{everbatim*}
+In the defining expression we could have used |i/10| but this gives
+less efficient internal form for the coefficients (the |10|'s end up
+in denominators).
+
+\begin{everbatim*}
+\begingroup
+  \def\PolToExprCmd#1{\PolDecToString{\xintREZ{#1}}}
+  \def\PolToExprTermPrefix#1{\newline\xintiiifSgn{#1}{}{+}{+}}
+  \def\PolToExprTimes{${}\cdot{}$}
+\ttfamily
+\PolToExpr{P}
+\endgroup
+\end{everbatim*}
+  % x^41
+  % -28.7*x^39
+  % +375.7117*x^37
+  % -2975.11006*x^35
+  % +15935.28150578*x^33
+  % -61167.527674162*x^31
+  % +173944.259366417394*x^29
+  % -373686.963560544648*x^27
+  % +613012.0665016658846445*x^25
+  % -771182.31133138163125495*x^23
+  % +743263.86672885754888959569*x^21
+  % -545609.076599482896371978698*x^19
+  % +301748.325708943677229642930528*x^17
+  % -123655.8987669450434698869844544*x^15
+  % +36666.1782054884005855608205864192*x^13
+  % -7607.85821367459445649518380016128*x^11
+  % +1053.15135918687298508885950223794176*x^9
+  % -90.6380005918141132650786081964032*x^7
+  % +4.33701563847327366842552218288128*x^5
+  % -0.0944770968420804735498178265088*x^3
+  % +0.00059190121813899276854174416896*x
+
+which shows coefficients with up to 36 significant digits...
+
+Stress test: not a hard challenge to \pkg{xint} + \pkg{polexpr}, but be a bit
+patient!
+
+\begin{everbatim*}
+  \PolDef{P}{mul((x-i*1e-1), i=-20..20)}%
+  \PolToSturm{P}{S}           % dutifully computes S_0, ..., S_{41}
+  % the [1] optional argument limits the search to interval (-10,10)
+  \PolSturmIsolateZeros[1]{S} % finds *exactly* (but a bit slowly) all 41 roots!
+  \PolPrintIntervals{S}       % nice, isn't it?
+% Unfortunately \PolPrintIntervals uses a non-breakable array environment
+% But see next section on how to customize \PolPrintIntervals and let it
+% allow pagebreaks
+\end{everbatim*}
+
+\begin{quote}
+  Release |0.5| has \emph{experimental} addition of optional argument |E| to
+  |\PolSturmIsolateZeros|. It instructs to search roots only in interval
+  |(-10^E, 10^E)|. Important: the extremities are \emph{assumed to not be roots}.
+  In this example, the |[1]| in |\PolSturmIsolateZeros[1]{S}| gives some speed
+  gain; without it, it turns out in this case that \pkg{polexpr} would have
+  started with |(-10^6, 10^6)| interval.
+
+  Please note that this feature may be removed or modified.
+\end{quote}
+
+\section{Roots of a Chebyshev polynomial}
+
+
+\begin{everbatim*}
+  \poldef T_0(x) := 1;
+  \poldef T_1(x) := x;
+  \catcode`@ 11
+  \count@ 2
+  \xintloop
+    \poldef T_\the\count@(x) :=
+            2x*T_\the\numexpr\count@-1\relax
+             - T_\the\numexpr\count@-2\relax;
+  \ifnum\count@<15
+  \advance\count@ 1
+  \repeat
+  \catcode`@ 12
+
+  $$T_{15} = \PolTypeset[X]{T_15}$$
+  \PolToSturm{T_15}{T_15}
+  \PolSturmIsolateZeros*{T_15}% "*" as we will want to confirm multiplicity one
+  % takes time (each next decimal digit is obtained by dichotomy)
+  \PolEnsureIntervalLengths{T_15}{-20}% ensure 20 decimal digits for each root
+\end{everbatim*}
+
+Here is now an example of customization.  Indeed
+|\PolPrintIntervals| default uses |array| and thus does not allow page breaks.
+And it uses |Left < Z < Right| as presentation of roots and we would like here
+rather |Z = decimal expansion...|.
+\begin{everbatim*}
+% 0.8.6 adds an internal patch which would allow usage of amsmath environments
+% like this:
+% \def\PolPrintIntervalsBeginEnv{\begin{align*}}
+% \def\PolPrintIntervalsEndEnv{\end{align*}}
+% (the problem was that align evaluates twice its contents so global variables
+%  need a reset at the end of first pass, which is what 0.8.6 took care of)
+%
+% Let's simply do this:
+\def\PolPrintIntervalsBeginEnv{\begingroup\leftskip3cm\relax}
+\def\PolPrintIntervalsEndEnv{\par\endgroup}
+%
+% The rows are separated by \PolPrintIntervalsRowSeparator which defaults to \\
+% with LaTeX and \cr with Plain. (prior to 0.8.6 it was hardcoded)
+\def\PolPrintIntervalsRowSeparator{\\[\jot]}
+% And we enter math mode manually at each row, copying pasting from package
+% defaults with some added mathon/mathoff:
+\def\PolPrintIntervalsKnownRoot{%
+   $\PolPrintIntervalsPrintMultiplicity\quad
+   \PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}=%
+   \PolPrintIntervalsPrintExactZero
+   $
+}%
+\def\PolPrintIntervalsUnknownRoot{%
+   $\PolPrintIntervalsPrintMultiplicity\quad
+   \xintifSgn{\PolPrintIntervalsTheLeftEndPoint}%
+     {\xintifSgn{\PolPrintIntervalsTheRightEndPoint}
+       {\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}=%
+        \PolPrintIntervalsPrintRightEndPoint\dots}%
+      {0>\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}>%
+        \PolPrintIntervalsPrintLeftEndPoint}%
+      {\impossibleA}}%
+    {\xintifSgn{\PolPrintIntervalsTheRightEndPoint}
+      {\impossibleB}%
+      {\impossibleC}%
+      {0<\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}<%
+        \PolPrintIntervalsPrintRightEndPoint}}%
+    {\xintifSgn{\PolPrintIntervalsTheRightEndPoint}
+      {\impossibleD}%
+      {\impossibleE}%
+      {\PolPrintIntervalsTheVar_{\PolPrintIntervalsTheIndex}=%
+        \PolPrintIntervalsPrintLeftEndPoint\dots}}%
+   $
+}%
+\PolPrintIntervals{T_15}
+\end{everbatim*}
+\end{document}
+
diff --git a/Master/texmf-dist/tex/generic/polexpr/polexpr.sty b/Master/texmf-dist/tex/generic/polexpr/polexpr.sty
index 04d189de8bb..7d3dbb36bda 100644
--- a/Master/texmf-dist/tex/generic/polexpr/polexpr.sty
+++ b/Master/texmf-dist/tex/generic/polexpr/polexpr.sty
@@ -2,7 +2,7 @@
 % License: LPPL 1.3c (author-maintained)
 % Usage: \input polexpr.sty   (Plain or other macro formats)
 %    or  \usepackage{polexpr} (LaTeX macro format)
-% Release 0.8.5 (2021/11/30) of polexpr.sty. This file inputs
+% Release 0.8.6 (2022/01/09) of polexpr.sty. This file inputs
 %   polexprcore.tex
 %   polexprexpr.tex
 %   polexprsturm.tex
@@ -54,12 +54,12 @@
 \XINTsetupcatcodes% (does \endlinechar13 in particular)
 \XINT_providespackage
 \ProvidesPackage{polexpr}%
-  [2021/05/27 v0.8.3 Polynomial expressions with rational coefficients (JFB)]%
+  [2022/01/09 v0.8.6 Polynomial expressions with rational coefficients (JFB)]%
 \begingroup
  \def\x#1/#2/#3 #4\xint:{#1#2#3}%
  \ifnum\expandafter\x\expanded{\csname ver@xintexpr.sty\endcsname}\xint:
        <20210527 % xint 1.4h
- \immediate\write128{! Package polexpr error: xintexpr too old, aborting input}%
+ \errmessage{Package polexpr error: xintexpr too old, aborting input}%
  \else\expandafter\xint_gobble_i
  \fi
 \endinput\endgroup
@@ -475,7 +475,7 @@
 }%
 %
 %% Euclidean division
-%  now based on the expandable routine from polexprcore.tex
+%  since 0.8 based on the expandable routine from polexprcore.tex
 %
 \def\PolDivide#1#2#3#4{% #3=quotient, #4=remainder of #1 by #2
     \POL@divide{#1}{#2}%
diff --git a/Master/texmf-dist/tex/generic/polexpr/polexprcore.tex b/Master/texmf-dist/tex/generic/polexpr/polexprcore.tex
index 3648bc8ccb2..8dcc893844b 100644
--- a/Master/texmf-dist/tex/generic/polexpr/polexprcore.tex
+++ b/Master/texmf-dist/tex/generic/polexpr/polexprcore.tex
@@ -1,7 +1,9 @@
-%% This file polexprcore.tex is part of the polexpr package (0.8.5, 2021/11/30)
-%% Core routines to match infix operators +, -, *, //, /:, ^, ** and some
-%% functions
-%% The atoms representing polynomials inside \xintexpr are
+%% filename: polexprcore.tex
+%% Part of the polexpr package (0.8.6, 2022/01/09)
+%%
+%% Core routines for infix operators +, -, *, //, /:, ^, ** and functions
+%%
+%% Memo: the atoms representing polynomials inside \xintexpr are
 %% - for constants: a numeric value  (indistinguishable. from scalars)
 %% - for degree at least 1: P<degree>.{c0}{c1}....{cN} with N = degree
 %% Auxiliaries
diff --git a/Master/texmf-dist/tex/generic/polexpr/polexprexpr.tex b/Master/texmf-dist/tex/generic/polexpr/polexprexpr.tex
index f3cc5ac28ef..bbc860cbca4 100644
--- a/Master/texmf-dist/tex/generic/polexpr/polexprexpr.tex
+++ b/Master/texmf-dist/tex/generic/polexpr/polexprexpr.tex
@@ -1,5 +1,7 @@
-%% This file polexprexpr.tex is part of the polexpr package (0.8.5, 2021/11/30)
-%% Extending \xintexpr syntax:
+%% filename: polexprexpr.tex
+%% Part of the polexpr package (0.8.6, 2022/01/09)
+%%
+%% Polynomial extensions to the \xintexpr syntax:
 %%
 %% 1. Authorize ' in variable and function names
 %%    This partially breaks infix operators 'and', 'or', 'xor', 'mod'
diff --git a/Master/texmf-dist/tex/generic/polexpr/polexprsturm.tex b/Master/texmf-dist/tex/generic/polexpr/polexprsturm.tex
index d590d48aa56..0edc8b8f0a2 100644
--- a/Master/texmf-dist/tex/generic/polexpr/polexprsturm.tex
+++ b/Master/texmf-dist/tex/generic/polexpr/polexprsturm.tex
@@ -1,7 +1,11 @@
-%% This file polexprsturm.tex is part of the polexpr package (0.8.5, 2021/11/30)
-%% Sturm Algorithm (polexpr 0.4)
+%% filename: polexprsturm.tex
+%% Part of the polexpr package (0.8.6, 2022/01/09)
+%%
+%% Implements the Sturm localization Algorithm
+%% Added at polexpr 0.4
+%%
 %% 0.5 uses primitive polynomials for faster evaluations afterwards
-%% 0.6 corrects misuse of \@ifstar! (mumble). \PolToSturm* was broken.
+%% 0.6 corrects misuse of \@ifstar (mumble). \PolToSturm* was broken.
 %% 0.6's \PolToSturm* defines both normalized and unnormalized, the
 %%       unnormalized using two underscores, so both are available
 %%       Sole difference is that \PolToSturm* also declares them as 
@@ -9,19 +13,20 @@
 %%       holding the coefficients in memory
 %% 0.6 fixes the case of a constant polynomial P which caused division
 %%     by zero error from P'.
-%% 0.8 - fixes 0.7.5 failure to have updated to xint 1.4 format the defined
-%%     \xintexpr variables holding the localization intervals extremities
-%%     - also, it uses the prem() in computing the Sturm chain, for a 3X
-%%     speed gain in the case of the "perturbed" first Wilkinson example
-%%
+%% 0.8 - fixes 0.7.5 compatibility bug with xint 1.4 internal format
+%%       regarding the defined \xintexpr variables holding the localization
+%%       intervals extremities
+%%     - also, it uses the prem() in computing the Sturm chain, with a 3X
+%%       speed gain in the case of the "perturbed" first Wilkinson example
+%% 0.8.6 has better a priori bounds for positive and negative roots
 \newcount\POL@count
 \newif\ifPOL@tosturm@makefirstprimitive\POL@tosturm@makefirstprimitivetrue
 \newif\ifPOL@isolz@nextwillneedrefine
 %%
 \def\PolToSturm{\POL@ifstar{\PolToSturm@@}{\PolToSturm@}}%
 \def\POL@aux@toint#1{\xintREZ{\xintNum{#1}}}% for polynomials with int. coeffs!
-%% Attention that some macros rely upon this one setting \POL@sturmname
-%% and \POL@sturm@N as it does
+%% Attention that some macros rely upon this one defining \POL@sturmname
+%% and \POL@sturm@N as it currently does
 \def\PolToSturm@#1#2{%
   \edef\POL@sturmname{#2}%
   % 0.6 uses 2 underscores (one before index, one after) to keep in memory
@@ -61,11 +66,9 @@
 }%
 \def\POL@tosturm@dosturm{%
   \POL@Diff@@one{\POL@sturmname _0_}{\POL@sturmname _1_}%
-  % re-utiliser \POL@varcoeffs directement?
   \POL@makeprimitive{\POL@sturmname _1_}% does not do \POL@newpol
   \POL@count\@ne
   \xintloop
-  % prior to 0.8, code was using here \POL@divide
     \POL@getprem{\POL@sturmname _\the\numexpr\POL@count-\@ne\relax _}%
                 {\POL@sturmname _\the\POL@count _}%
     \expandafter\POL@split\POL@R;\POL@degR\POL@polR
@@ -117,7 +120,8 @@
   \POL@sturmchain@getSV@at\POL@sturmchain@X
   #1\let#2\POL@sturmchain@SV
 }%
-\def\POL@sturmchain@getSV@at#1{% ATTENTION USES \POL@count
+% attention that this modifies current \POL@count value
+\def\POL@sturmchain@getSV@at#1{%
   \def\POL@sturmchain@SV{0}%
   \edef\POL@sturmchain@sign{\xintiiSgn{\POL@eval{\POL@sturmname _0}{#1}}}%
   \let\POL@isolz@lastsign\POL@sturmchain@sign
@@ -166,7 +170,6 @@
       {\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}}%
@@ -196,10 +199,6 @@
   \edef\POL@sturm@N{\@nameuse{PolSturmChainLength_\POL@sturmname}}%
   % isolate the roots (detects case of constant polynomial)
   \PolSturmIsolateZeros@{\POL@sturmname}%
-  % 0.8.4 fix: these declarations were formerly not executed in absence of roots!
-     % 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)
      \XINT_global
      \expandafter\let
          \csname POLuserpol@\POL@sturmname\POL@sqfnorr\expandafter\endcsname
@@ -219,7 +218,6 @@
      \begingroup\globaldefs\@ne
      \expandafter\POL@initarray\csname POL_ZM\POL@sturmname*\endcsname{1}%
      \endgroup
-     % 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
@@ -227,12 +225,10 @@
         {\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}}%
-% ATTENTION QUE LA CONVENTION DE SIGNE POUR \POL@findrat@E EST OPPOSÉE À CELLE
-% POUR LE CODE PLUS ANCIEN FAISANT "REFINE"
      \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
+     % factor them out from original (Sturm root) polynomial
      \ifnum\POL@findrat@E<7 %
          \def\POL@findrat@index{1}%
          \POL@findrat@loop@secondpass@direct
@@ -252,7 +248,6 @@
      \POL@newpol{\POL@sturmname\POL@norr}% with multiplicities
 }%
 \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_ZK\POL@sturmname*\xintiloopindex\endcsname
@@ -261,8 +256,6 @@
     \ifnum\xintiloopindex<\POL@isolz@NbOfRoots\space
     \repeat }%
   \begingroup\globaldefs\@ne
-  % 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
 }%
@@ -286,11 +279,6 @@
 }%
 \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@_
@@ -313,15 +301,13 @@
           \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@getE #1/1[#2]{#2}%
 \def\POL@findrat@loop@a{%
   % attention that the width may have been already smaller than 10^{-6}
   \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
@@ -364,9 +350,7 @@
   \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^9 au pire donc ok pour \numexpr
-  % we will drop integral part in our updating P
+% 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
@@ -382,7 +366,6 @@
   \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}%
@@ -474,7 +457,8 @@
        {\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
+% using ii version is an abuse
+  \xintiiifNeg{\POL@findrat@Lscaled}%
     {% 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]}%
@@ -572,13 +556,8 @@
   % 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_
@@ -640,7 +619,7 @@
   \fi
 }%
 \def\POL@isolzmult@defvar@M{%
-  % Attention that is used not only in ...GetMultiplicities@ but also
+  % Attention that this is used not only in ...GetMultiplicities@ but also
   % in FindRationalRoots
   \begingroup\xintglobaldefstrue
   % added at 0.7
@@ -739,7 +718,7 @@
   % #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 from a given polynomial)
-  % For reasons I have forgotten (no time now) this code **must** be used
+  % For reasons I have forgotten this code **must** be used
   % with a *normalized* Sturm chain.
   \edef\POL@sturmname{#2}%
   \edef\POL@sturmlength{\PolSturmChainLength{#2}}%
@@ -749,22 +728,23 @@
     \POL@isolz@getsignchanges@plusinf
     \POL@isolz@getsignchanges@minusinf
    \else
-    \edef\POL@isolz@E{\the\numexpr\xint_zapspaces #1 \xint_gobble_i\relax}%
-    \POL@sturmchain@getSV@at{1[\POL@isolz@E]}%
+    \edef\POL@isolz@E@pos{\the\numexpr\xint_zapspaces #1 \xint_gobble_i\relax}%
+    \let\POL@isolz@E@neg\POL@isolz@E@pos
+    \POL@sturmchain@getSV@at{1[\POL@isolz@E@pos]}%
     \let\POL@isolz@plusinf@SV  \POL@sturmchain@SV
     \let\POL@isolz@plusinf@sign\POL@sturmchain@sign
-    \POL@sturmchain@getSV@at{-1[\POL@isolz@E]}%
+    \POL@sturmchain@getSV@at{-1[\POL@isolz@E@neg]}%
     \let\POL@isolz@minusinf@SV  \POL@sturmchain@SV
     \let\POL@isolz@minusinf@sign\POL@sturmchain@sign
     \ifnum\POL@isolz@plusinf@sign=\z@
         \PackageError{polexpr}%
-{The polynomial #2 vanishes at set upper bound 10^\POL@isolz@E}%
-{Compile again with a bigger exponent in source. (X to abort).}%
+{The polynomial #2 vanishes at set upper bound 10^\POL@isolz@E@pos}%
+{Try again with a larger exponent. (X to abort).}%
     \fi
     \ifnum\POL@isolz@minusinf@sign=\z@
         \PackageError{polexpr}%
-{The polynomial #2 vanishes at set lower bound -10^\POL@isolz@E}%
-{Compile again with a bigger exponent in source. (X to abort).}%
+{The polynomial #2 vanishes at set lower bound -10^\POL@isolz@E@neg}%
+{Try again with a larger exponent. (X to abort).}%
     \fi
    \fi
    \edef\POL@isolz@NbOfRoots
@@ -792,7 +772,8 @@
 }%
 \def\POL@initarray#1#2{%
 % ATTENTION, if only one item, \xintAssignArray UNBRACES IT
-% so we use an \empty trick to avoid that. Maybe considered a bug of xinttools?
+% (is this to be considered as a bug of xinttools?)
+% We use an \empty trick to avoid that.
   \expandafter\xintAssignArray\expandafter\empty
     \romannumeral\xintreplicate{\POL@isolz@NbOfRoots}{{#2}}\to#1%
 }%
@@ -836,9 +817,18 @@
     \advance\POL@count\@ne
   \repeat
 }%
-% 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}%
+% This utility macro bounds positive roots (strictly) by a 10^Epos
+% and negative roots strictly by some -10^Eneg.
+% (prior to 0.8.6, an E was found with -10^E < all roots < 10^E)
+% To obtain Epos, the Cauchy bound  "1 + max_j {-a_j/lc(P)|}"
+% is used, where non-negative a_j/lc(P)'s are ignored. 
+% In case the a_j's all have same sign as lc(P) or vanish, there are
+% no positive roots.  And the macro in this case outputs an E=0 exponent.
+% But if at least one non-zero a_j has opposite sign to the leading coeff,
+% the produced E will be at least 1.
+% Thus if E=0 on exit, it is proof that there are no (positive) roots.
+\def\POL@isolz@updateE #1;%
+  {\unless\ifnum#1<\POL@isolz@E\space\edef\POL@isolz@E{\the\numexpr#1+\@ne}\fi}%
 \def\POL@isolz@getaprioribound{%
   \PolAssign{\POL@sturmname _0}\toarray\POL@arrayA
   \edef\POL@isolz@leading{\POL@arrayA{\POL@arrayA{0}}}%
@@ -849,19 +839,47 @@
      \expandafter\edef\csname POL@arrayA\the\POL@count\endcsname
         {\xintDiv{\POL@arrayA\POL@count}\POL@isolz@leading}%
   \repeat
-  \def\POL@isolz@E{1}% WE SEEK SMALLEST E SUCH HAT -10^E < roots < +10^E
+% We want an E such that 0 < positive roots < +10^E
+  \def\POL@isolz@E{0}%
   \advance\POL@count\m@ne
   \xintloop
   \ifnum\POL@count>\z@
-     \expandafter\POL@isolz@updateE
-     % use floating point to get decimal exponent
-     \romannumeral0\xintfloat[4]% should I use with [2] rather? (should work)
-       {\xintAdd{1/1[0]}{\xintAbs{\POL@arrayA\POL@count}}};%
+% only those coefficients with opposite sign to the leading coefficient
+% trigger an E update
+     \xintiiifSgn{\POL@arrayA\POL@count}%
+                 {\expandafter\POL@isolz@updateE
+     \the\numexpr\xintilogten{\xintAdd{1/1[0]}{\xintiiOpp{\POL@arrayA\POL@count}}};%
+                 }{}{}%
   \advance\POL@count\m@ne
   \repeat
-  % \ifxintverbose\xintMessage{polexpr}{Info}%
-  %   {Roots a priori bounded in absolute value by 10 to the \POL@isolz@E.}%
-  % \fi
+  \let\POL@isolz@E@pos\POL@isolz@E
+% We want an E such that 0 > negative roots > -10^E
+  \def\POL@isolz@E{0}%
+  \POL@count\POL@arrayA{0}\relax
+  \advance\POL@count\m@ne
+  \xintloop
+  \ifnum\POL@count>\@ne
+     \xintiiifSgn{\xintiiOpp{\POL@arrayA\POL@count}}%
+                 {\expandafter\POL@isolz@updateE
+     \the\numexpr\xintilogten{\xintAdd{1/1[0]}{\POL@arrayA\POL@count}};%
+                 }{}{}%
+     \advance\POL@count\m@ne
+     \xintiiifSgn{\POL@arrayA\POL@count}%
+                 {\expandafter\POL@isolz@updateE
+     \the\numexpr\xintilogten{\xintAdd{1/1[0]}{\xintiiOpp{\POL@arrayA\POL@count}}};%
+                 }{}{}%
+     \advance\POL@count\m@ne
+  \repeat
+  \ifnum\POL@count=\@ne
+     \xintiiifSgn{\xintiiOpp{\POL@arrayA\POL@count}}%
+                 {\expandafter\POL@isolz@updateE
+     \the\numexpr\xintilogten{\xintAdd{1/1[0]}{\POL@arrayA\POL@count}};%
+                 }{}{}%
+  \fi
+  \let\POL@isolz@E@neg\POL@isolz@E
+  \ifxintverbose
+     \xintMessage{polexpr}{Info}{Epos=\POL@isolz@E@pos, Eneg=\POL@isolz@E@neg.}%
+  \fi
 }%
 \def\POL@IsoRight@raw{\POL@IsoRight@Int/1[\POL@isolz@E]}%
 \def\POL@IsoLeft@raw {\POL@IsoLeft@Int/1[\POL@isolz@E]}%
@@ -875,8 +893,6 @@
        {\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.
   \global\POL@isolz@nextwillneedrefinefalse
   \def\POL@IsoRight@Int{0}%
   \POL@sturmchain@getSV@at\POL@IsoRight@raw
@@ -899,12 +915,13 @@
   % \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}%
-  \let\POL@isolz@@E\POL@isolz@E
+% 0.8.6 has separate initial E's for positive and negative roots
+  \let\POL@isolz@E\POL@isolz@E@neg
   \ifnum\POL@isolz@NbOfNegRoots>\z@
-% refactored at 0.7 to fix cases leading to an intervals with zero as end-point
+% refactored at 0.7 to fix cases leading to intervals having zero as end-point
     \POL@isolz@findroots@neg
   \fi
-  \let\POL@isolz@E\POL@isolz@@E
+  \let\POL@isolz@E\POL@isolz@E@pos
   \def\POL@IsoLeft@Int{0}%
   \let\POL@IsoLeftSV  \POL@IsoAtZeroSV  % véritable SV en zéro
   \let\POL@IsoLeftSign\POL@IsoAtZeroSign% véritable signe en zéro
@@ -973,10 +990,6 @@
   \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
@@ -985,7 +998,7 @@
     % 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 now 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
@@ -1101,14 +1114,6 @@
   \ifnum\POL@IsoRightSign=\POL@@IsoRightSign\space
   \repeat
   % now second root has been separated from the one at left end point
-% we update the storage of the root at left for it to have the same number
-% of digits in mantissa. No, I decided not to do that to avoid complications.
-  % \begingroup
-  %   \let\POL@IsoRight@Int\POL@IsoLeft@Int
-  %   \def\POL@IsoRightSign{0}%
-  %   \edef\POL@isolz@IntervalIndex{\the\numexpr\POL@isolz@IntervalIndex-\@ne}%
-  %   \POL@refine@storeleftandright
-  % \endgroup
   \edef\POL@@IsoRight@Int{\xintDSL{\xintInc{\xintDSR{\POL@IsoLeft@Int}}}}%
   \let\POL@IsoLeft@Int\POL@IsoRight@Int
   \let\POL@IsoLeftSign\POL@IsoRightSign
@@ -1128,7 +1133,7 @@
       % the IsoRightSign is now wrong but here we don't care
     \fi\fi
   \fi
-  % on exit, exact root found iff \POL@IsoRightSign is zero
+  % on exit, exact root has been found iff \POL@IsoRightSign is zero
   \POL@refine@storeleftandright
 \endgroup
 }%
@@ -1224,7 +1229,6 @@
     \fi
     \begingroup\xintglobaldefstrue
     % skip some overhead of \xintdefvar...
-    % Let me repeat: ATTENTION to change of internal format at xint 1.4
       \XINT_expr_defvar_one{\POL@sturmname L_\POL@isolz@IntervalIndex}%
          {{\POL@IsoLeft@rawout}}%
       \XINT_expr_defvar_one{\POL@sturmname R_\POL@isolz@IntervalIndex}%
@@ -1235,7 +1239,6 @@
     \endgroup
 }%
 %% \PolRefineInterval
-%% ATTENTION TO xint 1.4 INTERNAL CHANGES
 \def\POL@xintexprGetVar#1{\expandafter\expandafter\expandafter\xint_firstofone
                           \csname XINT_expr_varvalue_#1\endcsname}%
 % attention, also used by \POL@findrat@loop@a
@@ -1346,9 +1349,6 @@
 %
 %
 \def\PolIntervalWidth#1#2{%
-% le \xintRez est à cause des E positifs, car trailing zéros explicites
-% si je travaillais à partir des variables xintexpr directement ne devrait
-% pas être nécessaire, mais trop fragile par rapport à chgt internes possibles
      \romannumeral0\xintrez{\xintSub{\@nameuse{POL_ZR#1*}{#2}}%
                                     {\@nameuse{POL_ZL#1*}{#2}}}
 }%
@@ -1363,7 +1363,8 @@
 }%
 \def\POL@ensureintervallengths{%
    \POL@count\z@
-   % \POL@count used by \POL@sturmchain@getSV@at but latter not used
+   % attention that \POL@count would be modified by \POL@sturmchain@getSV@at
+   % but this latter macro not invoked by \POL@ensure@one
    \xintloop
      \advance\POL@count\@ne
      \edef\POL@isolz@IntervalIndex{\the\POL@count}%
@@ -1377,11 +1378,8 @@
    \edef\POL@sturmname{#1}%
    \edef\POL@ensure@targetE{\the\numexpr#3}%
    \edef\POL@isolz@IntervalIndex{\the\numexpr#2}%
-% 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_ZL\POL@sturmname*0\endcsname>\z@
-% je ne fais pas les \expandafter mais je préfèrerais ne pas être à l'intérieur
      \POL@ensure@one
    \fi
    \fi
@@ -1438,8 +1436,9 @@
 \catcode`_ 8 %
 \catcode`& 4 %
 \def\PolPrintIntervals{\POL@ifstar{\PolPrintIntervals@@}{\PolPrintIntervals@}}%
-% As explained in the docs, this is an example of customization so is not
-% itself customizable, apart from redefining it entirely!
+% As explained in the docs, the starred version is an example of customization
+% It is itself basically not easily customizable, except for this:
+\def\PolPrintIntervals@@arraystretch{2}% (the 2 was hardcoded prior to 0.8.6)
 \def\PolPrintIntervals@@{%
  \begingroup
    \def\POL@AfterPrintIntervals{\endgroup}%
@@ -1447,14 +1446,14 @@
    \let\PolPrintIntervalsUnknownRoot\POL@@PrintIntervalsUnknownRoot
    \let\PolPrintIntervalsKnownRoot\POL@@PrintIntervalsKnownRoot
 \ifdefined\array
-   \def\arraystretch{2}%
+   \let\arraystretch\PolPrintIntervals@@arraystretch
    \def\PolPrintIntervalsBeginEnv{\[\begin{array}{cl}}%\]
    \def\PolPrintIntervalsEndEnv{\end{array}\]}%
 \else
    \def\PolPrintIntervalsBeginEnv{$$\tabskip0pt plus 1000pt minus 1000pt
         \halign to\displaywidth\bgroup
-        \hfil\vrule height 2\ht\strutbox
-                    depth  2\dp\strutbox
+        \hfil\vrule height \PolPrintIntervals@@arraystretch\ht\strutbox
+                    depth  \PolPrintIntervals@@arraystretch\dp\strutbox
                     width  \z@
         $####$\tabskip6pt&$####$\hfil
         \tabskip0pt plus 1000pt minus 1000pt\cr}%$$
@@ -1476,6 +1475,10 @@
      \begingroup\edef\POL@tmp{\endgroup
        \unexpanded\expandafter{\PolPrintIntervalsBeginEnv}%
        \unexpanded\expandafter{\POL@PrintIntervals@Loop}%
+% This is added at 0.8.6 to allow usage of amsmath environment as they typeset
+% twice: we must prepare for a second execution.  Adds slight general overhead.
+       \gdef\noexpand\PolPrintIntervalsTheIndex{1}%
+       \noexpand\POL@PrintIntervals@DoDefs
        \unexpanded\expandafter{\PolPrintIntervalsEndEnv}%
      }\POL@tmp
    \fi
@@ -1485,7 +1488,7 @@
 }%
 \let\POL@AfterPrintIntervals\empty
 \let\PolPrintIntervalsNoRealRoots\empty
-\def\PolPrintIntervalsArrayStretch{1}%
+\def\PolPrintIntervalsArrayStretch{1}%   used by non-starred version
 \ifdefined\array
     \def\PolPrintIntervalsBeginEnv{\[\begin{array}{rcccl}}%
     \def\PolPrintIntervalsEndEnv{\end{array}\]}%
@@ -1558,7 +1561,6 @@
 }%
 \catcode`& 7 %
 \catcode`_ 11 %
-\def\POL@PrintIntervals@Loop#1{%
 \def\POL@PrintIntervals@Loop{%
   \POL@SturmIfZeroExactlyKnown\PolPrintIntervalsTheSturmName
                               \PolPrintIntervalsTheIndex
@@ -1568,10 +1570,15 @@
   \unless\ifnum\PolPrintIntervalsTheIndex>
                \@nameuse{POL_ZL\PolPrintIntervalsTheSturmName*0}
      \POL@PrintIntervals@DoDefs
-     \xint_afterfi{#1\POL@PrintIntervals@Loop}%
+     \xint_afterfi{\PolPrintIntervalsRowSeparator\POL@PrintIntervals@Loop}%
   \fi
-}}%
-\ifdefined\array\POL@PrintIntervals@Loop{\\}\else\POL@PrintIntervals@Loop{\cr}\fi
+}%
+% added at 0.8.6:
+\ifdefined\array
+  \def\PolPrintIntervalsRowSeparator{\\}%
+\else
+  \def\PolPrintIntervalsRowSeparator{\cr}%
+\fi
 \def\POL@PrintIntervals@DoDefs{%
   \xdef\PolPrintIntervalsTheLeftEndPoint{%
      \csname POL_ZL\PolPrintIntervalsTheSturmName*\PolPrintIntervalsTheIndex