diff options
Diffstat (limited to 'Master/texmf-dist/tex/generic/polexpr')
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 |