diff options
author | Karl Berry <karl@freefriends.org> | 2018-02-16 22:30:09 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2018-02-16 22:30:09 +0000 |
commit | 93ccafaede54490e77718a7d682343e94af262ec (patch) | |
tree | d0d24bd69b8b09fcd62656c4e31cb595b3446061 /Master/texmf-dist/tex | |
parent | 0f189d0ebeac95ecc1f5dab8d70d1a37633b8423 (diff) |
polexpr (16feb18)
git-svn-id: svn://tug.org/texlive/trunk@46653 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/tex')
-rw-r--r-- | Master/texmf-dist/tex/latex/polexpr/polexpr.sty | 1013 |
1 files changed, 946 insertions, 67 deletions
diff --git a/Master/texmf-dist/tex/latex/polexpr/polexpr.sty b/Master/texmf-dist/tex/latex/polexpr/polexpr.sty index 752efa3fa77..d4d45177744 100644 --- a/Master/texmf-dist/tex/latex/polexpr/polexpr.sty +++ b/Master/texmf-dist/tex/latex/polexpr/polexpr.sty @@ -1,14 +1,20 @@ % author: Jean-François Burnol % License: LPPL 1.3c (author-maintained) \ProvidesPackage{polexpr}% - [2018/01/18 v0.3.1 Polynomial expressions with rational coefficients (JFB)]% + [2018/02/16 v0.4 Polynomial expressions with rational coefficients (JFB)]% \RequirePackage{xintexpr}[2016/03/19]% xint 1.2g (or 1.2c 2015/11/16 at least) \edef\POL@restorecatcodes {\catcode`\noexpand\_ \the\catcode`\_ \catcode0 \the\catcode0\relax}% \catcode`\_ 11 \catcode0 12 %% AUXILIARIES -\newif\ifPOL@pol % (cf core algebra macros) +\newcount\POL@count +\newif\ifPOL@pol +\newif\ifxintveryverbose +\newif\ifpoltypesetall +\newif\ifPOL@sturm@normalize +\newif\ifPOL@isolz@nextwillneedrefine +\newif\ifpoltoexprall %% the main exchange structure (stored in macros \POLuserpol@<name>) %% is: degree.\empty{coeff0}{coeff1}....{coeffN} %% (degree=N except zero polynomial recognized from degree set to -1 @@ -22,19 +28,18 @@ {\def#3{#1}\expandafter\def\expandafter#4\expandafter{#2}}% % \def\POL@resultfromarray #1{% ATTENTION, **MUST** be executed with -% \count@ set to 1 + degree (thus \count@ = 0 for zero polynomial) - \edef\POL@result{% +% \count@ set to 1 + degree (\count@ = 0 for zero polynomial) + \edef\POL@result{\ifnum\count@>\z@ \the\numexpr\count@-\@ne.\noexpand\empty \xintiloop [1+1]% -% always done at least once with index 1, hence ok for zero polynomial \expandafter\POL@braceit\csname POL@array#1\xintiloopindex\endcsname \ifnum\xintiloopindex<\count@ - \repeat}% + \repeat + \else-1.\noexpand\empty{0/1[0]}\fi}% }% \def\POL@braceit#1{{#1}}% needed as \xintiloopindex can not "see" through braces -\newif\ifxintveryverbose \newcommand\PolDef[3][x]{\poldef #2(#1):=#3;}% \def\poldef{\edef\POL@restoresemicolon{\catcode59=\the\catcode59\relax}% \catcode59 12 \POL@defpol}% @@ -76,17 +81,23 @@ %% \def\POL@newpol#1{% \expandafter\POL@ifZero\csname POLuserpol@#1\endcsname - {\@namedef{XINT_expr_userfunc_#1}##1,{0/1[0]}% - \@namedef{XINT_flexpr_userfunc_#1}##1,{0[0]}}% + {\@namedef{XINT_expr_userfunc_#1}##1,{0/1[0]}}% {\POL@newpolhorner{#1}}% \unless\ifcsname XINT_expr_userfuncNE:#1\endcsname\POL@addtoextras{#1}\fi \expandafter\XINT_expr_defuserfunc \csname XINT_expr_func_#1\expandafter\endcsname \csname XINT_expr_userfunc_#1\endcsname + \expandafter\let\csname XINT_flexpr_func_#1\endcsname\@undefined + \ifxintverbose\POL@info{#1}\fi +}% +\def\POL@newfloatpol#1{% + \expandafter\POL@ifZero\csname POLuserpol@#1\endcsname + {\@namedef{XINT_flexpr_userfunc_#1}##1,{0[0]}}% + {\POL@newfloatpolhorner{#1}}% \expandafter\XINT_expr_defuserfunc \csname XINT_flexpr_func_#1\expandafter\endcsname \csname XINT_flexpr_userfunc_#1\endcsname - \ifxintverbose\POL@info{#1}\fi + \ifxintverbose\POL@floatinfo{#1}\fi }% \def\POL@info #1{% \xintMessage {polexpr}{Info}% @@ -94,9 +105,13 @@ associated to \string\XINT_expr_userfunc_#1\space whose meaning uses Horner scheme: \expandafter\meaning - \csname XINT_expr_userfunc_#1\endcsname^^J% - \@spaces And \string\XINT_flexpr_userfunc_#1\space for - the \string\xintfloatexpr\space parser has meaning + \csname XINT_expr_userfunc_#1\endcsname}% +}% +\def\POL@floatinfo #1{% + \xintMessage {polexpr}{Info}% + {Function #1 for the \string\xintfloatexpr\space parser is + associated to \string\XINT_flexpr_userfunc_#1\space + whose meaning uses Horner scheme: \expandafter\meaning \csname XINT_flexpr_userfunc_#1\endcsname}% }% @@ -114,6 +129,14 @@ \endgroup \expandafter\def\csname XINT_expr_userfunc_#1\expandafter\endcsname \expandafter##\expandafter1\expandafter,\expandafter{\POL@tmp{##1}}% +}% +\def\POL@newfloatpolhorner#1{% + %% redefine function to expand by Horner scheme. Is this useful? + %% perhaps bad idea for numerical evaluation of thing such as (1+x)^10? +% note: I added {0/1[0]} item to zero polynomial also to facilitate this + \expandafter\expandafter\expandafter\POL@split + \csname POLuserpol@#1\endcsname;\POL@var@deg\POL@var@coeffs + \edef\POL@var@coeffs{\xintRevWithBraces{\POL@var@coeffs}}% \begingroup \expandafter\POL@newpol@floathorner\POL@var@coeffs\relax \expandafter @@ -177,10 +200,14 @@ %\let\POL@original@redefinemacros\XINT_expr_redefinemacros % do locally \def\POL@redefinemacros{\POL@original_redefinemacros\POL@redefineextras}% \let\POL@redefineextras\@empty -% + + +\newcommand\PolGenFloatVariant[1]{\POL@newfloatpol{#1}}% + + \newcommand\PolLet[2]{\if=\noexpand#2\expandafter\xint_firstoftwo \else\expandafter\xint_secondoftwo\fi - \POL@@let\POL@let {#1}{#2}}% + \POL@@let\POL@let{#1}{#2}}% \def\POL@@let#1#2#3{\POL@let{#1}{#3}}% \def\POL@let#1#2{% \expandafter\let\csname POLuserpol@#1\expandafter\endcsname @@ -188,17 +215,20 @@ \unless\ifcsname XINT_expr_userfuncNE:#1\endcsname\POL@addtoextras{#1}\fi \expandafter\let\csname XINT_expr_userfunc_#1\expandafter\endcsname \csname XINT_expr_userfunc_#2\endcsname - \expandafter\let\csname XINT_flexpr_userfunc_#1\expandafter\endcsname - \csname XINT_flexpr_userfunc_#2\endcsname \expandafter\XINT_expr_defuserfunc \csname XINT_expr_func_#1\expandafter\endcsname \csname XINT_expr_userfunc_#1\endcsname - \expandafter\XINT_expr_defuserfunc - \csname XINT_flexpr_func_#1\expandafter\endcsname - \csname XINT_flexpr_userfunc_#1\endcsname \ifxintverbose\POL@info{#1}\fi }% - +\newcommand\PolGlobalLet[2]{\begingroup + \globaldefs\@ne + \if=\noexpand#2\expandafter\xint_firstoftwo + \else\expandafter\xint_secondoftwo\fi +% There is a potential problem related to \POL@addtoextras, (the local set-up +% will become the global one) but I will reconsider that another day + \POL@@globallet\POL@globallet {#1}{#2}}% +\def\POL@@globallet#1#2#3{\POL@globallet{#1}{#3}}% +\def\POL@globallet#1#2{\POL@let{#1}{#2}\endgroup}% \newcommand\PolAssign[1]{\def\POL@polname{#1}\POL@assign}% zap spaces in #1? \def\POL@assign#1\toarray#2{% @@ -224,13 +254,15 @@ }% -\newcommand\PolGet[1]{\def\POL@polname{#1}% zap spaces in #1? - \begingroup % closed in \POL@getfrom - \POL@getfrom}% -% attention au name clash proche avec \POL@get auxiliaire de \POL@add etc.. -\def\POL@getfrom#1\fromarray#2{% - \count@#2{0} % must be > 0, else could create infinite loop - % \ifnum\count@>\z@\else\InvalidArrayError_PolGet\fi +\newcommand\PolGet{} +\def\PolGet#1#2\fromarray#3{% + \begingroup % closed in \POL@getfromarray + \POL@getfromarray{#1}{#3}% + \POL@newpol{#1}% +}% +\def\POL@getfromarray#1#2{% + \count@=#2{0} %<- intentional space, + % must be > 0, else could create infinite loop \xintloop \edef\POL@tmp{#2{\count@}}% \xintiiifZero{\POL@tmp}% @@ -251,27 +283,30 @@ \expandafter \endgroup \expandafter - \def\csname POLuserpol@\POL@polname\expandafter\endcsname + \def\csname POLuserpol@#1\expandafter\endcsname \expandafter{\POL@result}% - \expandafter\POL@newpol\expandafter{\POL@polname}% }% -\newcommand\PolFromCSV[2]{\def\POL@polname{#1}% - \begingroup % closed in \POL@getfrom +\newcommand\PolFromCSV[2]{% + \begingroup % closed in \POL@getfromarray \xintAssignArray\xintCSVtoList{#2}\to\POL@arrayA - \POL@getfrom\fromarray\POL@arrayA + \POL@getfromarray{#1}\POL@arrayA + \POL@newpol{#1}% % semble un peu indirect et sous-optimal % mais je veux élaguer les coefficients nuls. Peut-être à revoir. }% -\newif\ifpoltypesetall \newcommand\PolTypesetCmdPrefix[1]{\xintiiifSgn{#1}{}{+}{+}}% \newcommand\PolTypesetCmd[1]{\xintifOne{\xintiiAbs{#1}}% {\ifnum\PolIndex=\z@\xintiiSgn{#1}\else - \xintiiifSgn{#1}{-}{}{}\fi}% - {\xintSignedFrac{#1}}}% + \xintiiifSgn{#1}{-}{}{}\fi + \let\PolIfCoeffIsPlusOrMinusOne\@firstoftwo}% + {\PolTypesetOne{#1}% + \let\PolIfCoeffIsPlusOrMinusOne\@secondoftwo}% + }% +\newcommand\PolTypesetOne{\xintSignedFrac}% \newcommand\PolTypesetMonomialCmd{% \ifcase\PolIndex\space % @@ -330,6 +365,10 @@ \newcommand\PolMapCoeffs[2]{% #1 = macro, #2 = name + \POL@mapcoeffs{#1}{#2}% + \POL@newpol{#2}% +}% +\def\POL@mapcoeffs#1#2{% \begingroup \def\POL@map@macro{#1}% \expandafter\expandafter\expandafter\POL@split @@ -352,7 +391,6 @@ \endgroup \expandafter \def\csname POLuserpol@#2\expandafter\endcsname\expandafter{\POL@result}% - \POL@newpol{#2}% }% \def\POL@map@loop.#1{\if\relax#1\expandafter\xint_gob_til_dot\fi \advance\count@\@ne @@ -362,11 +400,34 @@ \edef\index{\the\numexpr\index+\@ne}% \POL@map@loop.}% \def\POL@xintIrr#1{\xintIrr{#1}[0]}% -\newcommand\PolReduceCoeffs[1]{\PolMapCoeffs{\POL@xintIrr}{#1}}% +\def\POL@special@xintIrr#1{% + \expandafter\POL@special@xintIrr@i\romannumeral0\xintraw{#1}}% +\def\POL@special@xintIrr@i#1/#2[#3]{\xintIrr{#1/#2[0]}[#3]}% +\newcommand\PolReduceCoeffs{\@ifstar\POL@sreducecoeffs\POL@reducecoeffs}% +\def\POL@reducecoeffs#1{\PolMapCoeffs{\POL@xintIrr}{#1}}% +\def\POL@sreducecoeffs#1{\PolMapCoeffs{\POL@special@xintIrr}{#1}}% +\def\POL@special@xintIrr@skipraw#1{\POL@special@xintIrr@i#1}% used by ToSturm %% EUCLIDEAN DIVISION \newcommand\PolDivide[4]{% #3=quotient, #4=remainder of #1 by #2 + \POL@divide{#1}{#2}% + \expandafter\let\csname POLuserpol@#3\endcsname\POL@Q + \POL@newpol{#3}% + \expandafter\let\csname POLuserpol@#4\endcsname\POL@R + \POL@newpol{#4}% +}% +\newcommand\PolQuo[3]{% #3=quotient of #1 by #2 + \POL@divide{#1}{#2}% + \expandafter\let\csname POLuserpol@#3\endcsname\POL@Q + \POL@newpol{#3}% +}% +\newcommand\PolRem[3]{% #3=remainder of #1 by #2 + \POL@divide{#1}{#2}% + \expandafter\let\csname POLuserpol@#3\endcsname\POL@R + \POL@newpol{#3}% +}% +\newcommand\POL@divide[2]{% \begingroup \let\xintScalarSub\xintSub \let\XINT_fadd_C\POL_fadd_C @@ -386,17 +447,19 @@ \expandafter \endgroup \expandafter - \def\csname POLuserpol@#3\expandafter\expandafter\expandafter\endcsname + \def\csname POL@Q\expandafter\expandafter\expandafter\endcsname \expandafter\expandafter\expandafter{\expandafter\POL@Q\expandafter}% \expandafter - \def\csname POLuserpol@#4\expandafter\endcsname\expandafter{\POL@R}% - \POL@newpol{#3}% - \POL@newpol{#4}% + \def\csname POL@R\expandafter\endcsname\expandafter{\POL@R}% }% %% GCD \newcommand\PolGCD[3]{% sets #3 to the (unitary) G.C.D. of #1 and #2 + \POL@GCD{#1}{#2}{#3}% + \POL@newpol{#3}% +}% +\def\POL@GCD #1#2#3{% \begingroup \let\xintScalarSub\xintSub \let\XINT_fadd_C\POL_fadd_C @@ -424,7 +487,8 @@ \POL@gcd@exit AB}% {\ifnum\POL@degA<\POL@degB\space \let\POL@tmp\POL@B\let\POL@B\POL@A\let\POL@A\POL@tmp - \let\POL@tmp\POL@degB\let\POL@degB\POL@degA\let\POL@degA\POL@degB + \let\POL@tmp\POL@degB\let\POL@degB\POL@degA\let\POL@degA\POL@tmp + \let\POL@tmp\POL@polB\let\POL@polB\POL@polA\let\POL@polA\POL@tmp \fi \xintAssignArray\POL@polA\to\POL@arrayA \xintAssignArray\POL@polB\to\POL@arrayB @@ -434,7 +498,6 @@ \endgroup \expandafter\def\csname POLuserpol@#3\expandafter\endcsname \expandafter{\POL@result}% - \POL@newpol{#3}% }% \def\POL@normalize#1{% \expandafter\def\expandafter\POL@tmp\expandafter @@ -450,7 +513,7 @@ {\POL@normalize@leading}}[0]}% \advance\count@\m@ne \repeat -}% +}% \def\POL@gcd#1#2{% \POL@normalize{#2}% \edef\POL@degQ{\the\numexpr\csname POL@deg#1\endcsname @@ -498,7 +561,7 @@ \count@\numexpr\csname POL@deg#1\endcsname+\@ne\relax \POL@resultfromarray #1% }% - + %% TODO: BEZOUT @@ -533,7 +596,8 @@ \fi }% \def\POL@Diff@no #1#2{\POL@let{#2}{#1}}% -\def\POL@Diff@one #1#2{% +\def\POL@Diff@one #1#2{\POL@Diff@@one {#1}{#2}\POL@newpol{#2}}% +\def\POL@Diff@@one#1#2{% \expandafter\expandafter\expandafter\POL@split \csname POLuserpol@#1\endcsname;\POL@var@deg\POL@var@coeffs \ifnum\POL@var@deg<\@ne @@ -543,7 +607,6 @@ \expandafter\edef\csname POLuserpol@#2\endcsname {\the\numexpr\POL@var@deg-\@ne.\noexpand\empty\POL@var@coeffs}% \fi - \POL@newpol{#2}% }% % lazy way but allows to share with AntiDiff \def\POL@Iterate#1#2#3{% @@ -593,7 +656,8 @@ \fi }% \let\POL@AntiDiff@no\POL@Diff@no -\def\POL@AntiDiff@one #1#2{% +\def\POL@AntiDiff@one #1#2{\POL@AntiDiff@@one{#1}{#2}\POL@newpol{#2}}% +\def\POL@AntiDiff@@one#1#2{% \expandafter\expandafter\expandafter\POL@split \csname POLuserpol@#1\endcsname;\POL@var@deg\POL@var@coeffs \ifnum\POL@var@deg<\z@ @@ -603,7 +667,798 @@ \expandafter\edef\csname POLuserpol@#2\endcsname {\the\numexpr\POL@var@deg+\@ne.\noexpand\empty{0/1[0]}\POL@var@coeffs}% \fi - \POL@newpol{#2}% +}% + + +%% Sturm Algorithm (polexpr 0.4) +\newcommand\PolToSturm{% + \@ifstar{\POL@sturm@normalizefalse}{\POL@sturm@normalizetrue}% + \POL@ToSturm +}% +\newcommand\POL@ToSturm[2]{% + \edef\POL@sturmname{#2}% + \POL@let{\POL@sturmname _0}{#1}% + \POL@mapcoeffs{\POL@special@xintIrr}{\POL@sturmname _0}% + \POL@Diff@@one{\POL@sturmname _0}{\POL@sturmname _1}% + \POL@count\@ne + \xintloop + \POL@divide{\POL@sturmname _\the\numexpr\POL@count-\@ne}% + {\POL@sturmname _\the\POL@count}% + \expandafter\POL@split\POL@R;\POL@degR\POL@polR + \unless\ifnum\POL@degR=\m@ne + \edef\POL@polR{\xintApply{\POL@dooppandirr}{\POL@polR}}% + \advance\POL@count\@ne + \expandafter\edef\csname POLuserpol@\POL@sturmname _\the\POL@count\endcsname + {\POL@degR.\noexpand\empty\POL@polR}% + \repeat + \edef\POL@sturm@N{\the\POL@count}% + \ifPOL@sturm@normalize + \ifnum\PolDegree{\POL@sturmname _\POL@sturm@N}>\z@ + \xintloop + \advance\POL@count\m@ne + \POL@divide{\POL@sturmname _\the\POL@count}% + {\POL@sturmname _\POL@sturm@N}% + \expandafter\POL@split\POL@Q;\POL@degQ\POL@polQ + \edef\POL@polQ{\xintApply{\POL@special@xintIrr@skipraw}{\POL@polQ}}% + \edef\POL@Q{\POL@degQ.\noexpand\empty\POL@polQ}% + \expandafter\edef\csname POLuserpol@\POL@sturmname _\the\POL@count\endcsname + {\POL@degQ.\noexpand\empty\POL@polQ}% + \ifnum\POL@count>\z@ + \repeat + \@namedef{POLuserpol@\POL@sturmname _\POL@sturm@N}{0.\empty{1/1[0]}}% + \fi + \fi + \POL@count\z@ + \xintloop + \POL@newpol{\POL@sturmname _\the\POL@count}% + \unless\ifnum\POL@sturm@N=\POL@count + \advance\POL@count\@ne + \repeat + \expandafter\let\csname PolSturmChainLength_\POL@sturmname \endcsname\POL@sturm@N +}% +\def\POL@dooppandirr#1{\xintiiOpp{\POL@special@xintIrr@i#1}}% +\newcommand\PolSturmChainLength[1] + {\romannumeral`^^@\csname PolSturmChainLength_#1\endcsname}% +\newcommand\PolSetToSturmChainSignChangesAt[4][\global]{% + \edef\POL@sturmchain@X{\xintREZ{#4}}% + \edef\POL@sturmname{#3}% + \edef\POL@sturmlength{\PolSturmChainLength{\POL@sturmname}}% + \POL@sturmchain@getSV@at\POL@sturmchain@X + #1\let#2\POL@sturmchain@SV +}% +\def\POL@sturmchain@getSV@at#1{% ATTENTION USES \POL@count + \def\POL@sturmchain@SV{0}% + \edef\POL@isolz@lastsign{\xintiiSgn{\PolEval{\POL@sturmname _0}\At{#1}}}% + \let\POL@IsoRightSign\POL@isolz@lastsign % needed only for SturmIsolate etc... + \POL@count \z@ + \ifnum\POL@isolz@lastsign=\z@ + \edef\POL@isolz@lastsign + {\xintiiSgn{\PolEval{\POL@sturmname _1}\At{#1}}}% + \POL@count \@ne + \fi + \xintloop + \unless\ifnum\POL@sturmlength=\POL@count + \advance\POL@count \@ne + \edef\POL@isolz@newsign + {\xintiiSgn{\PolEval{\POL@sturmname _\the\POL@count}\At{#1}}}% + \ifnum\POL@isolz@newsign=\numexpr-\POL@isolz@lastsign\relax + \edef\POL@sturmchain@SV{\the\numexpr\POL@sturmchain@SV+\@ne}% + \let\POL@isolz@lastsign=\POL@isolz@newsign + \fi + \repeat +}% +\newcommand\PolSetToNbOfZerosWithin[5][\global]{% + \edef\POL@tmpA{\xintREZ{#4}}% + \edef\POL@tmpB{\xintREZ{#5}}% + \edef\POL@sturmname{#3}% + \edef\POL@sturmlength{\PolSturmChainLength{\POL@sturmname}}% + \POL@sturmchain@getSV@at\POL@tmpA + \let\POL@SVA\POL@sturmchain@SV + \POL@sturmchain@getSV@at\POL@tmpB + \let\POL@SVB\POL@sturmchain@SV + \ifnum\POL@SVA<\POL@SVB\space + #1\edef#2{\the\numexpr\POL@SVB-\POL@SVA}% + \else + #1\edef#2{\the\numexpr\POL@SVA-\POL@SVB}% + \fi +}% + + +\newcommand\PolSturmIsolateZeros[1]{% + % #1 name of Sturm chain (already pre-computed from a given polynomial) + \edef\POL@sturmname{#1}% + \edef\POL@sturmlength{\PolSturmChainLength{#1}}% + % Count number of sign changes at plus infinity in Sturm sequence + \def\POL@isolz@plusinf@SV{0}% + \edef\POL@isolz@lastsign{\xintiiSgn{\PolLeadingCoeff{#1_0}}}% + \let\POL@isolz@plusinf@sign\POL@isolz@lastsign + \POL@count\@ne + \xintloop + \edef\POL@isolz@newsign{\xintiiSgn{\PolLeadingCoeff{#1_\the\POL@count}}}% + \unless\ifnum\POL@isolz@newsign=\POL@isolz@lastsign + \edef\POL@isolz@plusinf@SV{\the\numexpr\POL@isolz@plusinf@SV+\@ne}% + \fi + \let\POL@isolz@lastsign=\POL@isolz@newsign + \ifnum\POL@sturmlength>\POL@count + \advance\POL@count\@ne + \repeat + % Count number of sign changes at minus infinity in Sturm sequence + \def\POL@isolz@minusinf@SV{0}% + \edef\POL@isolz@lastsign{\xintiiSgn{\PolLeadingCoeff{#1_0}}}% + \ifodd\PolDegree{#1_0} + \edef\POL@isolz@lastsign{\xintiiOpp{\POL@isolz@lastsign}}% + \fi + \let\POL@isolz@minusinf@sign\POL@isolz@lastsign + \POL@count\@ne + \xintloop + \edef\POL@isolz@newsign{\xintiiSgn{\PolLeadingCoeff{#1_\the\POL@count}}}% + \ifodd\PolDegree{#1_\the\POL@count} + \edef\POL@isolz@newsign{\xintiiOpp{\POL@isolz@newsign}}% + \fi + \unless\ifnum\POL@isolz@newsign=\POL@isolz@lastsign + \edef\POL@isolz@minusinf@SV{\the\numexpr\POL@isolz@minusinf@SV+\@ne}% + \fi + \let\POL@isolz@lastsign=\POL@isolz@newsign + \ifnum\POL@sturmlength>\POL@count + \advance\POL@count\@ne + \repeat + \edef\POL@isolz@NbOfRoots + {\the\numexpr\POL@isolz@minusinf@SV-\POL@isolz@plusinf@SV}% + \ifnum\POL@isolz@NbOfRoots=\z@ + \begingroup\globaldefs\@ne + \expandafter\xintAssignArray\expandafter\to\csname POL_ZeroInt#1L\endcsname + \expandafter\xintAssignArray\expandafter\to\csname POL_ZeroInt#1R\endcsname + \endgroup + \else + \begingroup\globaldefs\@ne + \expandafter\POL@isolz@initarray\csname POL_ZeroInt#1L\endcsname + \expandafter\POL@isolz@initarray\csname POL_ZeroInt#1R\endcsname + \endgroup + \expandafter\POL@isolz@getaprioribound + \fi +}% +\def\POL@isolz@initarray#1{% + \expandafter\xintAssignArray + \romannumeral\xintreplicate{\POL@isolz@NbOfRoots}{{0}}\to#1% +}% +% 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}% +\def\POL@isolz@getaprioribound{% + \PolAssign{\POL@sturmname _0}\toarray\POL@arrayA + \edef\POL@isolz@leading{\POL@arrayA{\POL@arrayA{0}}}% + \POL@count\z@ + \xintloop + \advance\POL@count\@ne + \ifnum\POL@arrayA{0}>\POL@count + \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 + \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}}};% + \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 + \POL@isolz@main +}% +\def\POL@IsoRight@raw{\POL@IsoRight@Int/1[\POL@isolz@E]}% +\def\POL@IsoLeft@raw {\POL@IsoLeft@Int/1[\POL@isolz@E]}% +\def\POL@IsoRight@rawout{% + \ifnum\POL@IsoRightSign=\z@\expandafter\xintREZ\fi\POL@IsoRight@raw}% +\def\POL@IsoLeft@rawout{% + \ifnum\POL@IsoRightSign=\z@ + \expandafter\xint_firstoftwo\else\expandafter\xint_secondoftwo + \fi{\xintREZ\POL@IsoRight@raw}% + {\POL@IsoLeft@Int/1[\POL@isolz@E]}}% +\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 + \let\POL@IsoRightSV\POL@sturmchain@SV + \let\POL@IsoAtZeroSV\POL@IsoRightSV + \let\POL@IsoAtZeroSign\POL@IsoRightSign + \ifnum\POL@IsoAtZeroSign=\z@ + \xdef\POL@isolz@IntervalIndex + {\the\numexpr\POL@isolz@minusinf@SV-\POL@IsoRightSV}% + \POL@refine@storeleftandright % store zero root + \edef\POL@IsoRightSV{\the\numexpr\POL@IsoRightSV+\@ne}% +% subtlety here if original polynomial had multiplicities, but ok. I checked! + \edef\POL@IsoRightSign % evaluated twice, but that's not so bad + {\xintiiOpp{\xintiiSgn{\PolEval{\POL@sturmname _1}\At{0/1[0]}}}}% + \fi + \def\POL@IsoLeft@Int{-1}% -10^E isn't a root! + \let\POL@IsoLeftSV\POL@isolz@minusinf@SV + \let\POL@IsoLeftSign\POL@isolz@minusinf@sign + \edef\POL@isolz@NbOfNegRoots{\the\numexpr\POL@IsoLeftSV-\POL@IsoRightSV}% + \gdef\POL@isolz@IntervalIndex{0}% + \begingroup + \let\POL@IsoAtZeroSV\POL@IsoRightSV % locally shifted if root at zero + \let\POL@IsoAtZeroSign\POL@IsoRightSign + \ifnum\POL@isolz@NbOfNegRoots>\z@ + \def\POL@IsoRight@Int{-1}% + \xintloop + \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}% + \POL@sturmchain@getSV@at\POL@IsoRight@raw + \let\POL@IsoRightSV\POL@sturmchain@SV + % would an \ifx test be quicker? (to be checked) + \ifnum\POL@IsoRightSV=\POL@IsoLeftSV\space + % no roots in-between, sign and SV kept + \repeat + \def\POL@IsoLeft@Int{-10}% + \let\POL@@IsoRightSign\POL@IsoRightSign % zero possible + \let\POL@@IsoRightSV\POL@IsoRightSV + \xintloop + \edef\POL@IsoRight@Int{\the\numexpr\POL@IsoLeft@Int+\@ne}% +% we could arguably do a more efficient dichotomy here + \POL@sturmchain@getSV@at\POL@IsoRight@raw + \let\POL@IsoRightSV\POL@sturmchain@SV + \POL@isolz@check + \ifnum\POL@isolz@IntervalIndex=\POL@isolz@NbOfNegRoots\space + \expandafter\xintbreakloop + \fi + \let\POL@IsoLeft@Int\POL@IsoRight@Int + \let\POL@IsoLeftSign\POL@IsoRightSign + \let\POL@IsoLeftSV\POL@IsoRightSV + \ifnum\POL@IsoRight@Int < -\tw@ + \repeat + \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfNegRoots\space + \def\POL@IsoRight@Int{-1}% + \let\POL@IsoRightSign\POL@@IsoRightSign + \let\POL@IsoRightSV\POL@@IsoRightSV + \POL@isolz@check + \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfNegRoots\space + \def\POL@IsoLeft@Int{-1}% + \let\POL@IsoLeftSign\POL@IsoRightSign + \let\POL@IsoLeftSV\POL@IsoRightSV + \def\POL@IsoRight@Int{0}% + \let\POL@IsoRightSV\POL@IsoAtZeroSV % altered if 0 was a root + \let\POL@IsoRightSign\POL@IsoAtZeroSign% id. +% this will recurse to locate roots with smaller decimal exponents + \POL@isolz@check % attention that this should not re-evaluate at 0 + \fi + \fi + \fi + \endgroup + \def\POL@IsoLeft@Int{0}% + \let\POL@IsoLeftSV\POL@IsoAtZeroSV + \let\POL@IsoLeftSign\POL@IsoAtZeroSign + \ifnum\POL@IsoLeftSign=\z@ + \xdef\POL@isolz@IntervalIndex{\the\numexpr\POL@isolz@IntervalIndex+\@ne}% + \global\POL@isolz@nextwillneedrefinetrue + \else + \global\POL@isolz@nextwillneedrefinefalse + \fi + \let\POL@@IsoRightSV=\POL@isolz@plusinf@SV + \let\POL@@IsoRightSign=\POL@isolz@plusinf@sign % 10^E not a root! + \edef\POL@isolz@NbOfPosRoots + {\the\numexpr\POL@IsoLeftSV-\POL@@IsoRightSV}% attention @@ + \ifnum\POL@isolz@NbOfPosRoots>\z@ + \def\POL@IsoRight@Int{1}% + \xintloop + \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}% + \POL@sturmchain@getSV@at\POL@IsoRight@raw + \let\POL@IsoRightSV\POL@sturmchain@SV + \ifnum\POL@IsoRightSV=\POL@@IsoRightSV\space + \let\POL@@IsoRightSign\POL@IsoRightSign % root here possible! + \repeat + \unless\ifnum\POL@IsoRightSV=\POL@IsoLeftSV\space + \POL@isolz@check % will recurse inside groups if needed + \fi + \def\POL@IsoLeft@Int{1}% + \let\POL@IsoLeftSV\POL@IsoRightSV + \let\POL@IsoLeftSign\POL@IsoRightSign + \xintloop +% we could arguably do a more efficient dichotomy here + \edef\POL@IsoRight@Int{\the\numexpr\POL@IsoLeft@Int+\@ne}% + \POL@sturmchain@getSV@at\POL@IsoRight@raw + \let\POL@IsoRightSV\POL@sturmchain@SV + \POL@isolz@check + \let\POL@IsoLeft@Int\POL@IsoRight@Int + \let\POL@IsoLeftSign\POL@IsoRightSign + \let\POL@IsoLeftSV\POL@IsoRightSV + \ifnum\POL@isolz@IntervalIndex=\POL@isolz@NbOfRoots\space + \expandafter\xintbreakloop + \fi + \ifnum\POL@IsoLeft@Int < \xint_c_ix + \repeat + \ifnum\POL@isolz@IntervalIndex<\POL@isolz@NbOfRoots\space + % get now the last, right most, root (or roots) + \def\POL@IsoRight@Int{10}% + \let\POL@IsoRightSign\POL@@IsoRightSign + \let\POL@IsoRightSV\POL@@IsoRightSV + \POL@isolz@check + \fi + \fi +}% +\def\POL@isolz@check{% \POL@IsoRightSign must be ready for use here +% \ifxintverbose +% \xintMessage{polexpr}{Info}% +% {\the\numexpr\POL@IsoLeftSV-\POL@IsoRightSV\relax\space roots +% in (\POL@IsoLeft@raw,\POL@IsoRight@raw] (E = \POL@isolz@E)}% +% \fi + \ifcase\numexpr\POL@IsoLeftSV-\POL@IsoRightSV\relax + % no root in ]left, right] + \global\POL@isolz@nextwillneedrefinefalse + \or + % exactly one root in ]left, right] + \xdef\POL@isolz@IntervalIndex{\the\numexpr\POL@isolz@IntervalIndex+\@ne}% + \ifnum\POL@IsoRightSign=\z@ + % if right boundary is a root, ignore previous flag + \global\POL@isolz@nextwillneedrefinefalse + \fi + % if left boundary is known to have been a root we refine interval + \ifPOL@isolz@nextwillneedrefine + \expandafter\expandafter\expandafter\POL@isolz@refine + \else + \POL@refine@storeleftandright + \ifnum\POL@IsoRightSign=\z@ + \global\POL@isolz@nextwillneedrefinetrue + \fi + \fi + \else + % more than one root, we need to recurse + \expandafter\POL@isolz@recursedeeper + \fi +}% +\def\POL@isolz@recursedeeper{% +% NOTE 2018/02/16. I SHOULD DO A REAL BINARY DICHOTOMY HERE WHICH ON AVERAGE +% SHOULD BRING SOME GAIN (LIKE WHAT IS ALREADY DONE FOR THE "refine" MACROS. +% THUS IN FUTURE THIS MIGHT BE REFACTORED. +\begingroup + \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}% + \edef\POL@@IsoRight@Int{\xintDSL{\POL@IsoRight@Int}}% + \let\POL@@IsoRightSign \POL@IsoRightSign + \let\POL@@IsoRightSV \POL@IsoRightSV + \edef\POL@IsoLeft@Int {\xintDSL{\POL@IsoLeft@Int}}% + \xintiloop[1+1] + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% + \POL@sturmchain@getSV@at\POL@IsoRight@raw + \let\POL@IsoRightSV\POL@sturmchain@SV + \POL@isolz@check + \let\POL@IsoLeft@Int\POL@IsoRight@Int + \let\POL@IsoLeftSV\POL@IsoRightSV + \let\POL@IsoLeftSign\POL@IsoRightSign% not used, actually + \ifnum\POL@IsoLeftSV=\POL@@IsoRightSV\space + \expandafter\xintbreakiloop + \fi + \ifnum\xintiloopindex < \xint_c_ix + \repeat + \let\POL@IsoRight@Int\POL@@IsoRight@Int + \let\POL@IsoRightSign\POL@@IsoRightSign + \let\POL@IsoRightSV \POL@@IsoRightSV + % if we exited the loop via breakiloop this is superfluous + % but it only costs one \ifnum + \POL@isolz@check +\endgroup +}% +\def\POL@isolz@refine{% + % starting point is first root = left < unique second root < right + % even if we hit exactly via refinement second root, we set flag false as + % processing will continue with original right end-point, which isn't a root + \global\POL@isolz@nextwillneedrefinefalse +\begingroup + \let\POL@@IsoRightSign\POL@IsoRightSign % already evaluated + \xintloop + \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}% + \edef\POL@IsoLeft@Int {\xintDSL{\POL@IsoLeft@Int}}% + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \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 + \ifnum\POL@IsoRightSign=\z@ % check if new Left is actually a root + \else + \edef\POL@IsoRight@Int{\xintDec{\POL@@IsoRight@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@@IsoRightSign\space + \POL@refine@doonce % we need to locate in interval (1, 9) in local scale + \else + \let\POL@IsoLeft@Int\POL@IsoRight@Int + \ifnum\POL@IsoRightSign=\z@ + \def\POL@IsoLeftSign{0}% + \else + \let\POL@IsoRight@Int\POL@@IsoRight@Int + % the IsoRightSign is now wrong but here we don't care + \fi\fi + \fi + % on exit, exact root found iff \POL@IsoRightSign is zero + \POL@refine@storeleftandright +\endgroup +}% +\def\POL@refine@doonce{% if exact root is found, always in IsoRight on exit +% NOTE: FUTURE REFACTORING WILL GET RID OF \xintiiAdd WHICH ARE A BIT COSTLY +% BUT BASICALLY NEEDED TO HANDLE BOTH NEGATIVE AND POSITIVE HERE. +% I WILL RE-ORGANIZE THE WHOLE THING IN FUTURE TO GET ROOTS STARTING FROM +% THE ORIGIN AND SIMPLY RE-LABEL THE NEGATIVE ONE AT THE END. 2018/02/16. + \let\POL@@IsoRight@Int\POL@IsoRight@Int % 9 + \let\POL@@IsoRightSign\POL@IsoRightSign + \edef\POL@IsoRight@Int{\xintiiAdd{4}{\POL@IsoLeft@Int}}% 5 + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int % 5 + \edef\POL@IsoRight@Int{\xintiiAdd{2}{\POL@IsoLeft@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int % 7 + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int % 8 + \let\POL@IsoRight@Int\POL@@IsoRight@Int % 9 + \let\POL@IsoRightSign\POL@@IsoRightSign % opposite of one at left + \fi % else 7, 8 with possible root at 8 + \else + \ifnum\POL@IsoRightSign=\z@ + \let\POL@IsoLeft@Int\POL@IsoRight@Int % root at 7 + \def\POL@IsoLeftSign{0}% + \else + \let\POL@@IsoRight@Int\POL@IsoRight@Int % 7 + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% 6 + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int % 6 + \let\POL@IsoRight@Int\POL@@IsoRight@Int % 7 + \let\POL@IsoRightSign\POL@@IsoRightSign + \fi % else 5, 6 with possible root at 6 + \fi\fi + \else + \ifnum\POL@IsoRightSign=\z@ + \let\POL@IsoLeft@Int\POL@IsoRight@Int % root at 5 + \def\POL@IsoLeftSign{0}% + \else + \let\POL@@IsoRight@Int\POL@IsoRight@Int % 5 + \edef\POL@IsoRight@Int{\xintiiAdd{2}{\POL@IsoLeft@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int % 3 + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% 4 + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int % 4 + \let\POL@IsoRight@Int\POL@@IsoRight@Int % 5 + \let\POL@IsoRightSign\POL@@IsoRightSign + \fi % else 3, 4 with possible root at 4 + \else + \ifnum\POL@IsoRightSign=\z@ + \let\POL@IsoLeft@Int\POL@IsoRight@Int % root at 3 + \def\POL@IsoLeftSign{0}% + \else + \let\POL@@IsoRight@Int\POL@IsoRight@Int % 3 + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% 2 + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\POL@IsoLeftSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int % 2 + \let\POL@IsoRight@Int\POL@@IsoRight@Int % 3 + \let\POL@IsoRightSign\POL@@IsoRightSign + \fi % else 1, 2 with possible root at 2 + \fi\fi + \fi\fi +}% +\def\POL@refine@storeleftandright{% + \expandafter + \xdef\csname POL_ZeroInt\POL@sturmname + L\POL@isolz@IntervalIndex\endcsname + {\PolDecToString{\POL@IsoLeft@rawout}}% + \expandafter + \xdef\csname POL_ZeroInt\POL@sturmname + R\POL@isolz@IntervalIndex\endcsname + {\PolDecToString{\POL@IsoRight@rawout}}% + \begingroup\globaldefs\@ne + \xintdefvar\POL@sturmname + L_\POL@isolz@IntervalIndex:=qfrac(\POL@IsoLeft@rawout);% + \xintdefvar\POL@sturmname + R_\POL@isolz@IntervalIndex:=qfrac(\POL@IsoRight@rawout);% + \endgroup +}% + + +%% \PolRefineInterval +\def\POL@xintexprGetVar#1{\expandafter\expandafter\expandafter + \XINT_expr_unlock\csname XINT_expr_var_#1\endcsname}% +\def\POL@set@IsoLeft@rawin{% + \edef\POL@IsoLeft@rawin + {\POL@xintexprGetVar{\POL@sturmname L_\POL@isolz@IntervalIndex}}% +}% +\def\POL@set@IsoRight@rawin{% + \edef\POL@IsoRight@rawin + {\POL@xintexprGetVar{\POL@sturmname R_\POL@isolz@IntervalIndex}}% +}% +\def\POL@set@IsoLeft@Int #1/1[#2]{% + \edef\POL@IsoLeft@Int{\xintDSH{\POL@isolz@E-#2}{#1}}% +}% +\newcommand\PolRefineInterval{\@ifstar\POL@srefine@start\POL@refine@start}% +\newcommand\POL@refine@start[3][1]{% + \edef\POL@isolz@IntervalIndex{\the\numexpr#3}% + \edef\POL@sturmname{#2}% + \expandafter\POL@refine@sharedbody\expandafter + {\expandafter\POL@refine@loop\expandafter{\the\numexpr#1}}% +}% +\def\POL@srefine@start#1#2{% + \edef\POL@isolz@IntervalIndex{\the\numexpr#2}% + \edef\POL@sturmname{#1}% + \POL@refine@sharedbody + {\let\POL@refine@left@next\POL@refine@main % we want to recurse if needed + \let\POL@refine@right@next\POL@refine@main % we want to recurse if needed + \POL@refine@main}% +}% +\def\POL@refine@sharedbody#1{% + \POL@set@IsoLeft@rawin + \edef\POL@IsoLeftSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoLeft@rawin}}}% + \ifnum\POL@IsoLeftSign=\z@ + % do nothing if that interval was already a singleton + \else + % else both end-points are not roots and there is a single one in-between + \POL@set@IsoRight@rawin + \edef\POL@IsoRightSign{\the\numexpr-\POL@IsoLeftSign}% + \edef\POL@isolz@E{\expandafter\POL@refine@getE + % je pense que le xintrez ici est superflu + \romannumeral0\xintrez{\xintSub{\POL@IsoRight@rawin}{\POL@IsoLeft@rawin}}}% + \expandafter\POL@set@IsoLeft@Int\POL@IsoLeft@rawin + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% + #1% + \POL@refine@storeleftandright + \fi +}% +\def\POL@refine@loop#1{% + \let\POL@refine@left@next \@empty % no recursion at end sub-intervals + \let\POL@refine@right@next\@empty + \xintiloop[1+1] + \POL@refine@main + \ifnum\POL@IsoRightSign=\z@ + \expandafter\xintbreakiloop + \fi + \ifnum\xintiloopindex<#1 + \repeat +}% +\def\POL@refine@main{% + \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}% + \edef\POL@IsoLeft@Int{\xintDSL{\POL@IsoLeft@Int}}% + \edef\POL@IsoRight@Int{\xintDSL{\POL@IsoRight@Int}}% + \let\POL@@IsoRight@Int\POL@IsoRight@Int + \let\POL@@IsoRightSign\POL@IsoRightSign + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\z@ + \let\POL@IsoLeft@Int\POL@IsoRight@Int % root at 1 + \def\POL@IsoLeftSign{0}% + \let\POL@next\@empty + \else + \ifnum\POL@IsoRightSign=\POL@@IsoRightSign\space + \let\POL@next\POL@refine@left@next % may be \@empty or \POL@refine@main for recursion + \let\POL@refine@right@next\@empty + \else + \let\POL@IsoLeft@Int\POL@IsoRight@Int + \edef\POL@IsoRight@Int{\xintDec{\POL@@IsoRight@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + \ifnum\POL@IsoRightSign=\z@ + \let\POL@IsoLeft@Int\POL@IsoRight@Int % root at 9 + \def\POL@IsoLeftSign{0}% + \let\POL@next\@empty + \else + \ifnum\POL@IsoRightSign=\POL@@IsoRightSign\space + \let\POL@next\POL@refine@doonce + \else + \let\POL@IsoLeft@Int\POL@IsoRight@Int + \let\POL@IsoRight@Int\POL@@IsoRight@Int + \let\POL@IsoRightSign\POL@@IsoRightSign + \let\POL@next\POL@refine@right@next + \let\POL@refine@left@next\@empty + \fi + \fi + \fi\fi + \POL@next +}% +% lacking pre-defined xintfrac macro here (such as an \xintRawExponent) +\def\POL@refine@getE#1[#2]{#2}% \xintREZ already applied, for safety + + +\newcommand\PolIntervalWidth[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_ZeroInt#1R}{#2}}% + {\@nameuse{POL_ZeroInt#1L}{#2}}} +}% + + +\newcommand\PolEnsureIntervalLengths[2]{% #1 = Sturm chain name, + % localize roots in intervals of length at most 10^{#2} + \POL@count\z@ + % \POL@count used by \POL@sturmchain@getSV@at but latter not used + \edef\POL@sturmname{#1}% + \edef\POL@ensure@targetE{\the\numexpr#2}% + \edef\POL@nbofroots{\csname POL_ZeroInt\POL@sturmname L\endcsname 0}% + \xintloop + \advance\POL@count\@ne + \edef\POL@isolz@IntervalIndex{\the\POL@count}% + \POL@ensure@one + \ifnum\POL@nbofroots>\POL@count + \repeat +}% +\newcommand\PolEnsureIntervalLength[3]{% #1 = Sturm chain name, + % #2 = index of interval + % localize roots in intervals of length at most 10^{#3} + \edef\POL@sturmname{#1}% + \edef\POL@ensure@targetE{\the\numexpr#3}% + \edef\POL@isolz@IntervalIndex{\the\numexpr#2}% + \POL@ensure@one +}% +\def\POL@ensure@one{% + \POL@set@IsoLeft@rawin + \POL@set@IsoRight@rawin + \edef\POL@ensure@delta{\xintREZ{\xintSub{\POL@IsoRight@rawin}{\POL@IsoLeft@rawin}}}% + \xintiiifZero{\POL@ensure@delta} + {} + {\edef\POL@isolz@E{\expandafter\POL@refine@getE\POL@ensure@delta}% + \expandafter\POL@set@IsoLeft@Int\POL@IsoLeft@rawin + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% + \ifnum\POL@isolz@E>\POL@ensure@targetE\space + \edef\POL@IsoLeftSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoLeft@raw}}}% + % at start left and right are not roots, and values of opposite signs + % \edef\POL@IsoRightSign{\the\numexpr-\POL@IsoLeftSign}% + \xintloop + \POL@ensure@Eloopbody % decreases E by one at each iteration + % if separation level is still too coarse we recurse at deeper level + \ifnum\POL@isolz@E>\POL@ensure@targetE\space + \repeat + % will check if right is at a zero, needs \POL@IsoRightSign set up + \POL@refine@storeleftandright + \fi + }% +}% +\def\POL@ensure@Eloopbody {% + \edef\POL@isolz@E{\the\numexpr\POL@isolz@E-\@ne}% + \edef\POL@IsoLeft@Int{\xintDSL{\POL@IsoLeft@Int}}% + % this will loop at most ten times + \xintloop + \edef\POL@IsoRight@Int{\xintInc{\POL@IsoLeft@Int}}% + \edef\POL@IsoRightSign + {\xintiiSgn{\PolEval{\POL@sturmname _0}\At{\POL@IsoRight@raw}}}% + % if we have found a zero at right boundary the \ifnum test will fail + % and we exit the loop + % else we exit the loop if sign at right boundary is opposite of + % sign at left boundary (the latter is +1 or -1, never 0) + % this is a bit wasteful if we go ten times to the right, because + % we know that there the sign will be opposite, evaluation was superfluous + \ifnum\POL@IsoLeftSign=\POL@IsoRightSign\space + \let\POL@IsoLeft@Int\POL@IsoRight@Int + \repeat + % check for case when we exited the inner loop because we actually + % found a zero, then we force exit from the main (E decreasing) loop + \ifnum\POL@IsoRightSign=\z@ + \expandafter\xintbreakloop + \fi +}% + + +\catcode`_ 8 +\newcommand\PolPrintIntervals[2][Z]{% + \POL@count \@nameuse{POL_ZeroInt#2L}{0} + \ifnum\POL@count=\z@ +% No real roots.\par + \else +% There are \the\POL@count\space distinct real roots:\par + \[\count@\POL@count + \global\POL@count\@ne + \begin{array}{rcccl} + \xintloop + \POL@SturmIfZeroExactlyKnown{#2}\POL@count + {% exact root + && + #1_{\the\POL@count}&=& + \POL@printintervals@prepare{#2R}% + \PolPrintIntervalsPrintExactZero + }% + {% interval with root in its strict interior + \POL@printintervals@prepare{#2L}% + \PolPrintIntervalsPrintLeftEndPoint&<& + #1_{\the\POL@count}&<& + \POL@printintervals@prepare{#2R}% + \PolPrintIntervalsPrintRightEndPoint + }% + \global\advance\POL@count\@ne + \unless\ifnum\POL@count>\count@ + \\% + \repeat + \end{array}\] + \fi +}% +\catcode`_ 11 +\newcommand\PolPrintIntervalsPrintExactZero {\PolPrintIntervalsTheEndPoint}% +\newcommand\PolPrintIntervalsPrintLeftEndPoint {\PolPrintIntervalsTheEndPoint}% +\newcommand\PolPrintIntervalsPrintRightEndPoint{\PolPrintIntervalsTheEndPoint}% +\def\POL@printintervals@prepare#1{% + \edef\PolPrintIntervalsTheIndex{\the\POL@count}% + \edef\PolPrintIntervalsTheEndPoint{\@nameuse{POL_ZeroInt#1}\POL@count}% + \xintiiifSgn{\POL@xintexprGetVar{#1_\PolPrintIntervalsTheIndex}} + {\let\PolIfEndPointIsPositive\xint_secondoftwo + \let\PolIfEndPointIsNegative\xint_firstoftwo + \let\PolIfEndPointIsZero\xint_secondoftwo} + {\let\PolIfEndPointIsPositive\xint_secondoftwo + \let\PolIfEndPointIsNegative\xint_secondoftwo + \let\PolIfEndPointIsZero\xint_firstoftwo} + {\let\PolIfEndPointIsPositive\xint_firstoftwo + \let\PolIfEndPointIsNegative\xint_secondoftwo + \let\PolIfEndPointIsZero\xint_secondoftwo}% +}% +\newcommand\POL@SturmIfZeroExactlyKnown[2]{% faster than public one, + % but does not check if #2 is in range + \romannumeral0\xintifeq{\POL@xintexprGetVar{#1L_\the\numexpr#2\relax}}% + {\POL@xintexprGetVar{#1R_\the\numexpr#2\relax}}% +}% + + +\newcommand\PolSturmIfZeroExactlyKnown[2]{% + \romannumeral0\xintifeq{\PolSturmIsolatedZeroLeft{#1}{#2}}% + {\PolSturmIsolatedZeroRight{#1}{#2}}% +}% +\newcommand\PolSturmIsolatedZeroLeft[2]{% + \romannumeral`^^@\csname POL_ZeroInt#1L\endcsname{#2}}% +\newcommand\PolSturmIsolatedZeroRight[2]{% + \romannumeral`^^@\csname POL_ZeroInt#1R\endcsname{#2}}% +\newcommand\PolSturmNbOfIsolatedZeros[1]{% + \romannumeral`^^@\csname POL_ZeroInt#1L0\endcsname +}% + +%% \PolDecToString (should become an xintfrac macro at some point) +\newcommand\PolDecToString[1]{\romannumeral0\expandafter + \POL@dectostring\romannumeral0\xintraw{#1}}% +\def\POL@dectostring #1/#2[#3]{\xintiiifZero {#1}% + \POL@dectostring@z + {\if1\XINT_isOne{#2}\expandafter\POL@dectostring@a + \else\expandafter\POL@dectostring@b + \fi}% + #1/#2[#3]% +}% +\def\POL@dectostring@z#1[#2]{ 0}% +\def\POL@dectostring@a#1/#2[#3]{% + \ifnum#3<\z@\xint_dothis{\xinttrunc{-#3}{#1[#3]}}\fi + \xint_orthat{\xintiie{#1}{#3}}% +}% +\def\POL@dectostring@b#1/#2[#3]{% just to handle this somehow + \ifnum#3<\z@\xint_dothis{\xinttrunc{-#3}{#1[#3]}/#2}\fi + \xint_orthat{\xintiie{#1}{#3}/#2}% +}% + + +\newcommand\PolMakeMonic[1]{% + \edef\POL@leadingcoeff{\PolLeadingCoeff{#1}}% + \edef\POL@leadingcoeff@inverse{\xintDiv{1/1[0]}{\POL@leadingcoeff}}% + \PolMapCoeffs{\xintMul{\POL@leadingcoeff@inverse}}{#1}% }% @@ -635,7 +1490,7 @@ \def\POL@add@b{\POL@get\POL@B\POL@add@c}% \def\POL@add@c{% \global\POL@poltrue - \POL@ifZero\POL@A + \POL@ifZero\POL@A {\let\POL@result\POL@B}% {\POL@ifZero\POL@B {\let\POL@result\POL@A}% @@ -678,7 +1533,7 @@ \def\POL@mul@b{\POL@get\POL@B\POL@mul@c}% \def\POL@mul@c{% \global\POL@poltrue - \POL@ifZero\POL@A + \POL@ifZero\POL@A {\def\POL@result{-1.\empty{0/1[0]}}}% {\POL@ifZero\POL@B {\def\POL@result{-1.\empty{0/1[0]}}}% @@ -731,7 +1586,7 @@ \ifnum\count\tw@<\count@ \repeat \expandafter\let\csname POL@arrayC\the\count@\endcsname\POL@tmp -}% +}% \def\POL@@mul@phaseIIloopbody{% \advance\count@\@ne \def\POL@tmp{0[0]}% @@ -748,7 +1603,7 @@ {\@nameuse{POL@arrayB\the\numexpr\count@+\@ne-\count\tw@}}% }% }% - \repeat + \repeat \expandafter\let\csname POL@arrayC\the\count@\endcsname\POL@tmp }% \def\POL@@mul@phaseIIIloopbody{% @@ -766,7 +1621,7 @@ }% }% \ifnum\@nameuse{POL@arrayA0}>\count\tw@ - \repeat + \repeat \expandafter\let\csname POL@arrayC\the\count@\endcsname\POL@tmp }% @@ -848,7 +1703,7 @@ \ifodd\POL@pow@exp\space \expandafter\POL@@pow@odd \the\numexpr(\POL@pow@exp+\@ne)/\tw@-\@ne\expandafter.% - \else + \else \expandafter\POL@@pow@even \the\numexpr(\POL@pow@exp+\@ne)/\tw@-\@ne\expandafter.% \fi @@ -915,7 +1770,7 @@ \xintScalarDiv{\csname POL@arrayR\the\count@\endcsname}% {\POL@B@leading}}% \expandafter\let\csname POL@arrayQ\the\count\tw@\endcsname - \POL@@div@ratio + \POL@@div@ratio \advance\count@\m@ne \advance\count\tw@\m@ne \count4 \count@ @@ -956,17 +1811,41 @@ %% EXPANDABLE MACROS -\newcommand\PolEval{}% -\def\PolEval#1#2\At#3{\romannumeral`^^@\xinttheexpr #1(#3)\relax}% +\def\Pol@Eval@fork#1\At#2#3\krof{#2}% +\newcommand\PolEval[3]{\romannumeral`^^@\Pol@Eval@fork + #2\PolEvalAt + \At\PolEvalAtExpr\krof {#1}{#3}% +}% +\newcommand\PolEvalAt[2] + {\xintpraw{\csname XINT_expr_userfunc_#1\endcsname{#2},}}% +\newcommand\PolEvalAtExpr[2]{\xinttheexpr #1(#2)\relax}% +% +\newcommand\PolEvalReduced[3]{\romannumeral`^^@\Pol@Eval@fork + #2\PolEvalReducedAt + \At\PolEvalReducedAtExpr\krof {#1}{#3}% +}% +\newcommand\PolEvalReducedAt[2]{% + \xintpraw % in order not to print denominator if the latter equals 1 + {\xintIrr{\csname XINT_expr_userfunc_#1\endcsname{#2},}[0]}% +}% +\newcommand\PolEvalReducedAtExpr[2]{% + \xintpraw + {\xintIrr{\romannumeral`^^@\xintthebareeval#1(#2)\relax}[0]}% +}% % -\newcommand\PolEvalReduced{}% -\def\PolEvalReduced#1#2\At#3{% - \romannumeral0\xintpraw % only serves to not print denominator if = 1 - {\xintIrr{\romannumeral`^^@\xintthebareeval#1(#3)\relax}[0]}% +\newcommand\PolFloatEval[3]{\romannumeral`^^@\Pol@Eval@fork + #2\PolFloatEvalAt + \At\PolFloatEvalAtExpr\krof {#1}{#3}% }% +\newcommand\PolFloatEvalAt[2] + {\xintpfloat{\csname XINT_flexpr_userfunc_#1\endcsname{#2},}}% +\newcommand\PolFloatEvalAtExpr[2]{\xintthefloatexpr #1(#2)\relax}% % -\newcommand\PolFloatEval{}% -\def\PolFloatEval#1#2\At#3{\romannumeral`^^@\xintthefloatexpr #1(#3)\relax}% +\newcommand\PolLeadingCoeff[1]{% + \romannumeral`^^@\expandafter\expandafter\expandafter\xintlastitem + \expandafter\expandafter\expandafter + {\csname POLuserpol@#1\endcsname}% +}% % \newcommand\PolNthCoeff[2]{\romannumeral`^^@% \expandafter\POL@nthcoeff @@ -990,11 +1869,10 @@ \newcommand\PolToCSV[1]{\romannumeral0\xintlistwithsep{, }{\PolToList{#1}}}% -\newif\ifpoltoexprall \newcommand\PolToExprCmd[1]{\xintPRaw{\xintRawWithZeros{#1}}}% \newcommand\PolToFloatExprCmd[1]{\xintFloat{#1}}% \let\PolToExprTermPrefix\PolTypesetCmdPrefix -\newcommand\PolToExprOneTerm[2]{% +\newcommand\PolToExprOneTermStyleA[2]{% \ifnum#2=\z@ \PolToExprCmd{#1}% \else @@ -1007,6 +1885,7 @@ \else\PolToExprVar^\xintiiAbs{#2}% \fi }% +\let\PolToExprOneTerm\PolToExprOneTermStyleA \newcommand\PolToExprOneTermStyleB[2]{% \ifnum#2=\z@ \xintNumerator{#1}% |