% Arbitrary precision numbers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2014, 2015, 2016, 2018 Petr Olsak
% See the documentation apnum.pdf or apnum.d for more information
\def\apVERSION{1.7 <Apr 2018>}
\message{The Arbitrary Precision Numbers, \apVERSION}
%%%%%%%%%%%% Internal registers, sec. 2.1 in apnum.pdf
\newcount\apnumA \newcount\apnumB \newcount\apnumC \newcount\apnumD
\newcount\apnumE \newcount\apnumF \newcount\apnumG \newcount\apnumH
\newcount\apnumO \newcount\apnumP \newcount\apnumL
\newcount\apnumX \newcount\apnumY \newcount\apnumZ
\newcount\apSIGNa \newcount\apSIGNb \newcount\apEa \newcount\apEb
\newif\ifapX
\newcount\apSIGN
\newcount\apE
\newcount\apTOT \apTOT=0
\newcount\apFRAC \apFRAC=20
\newcount\apEX \apEX=10
\apnumZ=\catcode`\@ \catcode`\@=12
%%%%%%%%%%%% Evaluation of the expression, sec. 2.2 in apnum.pdf
\def\evaldef{\relax \apEVALa{\apEadd\OUT}}
\def\evalmdef{\relax \apEVALa{}}
\def\apEVALa#1#2#3{\begingroup \apnumA=0 \apnumE=1 \apEVALb#3\limits \tmpb \apEND #1\let#2=\OUT}
\def\apEVALb{\def\tmpa{}\apEVALc}
\def\apEVALc#1{%
\ifx+#1\apEVALd \apEVALc \fi
\ifx-#1\edef\tmpa{\tmpa-}\apEVALd\apEVALc \fi
\ifx(#1\apEVALd \apEVALe \fi
\ifx\the#1\apEVALd \apEVALf\the\fi
\ifx\number#1\apEVALd \apEVALf\number\fi
\apTESTdigit#1\iftrue
\ifx E#1\let\tmpb=\tmpa \expandafter\apEVALd\expandafter\apEVALk
\else \edef\tmpb{\tmpa#1}\expandafter\apEVALd\expandafter\apEVALn\fi\fi
\edef\tmpb{\tmpa\noexpand#1}\expandafter
\futurelet\expandafter\apNext\expandafter\apEVALg\romannumeral-`\.%
}
\def\apEVALd#1\fi#2-`\.{\fi#1}
\def\apEVALe{%
\ifx\tmpa\empty \else \ifnum\tmpa1<0 \def\tmpb{-1}\apEVALp \apMUL 4\fi\fi
\advance\apnumA by4
\apEVALb
}
\def\apEVALf#1#2{\expandafter\def\expandafter\tmpb\expandafter{\tmpa#1#2}\apEVALo}
\def\apEVALg{\ifx\apNext \bgroup \expandafter\apEVALh \else \expandafter\apEVALo \fi}
\def\apEVALh#1{\expandafter\def\expandafter\tmpb\expandafter{\tmpb{#1}}\expandafter
\futurelet\expandafter\apNext\expandafter\apEVALg\romannumeral-`\.}
\def\apEVALk{\afterassignment\apEVALm\apE=}
\def\apEVALm{\edef\tmpb{\tmpb E\the\apE}\apEVALo}
\def\apEVALn#1{\apTESTdigit#1%
\iftrue \ifx E#1\afterassignment\apEVALm\expandafter\expandafter\expandafter\apE
\else\edef\tmpb{\tmpb#1}\expandafter\expandafter\expandafter\apEVALn\fi
\else \expandafter\apEVALo\expandafter#1\fi
}
\def\apEVALo#1{\let\apNext=\apEVALb
\ifx+#1\apEVALp \apPLUS 1\fi
\ifx-#1\apEVALp \apMINUS 1\fi
\ifx*#1\apEVALp \apMUL 2\fi
\ifx/#1\apEVALp \apDIV 2\fi
\ifx^#1\apEVALp \apPOWx 3\fi
\ifx)#1\advance\apnumA by-4 \let\apNext=\apEVALo \let\tmpa=\relax
\ifnum\apnumA<0 \apEVALerror{many brackets ")"}\fi
\fi
\ifx\limits#1%
\ifnum\apnumA>0 \apEVALerror{missing bracket ")"}\let\tmpa=\relax
\else \apEVALp\END 0\let\apNext=\relax \fi
\fi
\ifx\tmpa\relax \else \apEVALerror{unknown operator "\string#1"}\fi
\apnumE=0 \apNext
}
\def\apEVALp#1#2{%
\apnumB=#2 \advance\apnumB by\apnumA
\toks0=\expandafter{\expandafter{\tmpb}{#1}}%
\expandafter\apEVALpush\the\toks0\expandafter{\the\apnumB}% {value}{op}{priority}
\let\tmpa=\relax
}
\def\apEVALstack{{}{}{0}.}
\def\apEVALpush#1#2#3{% value, operator, priority
\toks0={{#1}{#2}{#3}}%
\expandafter\def\expandafter\apEVALstack\expandafter{\the\toks0\apEVALstack}%
\expandafter\apEVALdo\apEVALstack@%
}
\def\apEVALdo#1#2#3#4#5#6#7@{%
\apnumB=#3 \ifx#2\apPOWx \advance\apnumB by1 \fi
\ifnum\apnumB>#6\else
\ifnum#6=0 \def\tmpb{#1}%\toks0={#1}\message{RESULT: \the\toks0}
\ifnum\apnumE=1 \def\tmpb{\apPPn{#1}}\fi
\else \def\apEVALstack{#7}\apEVALpush{#5{#4}{#1}}{#2}{#3}%
\fi\fi
}
\def\apEVALerror#1{\message{\noexpand\evaldef ERROR: #1.}%
\def\OUT{0}\apE=0\apSIGN=0\def\apNext##1\apEND{\apEND}%
}
\def\apTESTdigit#1#2{%
\ifx E#1\apXtrue \else
\ifcat.\noexpand#1%
\ifx.#1\apXtrue \else
\ifnum`#1<`0 \apXfalse\else
\ifnum`#1>`9 \apXfalse\else \apXtrue\fi
\fi\fi
\else \apXfalse
\fi\fi
\ifapX
}
%%%%%%%%%%%% Preparation of the parameter, sec. 2.3 in apnum.pdf
\def\apPPa#1#2{\expandafter\apPPb#2@#1}
\def\apPPb{\def\tmpc{}\apSIGN=1 \apE=0 \expandafter\expandafter\expandafter\apPPc}
\def\apPPc#1{%
\ifx+#1\apPPd \fi
\ifx-#1\apSIGN=-\apSIGN \apPPd \fi
\ifx\relax#1\apPPe \fi
\apPPg#1%
}
\def\apPPd#1\apPPg#2{\fi\expandafter\expandafter\expandafter\apPPc}
\def\apPPe#1\apPPg#2#3@{\fi
\begingroup\apE=0 #3% execution of the parameter in the group
\edef\tmpb{\apE=\the\apE\relax\noexpand\apPPf\OUT@}\expandafter\endgroup\tmpb
}
\def\apPPf#1{\ifx-#1\apSIGN=-\apSIGN \expandafter\apPPg\else\expandafter\apPPg\expandafter#1\fi}
\def\apPPg#1{%
\ifx.#1\def\tmpc{.}\apPPh\fi
\ifx\tmpc\empty\else\edef\tmpc{\tmpc#1}\fi
\ifx0#1\apPPh\fi
\ifx\tmpc\empty\edef\tmpc{#1}\fi
\ifx@#1\def\tmpc{@}\apSIGN=0 \fi
\expandafter\apPPi\tmpc
}
\def\apPPh#1\apPPi\tmpc{\fi\apPPg}
\def\apPPi{\ifnum\apE=0 \expandafter\apPPk \else \expandafter\apPPj \fi}
\def\apPPj#1@#2{\def#2{#1}}
\def\apPPk#1@#2{\ifx@#1@\apSIGN=0 \def#2{0}\else \apPPl#1E@#2\fi}
\def\apPPl#1E#2@#3{%
\ifx@#1@\def#3{1}\else\def#3{#1}\fi
\ifx@#2@\else \afterassignment\apPPm \apE=#2\fi
}
\def\apPPm E{}
\def\apPPn#1{\expandafter\apPPb#1@\OUT
\ifnum\apSIGN=0 \def\OUT{0}\fi
\ifnum\apSIGN<0 \edef\OUT{-\OUT}\fi
}
\def\apPPab#1#2#3{%
\expandafter\apPPb#2@\tmpa \apSIGNa=\apSIGN \apEa=\apE
\expandafter\apPPb#3@\tmpb \apSIGNb=\apSIGN \apEb=\apE
#1%
}
\def\apPPs#1#2#3{\def\tmpc{#3}\expandafter\apPPt\expandafter#1#2.@#2}
\def\apPPt#1#2{%
\ifx-#2\apnumG=-1 \def\apNext{#1}%
\else \ifx0#2\apnumG=0 \def\apNext{\apPPu#1}\else \apnumG=1 \def\apNext{#1#2}\fi\fi
\apNext
}
\def\apPPu#1#2.@#3{\ifx@#2@\apnumG=0 \ifx#1\apROUNDa\def\XOUT{}\fi
\else\def\apNext{\apPPt#1#2.@#3}\expandafter\apNext\fi
}
%%%%%%%%%%%% Addition and Subtraction, sec. 2.4 in apnum.pdf
\def\apPLUS{\relax \apPPab\apPLUSa}
\def\apMINUS#1#2{\relax \apPPab\apPLUSa{#1}{-#2}}
\def\apPLUSa{%
\ifnum\apEa=\apEb \apE=\apEa \else \apPLUSxE \fi
\apDIG\tmpa\relax \apnumA=\apnumD % digits before decimal point
\apDIG\tmpb\relax \apnumB=\apnumD
\apIVmod \apnumA \apnumE \advance\apnumA by-\apnumE % digits in the first Digit
\apIVmod \apnumB \apnumF \advance\apnumB by-\apnumF
\apnumC=\apnumB \advance\apnumC by-\apnumA % difference between Digits
\ifnum\apSIGNa<0 \def\apPLUSxA{-}\else \def\apPLUSxA{}\fi
\ifnum\apSIGNb<0 \def\apPLUSxB{-}\else \def\apPLUSxB{}\fi
\apSIGN=0 % \apSIGN=0 means that we are doing subtraction
\ifx\apPLUSxA\empty \ifx\apPLUSxB\empty \apSIGN=1 \fi\fi
\if\apPLUSxA-\relax \if\apPLUSxB-\relax \apSIGN=-1 \def\apPLUSxA{}\def\apPLUSxB{}\fi\fi
\ifnum\apnumC>0 \apPLUSg \apPLUSb \tmpb\apnumF \tmpa\apnumE \apnumB % first pass
\else \apnumC=-\apnumC \apPLUSb \tmpa\apnumE \tmpb\apnumF \apnumA
\fi
\ifnum\apnumG=0 \def\OUT{0}\apSIGN=0 \apE=0 \else
\ifnum\apSIGN=0 \apSIGN=\apnumG \let\apNext=\apPLUSm \else \let\apNext=\apPLUSp \fi
\apnumX=0 \edef\OUT{\expandafter}\expandafter \apNext \OUT@% second pass
\ifnum\apnumD<1 % result in the form .000123
\apnumZ=-\apnumD
\def\tmpa{.}%
\ifnum\apnumZ>0 \apADDzeros\tmpa \fi % adding dot and left zeros
\edef\OUT{\ifnum\apSIGN<0-\fi\tmpa\OUT}%
\else
\edef\OUT{\expandafter}\expandafter\apPLUSy \OUT@% removing left zeros
\fi\fi
}
\def\apPLUSb#1#2#3#4#5{%
\edef\tmpd{\ifcase#4\or{}{}{}\or{}{}\or{}\fi#3}%
\edef\tmpc{\ifcase#2\or{}{}{}\or{}{}\or{}\fi}%
\let\apNext=\apPLUSc \apnumD=#5\advance\apnumD by4 \apnumG=0 \apnumZ=0 \def\OUT{}%
\expandafter\expandafter\expandafter\apPLUSc\expandafter\tmpc#1\apNL\apNL\apNL\apNL@%
}
\def\apPLUSc#1#2#3#4{\apnumY=\apPLUSxA#1#2#3#4\relax
\ifx\apNL#4\let\apNext=\apPLUSd\fi
\ifx\apNL#1\relax \ifx\tmpd\empty \expandafter\expandafter\expandafter\apPLUSf \fi\fi
\apPLUSe
}
\def\apPLUSd{\apnumY=0 \ifx\tmpd\empty \expandafter\apPLUSf \else\expandafter \apPLUSe\fi}
\def\apPLUSe{%
\ifnum\apnumC>0 \advance\apnumC by-4
\else \apIVread\tmpd \advance\apnumY by\apPLUSxB\apnumX \fi
\ifnum\apnumZ=0 \apPLUSh \fi
\edef\OUT{{\the\apnumY}\OUT}%
\advance\apnumD by-4
\apNext
}
\def\apPLUSf#1@{}
\def\apPLUSg{\let\tmpc=\apPLUSxA \let\apPLUSxA=\apPLUSxB \let\apPLUSxB=\tmpc}
\def\apPLUSh{\apnumZ=\apnumY
\ifnum\apnumY=0 \else \ifnum\apnumY<0 \apnumG=-1 \apnumY=-\apnumY \apPLUSg \else\apnumG=1 \fi\fi
}
\def\apPLUSm#1{%
\ifx@#1\else
\apnumA=#1 \advance\apnumA by-\apnumX
\ifnum\apnumA<0 \advance\apnumA by\apIVbase \apnumX=1 \else \apnumX=0 \fi
\apPLUSw
\expandafter\apPLUSm
\fi
}
\def\apPLUSp#1{%
\ifx@#1\ifnum\apnumX>0 \apnumA=1 \apPLUSw \fi % .5+.5=.1 bug fixed
\else
\apnumA=\apnumX \advance\apnumA by#1
\ifnum\apnumA<\apIVbase \apnumX=0 \else \apnumX=1 \advance\apnumA by-\apIVbase \fi
\apPLUSw
\expandafter\apPLUSp
\fi
}
\def\apPLUSw{%
\ifnum\apnumD=0 \ifx\OUT\empty \def\OUT{\empty}\else \edef\OUT{.\OUT}\fi \fi
\advance\apnumD by4
\ifx\OUT\empty \edef\tmpa{\apIVwrite\apnumA}\edef\OUT{\apREMzerosR\tmpa}%
\else \edef\OUT{\apIVwrite\apnumA\OUT}\fi
}
\def\apPLUSy#1{\ifx0#1\expandafter\apPLUSy\else \expandafter\apPLUSz\expandafter#1\fi}
\def\apPLUSz#1@{\edef\OUT{\ifnum\apSIGN<0-\fi#1}}
\def\apPLUSxE{%
\apnumE=\apEa \advance\apnumE by-\apEb
\ifnum\apEa>\apEb \apPPs\apROLLa\tmpb{-\apnumE}\apE=\apEa
\else \apPPs\apROLLa\tmpa{\apnumE}\apE=\apEb \fi
}
%%%%%%%%%%%% Multiplication, sec. 2.5 in apnum.pdf
\def\apMUL{\relax \apPPab\apMULa}
\def\apMULa{%
\apE=\apEa \advance\apE by\apEb
\apSIGN=\apSIGNa \multiply\apSIGN by\apSIGNb
\ifnum\apSIGN=0 \def\OUT{0}\apE=0 \else
\apDIG\tmpa\apnumA \apnumX=\apnumA \advance\apnumA by\apnumD
\apDIG\tmpb\apnumB \advance\apnumX by\apnumB \advance\apnumB by\apnumD
\apnumD=\apnumX % \apnumD = the number of digits after decimal point in the result
\apIVmod \apnumA \apnumF % \apnumF = digits in the first Digit of \tmpa
\edef\tmpc{\ifcase\apnumF\or{}{}{}\or{}{}\or{}\fi}\def\OUT{}%
\expandafter\expandafter\expandafter \apMULb \expandafter \tmpc \tmpa @@@@%
\edef\OUT{*.\OUT}%
\apIVmod \apnumB \apnumF % \apnumF = digits in the first Digit of \tmpb
\edef\tmpc{\ifcase\apnumF\or{}{}{}\or{}{}\or{}\fi}\def\tmpa{}%
\expandafter\expandafter\expandafter \apMULc \expandafter \tmpc \tmpb @@@@%
\expandafter\apMULd \tmpa@%
\expandafter\apMULg \OUT
\edef\tmpa{\ifnum\apSIGN<0-\fi}%
\ifnum\apnumD>0 \apnumZ=\apnumD \edef\tmpa{\tmpa.}\apADDzeros\tmpa \fi
\ifx\tmpa\empty \else \edef\OUT{\tmpa\OUT}\fi
\fi
}
\def\apMULb#1#2#3#4{\ifx@#4\else
\ifx\OUT\empty \edef\OUT{{#1#2#3#4}*}\else\edef\OUT{{#1#2#3#4}0\OUT}\fi
\expandafter\apMULb\fi
}
\def\apMULc#1#2#3#4{\ifx@#4\else \edef\tmpa{{#1#2#3#4}\tmpa}\expandafter\apMULc\fi}
\def\apMULd#1{\ifx@#1\else
\apnumA=#1 \expandafter\apMULe \OUT
\expandafter\apMULd
\fi
}
\def\apMULe#1*#2{\apnumX=0 \def\OUT{#1{#2}*}\def\apOUTl{}\apnumO=1 \apnumL=0 \apMULf}
\def\apMULf#1#2{%
\advance\apnumO by-1 \ifnum\apnumO=0 \apOUTx \fi
\apnumB=#1 \multiply\apnumB by\apnumA \advance\apnumB by\apnumX
\ifx*#2%
\ifnum\apnumB<\apIVbase
\edef\OUT{\OUT\expandafter\apOUTs\apOUTl.,\ifnum\the\apnumB#1=0 \else{\the\apnumB}{#1}\fi*}%
\else \apIVtrans
\expandafter \edef\csname apOUT:\apOUTn\endcsname
{\csname apOUT:\apOUTn\endcsname{\the\apnumB}{#1}}%
\apMULf0*\fi
\else \advance\apnumB by#2
\ifnum\apnumB<\apIVbase \apnumX=0 \else \apIVtrans \fi
\expandafter
\edef\csname apOUT:\apOUTn\endcsname{\csname apOUT:\apOUTn\endcsname{\the\apnumB}{#1}}%
\expandafter\apMULf \fi
}
\def\apMULg#1{\def\OUT{}\apMULh}
\def\apMULh#1{\ifx*#1\expandafter\apMULi
\else \apnumA=#1 \apMULo4{\apIVwrite\apnumA}%
\expandafter\apMULh
\fi
}
\def\apMULi#1#2#3{\apnumA=#1
\ifx*#3\apMULo{\apNUMdigits\tmpa}{\the\apnumA}\expandafter\apMULj
\else \apMULo4{\apIVwrite\apnumA}\expandafter\apMULi
\fi{#3}%
}
\def\apMULj#1{}
\def\apMULo#1#2{\edef\tmpa{#2}%
\advance\apnumD by-#1
\ifnum\apnumD<1 \ifnum\apnumD>-4 \apMULt\fi\fi
\edef\OUT{\tmpa\OUT}%
}
\def\apMULt{\edef\tmpa{\apIVdot{-\apnumD}\tmpa}\edef\tmpa{\tmpa}}
%%%%%%%%%%%% Division, sec. 2.6 in apnum.pdf
\def\apDIV{\relax \apPPab\apDIVa}
\def\apDIVa{%
\ifnum\apSIGNb=0 \apERR{Dividing by zero}\else
\apSIGN=\apSIGNa \multiply\apSIGN by\apSIGNb
\ifnum\apSIGNa=0 \def\OUT{0}\def\XOUT{0}\apE=0 \apSIGN=0 \else
\apE=\apEa \advance\apE by-\apEb
\apDIG\tmpb\relax \apnumB=\apnumD
\apDIG\tmpa\relax \apnumH=\apnumD
\advance\apnumD by-\apnumB % \apnumD = num. of digits before decimal point in the result
\apDIVcomp\tmpa\tmpb % apXtrue <=> A>=B, i.e 1 digit from A/B
\ifapX \advance\apnumD by1 \advance\apnumH by1 \fi
\apnumC=\apTOT
\ifnum\apTOT<0 \apnumC=-\apnumC
\ifnum\apnumD>\apnumC \apnumC=\apnumD \fi
\fi
\ifnum\apTOT=0 \apnumC=\apFRAC \advance\apnumC by\apnumD
\else \apnumX=\apFRAC \advance\apnumX by\apnumD
\ifnum\apnumC>\apnumX \apnumC=\apnumX \fi
\fi
\ifnum\apnumC>0 % \apnumC = the number of digits in the result
\advance\apnumH by-\apnumC % \apnumH = the position of decimal point in the remainder
\apIVmod \apnumC \apnumF % \apnumF = the number of digits in the first Digit
\apIVread\tmpb \apnumB=\apnumX % \apnumB = partial divisor
\apnumX=\apnumF \ifapX \advance\apnumX by-1 \fi
\apIVreadX\apnumX\tmpa
\apnumA=\apnumX % \apnumA = first Digit of the partial dividend
\apIVread\tmpa % \apnumX = second Digit of the partial dividend
\edef\apDIVxA{\the\apnumA\apIVwrite\apnumX}% first partial dividend
\edef\apDIVxB{\the\apnumB}% partial divisor
\edef\XOUT{{\apDIVxB}{\the\apnumX}@{\the\apnumA}}% the \XOUT is initialized
\edef\OUT{\ifnum\apSIGN<0-\fi}%
\ifnum\apnumD<0 \edef\OUT{\OUT.}\apnumZ=-\apnumD \apADDzeros\OUT \fi
\apnumE=1 \apnumZ=0
\let\apNext=\apDIVg \apNext % <--- the main calculation loop is here
\ifnum\apnumD>0 \apnumZ=\apnumD \apADDzeros\OUT \fi
\ifnum\apnumE=0 \def\XOUT{0}\else % extracting remainder from \XOUT
\edef\XOUT{\expandafter}\expandafter\apDIVv\XOUT
\def\tmpc{\apnumH}\apnumG=\apSIGNa \expandafter\apROLLa\XOUT.@\XOUT
\fi
\else
\def\OUT{0}\def\XOUT{0}\apE=0 \apSIGN=0
\fi\fi\fi
}
\def\apDIVcomp#1#2{%
\expandafter\def\expandafter\tmpc\expandafter{#1\apNL\apNL\apNL\apNL\apNL\apNL\apNL\apNL@}%
\expandafter\def\expandafter\tmpd\expandafter{#2\apNL\apNL\apNL\apNL\apNL\apNL\apNL\apNL@}%
\def\apNext{\expandafter\expandafter\expandafter\apDIVcompA\expandafter\tmpc\tmpd}%
\apXtrue \apNext
}
\def\apDIVcompA#1#2#3#4#5#6#7#8#9@{%
\ifx#8\apNL \def\tmpc{0000000\apNL@}\else\def\tmpc{#9@}\fi
\apnumX=#1#2#3#4#5#6#7#8\relax
\apDIVcompB
}
\def\apDIVcompB#1#2#3#4#5#6#7#8#9@{%
\ifnum\apnumX<#1#2#3#4#5#6#7#8 \let\apNext=\relax \apXfalse \else
\ifnum\apnumX>#1#2#3#4#5#6#7#8 \let\apNext=\relax \apXtrue
\fi\fi
\ifx\apNext\relax\else
\ifx#8\apNL \def\tmpd{0000000\apNL@}\ifx\tmpc\tmpd\let\apNext=\relax\fi
\else\def\tmpd{#9@}\fi
\fi
\apNext
}
\def\apDIVg{%
\ifx\tmpb\empty
\ifx\tmpa\empty \def\apNext{\apDIVi!}\let\apNexti=\apDIVi
\else \def\apNext{\expandafter\apDIVh\tmpa\apNL\apNL\apNL\apNL!}\let\apNexti=\apDIVh
\fi\fi
\ifx\apNext\apDIVg
\apIVread\tmpa \apnumA=\apnumX
\apIVread\tmpb
\edef\XOUT{{\the\apnumX}{\the\apnumA}\XOUT}%
\fi
\apNext
}
\def\apDIVh#1#2#3#4{\apnumZ=#1#2#3#4
\ifx\apNL#4\let\apNexti=\apDIVi\fi
\apDIVi
}
\def\apDIVi{%
\ifnum\apnumE=0 \apnumC=0 \fi
\ifnum\apnumC>0
\expandafter\apDIVp\XOUT
\advance\apnumC by-4
\apnumZ=0
\expandafter\apNexti
\else
\expandafter\apDIVj
\fi
}
\def\apDIVj#1!{}
\def\apDIVp{%
\apnumA=\apDIVxA \divide\apnumA by\apDIVxB
\def\apOUTl{}\apnumO=1 \apnumL=0
\apnumX=0 \apnumB=0 \apnumE=0
\let\apNext=\apDIVq \apNext 0\apnumZ
}
\def\apDIVq#1#2#3{% B A B
\advance\apnumO by-1 \ifnum\apnumO=0 \apOUTx \fi
\apnumY=\apnumB
\apnumB=#1\multiply\apnumB by-\apnumA
\advance\apnumB by#2\advance\apnumB by-\apnumX
\ifnum\apnumB<0 \apnumX=\apnumB \advance\apnumX by1
\divide\apnumX by-\apIVbase \advance\apnumX by1
\advance\apnumB by\the\apnumX 0000
\else \apnumX=0 \fi
\expandafter
\edef\csname apOUT:\apOUTn\endcsname{\csname apOUT:\apOUTn\endcsname{#3}{\the\apnumB}}%
\ifnum\apnumE<\apnumB \apnumE=\apnumB \fi
\ifx@#3\let\apNext=\apDIVr \fi
\apNext{#3}%
}
\def\apDIVr#1#2{%
\ifnum\apnumX=#2 % the calculated Digit is OK, we save it
\edef\XOUT{\expandafter\apOUTs\apOUTl.,}%
\edef\tmpa{\ifnum\apnumF=4 \expandafter\apIVwrite\else \expandafter\the\fi\apnumA}%
\ifnum\apnumD<\apnumF \ifnum\apnumD>-1 \apDIVt \fi\fi %adding dot
\ifx\apNexti\apDIVh \apnumE=1 \fi
\ifnum\apnumE=0 \apDIVu % removing zeros
\advance\apnumD by-\apNUMdigits\tmpa \relax
\else \advance\apnumD by-\apnumF \apnumF=4 \fi
\edef\OUT{\OUT\tmpa}% save the Digit
\edef\apDIVxA{\the\apnumB\apIVwrite\apnumY}% next partial dvividend
\else % we need do correction and run the remainder calculation again
\advance\apnumA by-1 \apnumX=0 \apnumB=0 \apnumE=0
\def\apOUTl{}\apnumO=1 \apnumL=0
\def\apNext{\let\apNext=\apDIVq
\expandafter\apNext\expandafter0\expandafter\apnumZ\XOUT}%
\expandafter\apNext
\fi
}
\def\apDIVt{\edef\tmpa{\apIVdot\apnumD\tmpa}\edef\tmpa{\tmpa}}
\def\apDIVu{\edef\tmpa{\apREMzerosR\tmpa}\edef\tmpa{\apREMdotR\tmpa}}
\def\apDIVv#1#2{\apnumX=#2
\ifx@#1\apDIVw{.\apIVwrite\apnumX}\else\apDIVw{\apIVwrite\apnumX}\expandafter\apDIVv\fi
}
\def\apDIVw#1{%
\ifx\XOUT\empty \ifnum\apnumX=0
\else \edef\tmpa{#1}\edef\XOUT{\apREMzerosR\tmpa\XOUT}%
\fi
\else \edef\XOUT{#1\XOUT}\fi
}
%%%%%%%%%%%% Power to the integer, sec. 2.7 in apnum.pdf
\def\apPOW{\relax \apPPab\apPOWa} \let\apPOWx=\apPOW % for usage as ^ operator
\def\apPOWa{%
\ifnum\apSIGNa=0 \def\OUT{0}\apSIGN=0 \apE=0 \else
\ifnum\apSIGNb=0 \def\OUT{1}\apSIGN=1 \apE=0 \else
\apDIG\tmpb\apnumB
\ifnum\apnumB>0 \apERR{POW: non-integer exponent is not implemented yet}\apPOWe\fi
\ifnum\apEb=0 \else \apERR{POW: the E notation of exponent isn't allowed}\apPOWe\fi
\ifnum\apnumD>8 \apERR{POW: too big exponent.
Do you really need about 10^\the\apnumD\space digits in output?}\apPOWe\fi
\apE=\apEa \multiply\apE by\tmpb\relax
\apSIGN=\apSIGNa
\ifodd\tmpb \else \apSIGN=1 \fi
\apDIG\tmpa\apnumA \apnumC=\apnumA \advance\apnumC by\apnumD
\apnumD=\apnumA \multiply\apnumD by\tmpb
\apIVmod \apnumC \apnumA
\edef\tmpc{\ifcase\apnumA\or{}{}{}\or{}{}\or{}\fi}\def\OUT{}%
\expandafter\expandafter\expandafter \apMULb \expandafter \tmpc \tmpa @@@@%
\edef\OUT{*.\OUT}% \OUT := \tmpa in interleaved format
\def\tmpc{*.1*}%
\apnumE=\tmpb\relax \apPOWb
\expandafter\apPOWg \tmpc % \OUT := \tmpc in human raedable form
\ifnum\apnumD=0 \ifnum \apSIGN<0 \edef\OUT{-\OUT}\fi
\else \def\tmpc{-\apnumD}\apnumG=\apSIGN \expandafter\apROLLa\OUT.@\OUT\fi
\ifnum\apSIGNb<0 \apPPab\apDIVa 1\OUT \fi
\relax
\fi\fi
}
\def\apPOWb{%
\ifodd\apnumE \def\tmpb{}\expandafter\apPOWd\OUT
\let\tmpd=\OUT \let\OUT=\tmpc
\expandafter\apMULd \tmpb@\expandafter\apPOWn\OUT@%
\let\tmpc=\OUT \let\OUT=\tmpd
\fi
\divide\apnumE by2
\ifnum\apnumE>0 \expandafter\apPOWt\OUT \expandafter\apPOWn\OUT@%
\expandafter\apPOWb
\fi
}
\def\apPOWd#1#2{% \apPOWd <spec format> => \tmpb (in simple reverse format)
\ifx*#1\expandafter\apPOWd \else
\edef\tmpb{\tmpb{#1}}%
\ifx*#2\else \expandafter\expandafter\expandafter\apPOWd\fi
\fi
}
\def\apPOWe#1\relax{\fi}
\def\apPOWg#1#2{\def\OUT{}\apPOWh} % conversion to the human readable form
\def\apPOWh#1#2{\apnumA=#1
\ifx*#2\edef\OUT{\the\apnumA\OUT}\else \edef\OUT{\apIVwrite\apnumA\OUT}\expandafter\apPOWh\fi
}
\def\apPOWn#1{\def\OUT{*}\apPOWna}
\def\apPOWna#1{\ifx*#1\expandafter\apPOWnn\else \edef\OUT{\OUT0{#1}}\expandafter\apPOWna\fi}
\def\apPOWnn#1#2{\ifx*#1\edef\OUT{\OUT*}\else\edef\OUT{\OUT0{#1}}\expandafter\apPOWnn\fi}
\def\apPOWt#1#2{\apPOWu} % power to two
\def\apPOWu#1#2{\apnumA=#1
\expandafter\apPOWv\OUT
\ifx*#2\else \expandafter\apPOWu\fi
}
\def\apPOWv#1*#2#3#4{\def\apOUTl{}\apnumO=1 \apnumL=0
\apnumB=\apnumA \multiply\apnumB by\apnumB \multiply\apnumA by2
\ifx*#4\else\advance\apnumB by#4 \fi
\ifx\apnumB<\apIVbase \apnumX=0 \else \apIVtrans \fi
\edef\OUT{#1{#2}{\the\apnumB}*}%
\ifx*#4\apMULf0*\else\expandafter\apMULf\fi
}
%%%%%%%%%%%% ROLL, ROUND and NORM macros, sec. 2.8 in apnum.pdf
\def\apROLL{\apPPs\apROLLa}
\def\apROLLa{\apnumA=\tmpc\relax \ifnum\apnumA<0 \expandafter\apROLLc\else \expandafter\apROLLg\fi}
\def\apROLLc{\edef\tmpc{}\edef\tmpd{\ifnum\apnumG<0-\fi}\apnumB=0 \apROLLd}
\def\apROLLd#1{%
\ifx.#1\expandafter\apROLLe
\else \edef\tmpc{\tmpc#1}%
\advance\apnumB by1
\expandafter\apROLLd
\fi
}
\def\apROLLe#1{\ifx@#1\edef\tmpc{\tmpc.@}\else\edef\tmpc{\tmpc#1}\fi
\advance\apnumB by\apnumA
\ifnum\apnumB<0
\apnumZ=-\apnumB \edef\tmpd{\tmpd.}\apADDzeros\tmpd
\expandafter\expandafter\expandafter\apROLLf\expandafter\tmpc
\else
\apnumA=\apnumB
\expandafter\expandafter\expandafter\apROLLi\expandafter\tmpc
\fi
}
\def\apROLLf#1.@#2{\edef#2{\tmpd#1}}
\def\apROLLg#1{\edef\tmpd{\ifnum\apnumG<0-\fi}\ifx.#1\apnumB=0 \else\apnumB=1 \fi \apROLLh#1}
\def\apROLLh#1{\ifx.#1\expandafter\apROLLi\else \edef\tmpd{\tmpd#1}\expandafter\apROLLh\fi}
\def\apROLLi#1{\ifx.#1\expandafter\apROLLi\else
\ifnum\apnumA>0 \else \apROLLj \apROLLk#1\fi
\ifx@#1\apROLLj \apROLLi0@\fi
\advance\apnumA by-1
\ifx0#1\else \apnumB=1 \fi
\ifnum\apnumB>0 \edef\tmpd{\tmpd#1}\fi
\expandafter\apROLLi\fi
}
\def\apROLLj#1\fi#2\apROLLi\fi{\fi\fi#1}
\def\apROLLk#1{\ifx@#1\expandafter\apROLLo\expandafter@\else
\def\tmpc{}\apnumB=0 \expandafter\apROLLn\expandafter#1\fi
}
\def\apROLLn#1{%
\ifx.#1\ifnum\apnumB>0 \edef\tmpd{\tmpd.\tmpc}\fi \expandafter\apROLLo
\else \edef\tmpc{\tmpc#1}\advance\apnumB by#1 \expandafter\apROLLn
\fi
}
\def\apROLLo@#1{\let#1=\tmpd}
\def\apROUND{\apPPs\apROUNDa}
\def\apROUNDa{\apnumD=\tmpc\relax
\ifnum\apnumD<0 \expandafter\apROUNDe
\else \expandafter\apROUNDb
\fi
}
\def\apROUNDb#1.{\edef\tmpc{#1}\apnumX=0 \def\tmpd{}\let\apNext=\apROUNDc \apNext}
\def\apROUNDc#1{\ifx@#1\def\apNext{\apROUNDd.@}%
\else \advance\apnumD by-1
\ifnum\apnumD<0 \def\apNext{\apROUNDd#1}%
\else \ifx.#1\else \advance\apnumX by#1 \edef\tmpd{\tmpd#1}\fi
\fi
\fi \apNext
}
\def\apROUNDd#1.@#2{\def\XOUT{#1}\edef\XOUT{\apREMzerosR\XOUT}%
\ifnum\apnumX=0 \def\tmpd{}\fi
\ifx\tmpd\empty
\ifx\tmpc\empty \def#2{0}%
\else \edef#2{\ifnum\apnumG<0-\fi\tmpc}\fi
\else\edef#2{\ifnum\apnumG<0-\fi\tmpc.\tmpd}\fi
}
\def\apROUNDe#1.@#2{\apnumC=\apnumD
\apPPs\apROLLa#2{\apnumC}\apPPs\apROUNDa#2{0}\apPPs\apROLLa#2{-\apnumC}%
}
\def\apNORM{\apPPs\apNORMa}
\def\apNORMa#1.@#2{\ifnum\apnumG<0 \def#2{#1}\fi \expandafter\apNORMb\expandafter#2\tmpc@}
\def\apNORMb#1#2#3@{%
\ifx.#2\apnumC=#3\relax \apDIG#1\apnumA \apNORMc#1%
\else \apnumC=#2#3\relax \apDIG#1\relax \apNORMd#1%
\fi
}
\def\apNORMc#1{\advance\apE by-\apnumA \advance\apE by\apnumC
\def\tmpc{-\apnumC}\expandafter\apROLLa#1.@#1%
}
\def\apNORMd#1{\advance\apE by\apnumD \advance\apE by-\apnumC
\def\tmpc{\apnumC}\expandafter\apROLLa\expandafter.#1.@#1%
}
\def\apEadd#1{\ifnum\apE=0 \else\edef#1{#1E\ifnum\apE>0+\fi\the\apE}\apE=0 \fi}
\def\apEnum#1{\ifnum\apE=0 \else\apROLL#1\apE \apE=0 \fi}
%%%%%%%%%%%% Miscelaneous macros, sec. 2.9 in apnum.pdf
\def\apEND{\global\let\apENDx=\OUT
\edef\tmpb{\apSIGN=\the\apSIGN \apE=\the\apE}%
\expandafter\endgroup \tmpb \let\OUT=\apENDx
}
\def\apDIG#1#2{\ifx\relax#2\def\tmpc{}\else #2=0 \def\tmpc{\advance#2 by1 }\fi
\apnumD=0 \expandafter\apDIGa#1..@#1%
}
\def\apDIGa#1{\ifx.#1\csname apDIG\ifnum\apnumD>0 c\else b\fi\expandafter\endcsname
\else \advance\apnumD by1 \expandafter\apDIGa\fi}
\def\apDIGb#1{%
\ifx0#1\advance\apnumD by-1 \tmpc \expandafter\apDIGb
\else \expandafter\apDIGc \expandafter#1\fi
}
\def\apDIGc#1.{\def\tmpd{#1}%
\ifx\tmpc\empty \let\apNext=\apDIGe
\else \def\apNext{\expandafter\apDIGd\tmpd@}%
\fi \apNext
}
\def\apDIGd#1{\ifx@#1\expandafter\apDIGe \else \tmpc \expandafter\apDIGd \fi}
\def\apDIGe#1@#2{%
\ifx@#1@\else % #1=empty <=> the param has no dot, we need to do nothing
\ifnum\apnumD>0 \edef#2{\expandafter\apDIGf#2@}% the dot plus digits before dot
\else \let#2=\tmpd % there are only digits after dot, use \tmpd
\fi\fi
}
\def\apDIGf#1.#2@{#1#2}
\def\apNL{0}
\def\apIVread#1{\expandafter\apIVreadA#1\apNL\apNL\apNL\apNL\apNL@#1}
\def\apIVreadA#1#2#3#4#5\apNL#6@#7{\apnumX=#1#2#3#4\relax \def#7{#5}}
\def\apIVreadX#1#2{\edef\tmpc{\ifcase#1{}{}{}0\or{}{}{}\or{}{}\or{}\fi}%
\expandafter\expandafter\expandafter\apIVreadA\expandafter\tmpc#2\apNL\apNL\apNL\apNL\apNL@#2%
}
\def\apIVwrite#1{\ifnum#1<1000 0\ifnum#1<100 0\ifnum#1<10 0\fi\fi\fi\the#1}
\mathchardef\apIVbase=10000
\def\apIVtrans{\apnumX=\apnumB \divide\apnumB by\apIVbase \multiply\apnumB by-\apIVbase
\advance\apnumB by\apnumX \divide\apnumX by\apIVbase
}
\def\apIVmod#1#2{#2=#1\divide#2by4 \multiply#2by-4 \advance#2by#1\relax
\ifnum#2>0 \else \advance#2by4 \fi
}
\def\apIVdot#1#2{\noexpand\apIVdotA\ifcase#1....\or...\or..\or.\fi #2....@}
\def\apIVdotA#1#2#3#4#5.#6@{\ifx.#1\else#1\fi
\ifx.#2\else#2\fi \ifx.#3\else#3\fi \ifx.#4\else#4\fi\ifx.#5.\else.#5\fi
}
\def\apNUMdigits#1{\expandafter\apNUMdigitsA#1@@@@!}
\def\apNUMdigitsA#1#2#3#4#5!{\ifx@#4\ifx@#3\ifx@#2\ifx@#10\else1\fi \else2\fi \else3\fi \else4\fi}
\def\apADDzeros#1{\edef#1{#10}\advance\apnumZ by-1
\ifnum\apnumZ>0 \expandafter\apADDzeros\expandafter#1\fi
}
\def\apREMzerosR#1{\expandafter\apREMzerosRa#1@0@!}
\def\apREMzerosRa#10@#2!{\ifx!#2!\apREMzerosRb#1\else\apREMzerosRa#1@0@!\fi}
\def\apREMzerosRb#1@{#1}
\def\apREMdotR#1{\expandafter\apREMdotRa#1@.@!}
\def\apREMdotRa#1.@#2!{\ifx!#2!\apREMzerosRb#1\else#1\fi}
\def\apREMfirst#1{\expandafter\apREMfirsta#1@#1}
\def\apREMfirsta#1#2@#3{\def#3{#2}}
\def\apOUTx{\apnumO=7
\edef\apOUTn{\the\apnumL}\edef\apOUTl{\apOUTl\apOUTn,}%
\expandafter\def\csname apOUT:\apOUTn\endcsname{}%
\advance\apnumL by1
}
\def\apOUTs#1,{\ifx.#1\else\csname apOUT:#1\expandafter\endcsname\expandafter\apOUTs\fi}
\def\apINIT{\begingroup \let\do=\apEVALxdo \let\localcounts=\apCOUNTS}
\def\apCOUNTS#1{\ifx;#1\else
\advance\count10 by1 \countdef#1=\count10
\expandafter\apCOUNTS\fi
}
\def\apEVALxdo#1=#2;{#2\let#1=\OUT}
\def\apRETURN#1\apEND{\fi\apEND}
\def\apERR#1{\errmessage{#1}}
{\lccode`\?=`\p \lccode`\!=`\t \lowercase{\gdef\apNOPT#1?!{#1}}}
\def\loop#1\repeat{\def\body{#1\relax\expandafter\body\fi}\body}
%%%%%%%%%%%% Function-like macros, sec. 2.10 in apnum.pdf
\def\ABS#1{\relax % mandatory \relax for "function-like" macros
\evalmdef\OUT{#1}% % evaluation of the input parameter
\ifnum\apSIGN<0 % if (input < 0)
\apSIGN=1 % sign = 1
\apREMfirst\OUT % remove first "minus" from OUT
\fi % fi
}
\def\SGN#1{\relax \evaldef\OUT{#1}\edef\OUT{\the\apSIGN}\apE=0 }
\def\iDIV#1#2{\relax \apINIT % calculation in group
\evalmdef\apAparam{#1}\apEnum\apAparam
\evalmdef\apBparam{#2}\apEnum\apAparam % evaluation of the parameters
\apTOT=0 \apFRAC=0 \apDIV\apAparam\apBparam % integer division
\apEND % end of group
}
\def\iMOD#1#2{\relax \apINIT % calculation in group
\evalmdef\apAparam{#1}\apEnum\apAparam
\evalmdef\apBparam{#2}\apEnum\apBparam % evaluation of the parameters
\apTOT=0 \apFRAC=0 \apDIV\apAparam\apBparam % integer division
\let\OUT=\XOUT % remainder is the output
\apEND % end of group
}
\def\iFLOOR#1{\relax \evalmdef\OUT{#1}\apEnum\OUT \apROUND\OUT0%
\ifnum\apSIGN<0 \ifx\XOUT\empty \else \apPLUS\OUT{-1}\fi\fi
\def\tmp{0}\ifx\tmp\OUT \apSIGN=0 \fi
}
\def\iFRAC#1{\relax
\evalmdef\OUT{#1}\apEnum\OUT \apROUND\OUT0% % preparing the parameter
\ifx\XOUT\empty \def\OUT{0}\apSIGN=0 % empty fraction part means zero
\else \ifnum\apSIGN<0
\edef\XOUT{-.\XOUT}\apPLUS1\XOUT % OUT = 1 - .\XOUT
\else \edef\OUT{.\XOUT}\apSIGN=1 % else OUT = .\XOUT
\fi \fi
}
\def\FAC#1{\relax \apINIT % "function-like" in the group, FAC = factorial
\evalmdef\OUT{#1}\apEnum\OUT % preparing the parameter
\localcounts \N;% % local \newcount
\ifnum\apSIGN<0 \apERR{\string\FAC: argument {\OUT} cannot be negative}\apRETURN\fi
\let\tmp=\OUT \apROUND\tmp0% % test, if parameter is integer
\ifx\XOUT\empty \else \apERR{\string\FAC: argument {\OUT} must be integer}\apRETURN\fi
\N=\OUT\relax % N = param (error here if it is an big integer)
\ifnum\N=0\def\OUT{1}\apSIGN=1 \fi % special definition for factorial(0)
\loop \ifnum \N>2 \advance\N by-1 % loop if (N>2) N--
\apMUL{\OUT}{\the\N}\repeat % OUT = OUT * N , repeat
\apEND % end of group
}
\def\BINOM#1#2{\relax \apINIT % BINOM = {#1 \choose #2} ...
\evalmdef\apAparam{#1}\apEnum\apAparam
\evalmdef\apBparam{#2}\apEnum\apBparam % preparation of the parameters
\localcounts \A \B \C ;% % local \newcounts
\let\OUT=\apBparam \apROUND\OUT0% % test if B is integer
\ifx\XOUT\empty\else\apERR{\string\BINOM: second arg. {\apBparam} must be integer}\apRETURN\fi
\let\OUT=\apAparam \apROUND\OUT0% % test if A is integer
\ifx\XOUT\empty % A is integer:
\A=\apAparam \B=\apBparam % A = #1, B = #2
\C=\A \advance\C by-\B % C = A - B
\ifnum\C>\B \C=\B \fi % if (C > B) C = B fi
\ifnum\A<0 \C=\B % if (A < 0) C = B fi
\else \ifnum\A<\B \def\OUT{0}\apSIGN=0 % if (0 <= A < B) OUT = 0 return
\expandafter\expandafter\expandafter \apRETURN \fi\fi
\def\step{\advance\A by-1 \apMUL\OUT{\the\A}}%
\else \C=\apBparam % A is not integer
\def\step{\let\apBparam\OUT \do\apAparam=\apPLUS\apAparam{-1};%
\let\OUT=\apBparam \apMUL\OUT\apAparam}%
\fi
\ifnum\C=0 \def\OUT{1}\apSIGN=1 \apRETURN\fi
\do\D=\FAC{\the\C};% % D = C!
\let\OUT=\apAparam % OUT = #1
\loop \advance\C by-1 % loop C--
\ifnum\C>0 \step \repeat % if (C > 0) A--, OUT = OUT * A, repeat
\apDIV{\OUT}{\D}% % OUT = OUT / D
\apEND
}
\def\SQRT#1{\relax \apINIT % OUT = SQRT(#1) ...
\evalmdef\A{#1}% % parameter preparation
\localcounts \M \E ;% % local counters
\E=\apE \apE=0
\ifnum\apSIGN=0 \apRETURN\fi % SQRT(0) = 0 (OUT is set to 0 by previous \evaldef)
\ifnum\apSIGN<0 \apERR{\string\SQRT: argument {\A} is out of range}\apRETURN\fi
\ifodd\E \apROLL\A{-1}\advance\E by1 \fi % we need the E representation with even exponent
\let\B=\A \let\C=\A
\apDIG\C\relax \M=\apnumD % M is the number of digits before decimal point
\advance\M by-2 \ifodd\M \advance\M by1 \fi % M = M - 2 , M must be even
\ifx\apSQRTxo\undefined % we need to calculate Xo
\ifnum\M=0 \else \apROLL\B{-\M}\divide\M by2 \fi % shift decimal point by -M, M = M / 2
\apSQRTr\B \let\Xn=\OUT % Xn = estimate of SQRT
\ifnum\M<0 \let\A=\B \fi % if (A < 1) calculate with B where decimal point is shifted
\ifnum\M>0 \apROLL\Xn \M \fi % if (A >= 100) shift the decial point of initial guess
\else \let\Xn=\apSQRTxo \fi
\loop % loop ... Newton's method
\apDIV{\apPLUS{\Xn}{\apDIV{\A}{\Xn}}}{2}% % OUT = (Xn + A/Xn) / 2
\ifx\OUT\Xn \else % if (OUT != Xn)
\let\Xn=\OUT \repeat % Xn = OUT, repeat
\ifnum\M<0 \apROLL\OUT\M \fi % shift the decimal point by M back
\apE=\E \divide\apE by2 % correct the E exponent
\apEND
}
\def\apSQRTr#1{\dimen0=#1pt \apnumB=1 \apnumC=1 \apSQRTra}
\def\apSQRTra{\advance\apnumB by2 \advance\apnumC by\apnumB % B = difference, C = x_i
\ifnum\apnumC>100 \def\OUT{10}\else
\ifdim\dimen0<\apnumC pt \apSQRTrb \else
\expandafter\expandafter\expandafter\apSQRTra\fi\fi
}
\def\apSQRTrb{% x = dimen0, B = x_i - x_{i-1}, C = x_i = i
\ifdim\dimen0<4pt
\ifdim\dimen0>2pt \dimen1=4pt \advance\dimen1 by-\dimen0 \divide\dimen1 by2
\else \dimen1=\dimen0 \advance\dimen1 by-1pt \fi
\dimen1=.080884\dimen1 % dimen1 = additional linear correction
\else \dimen1=0pt \fi
\advance\apnumC by-\apnumB % C = x_{i-1}
\advance\dimen0 by-\apnumC pt % dimen0 = (x - x_{i-1})
\divide\dimen0 by\apnumB % dimen0 = (x - x_{i-1}) / difference
\divide\apnumB by2 % B = i-1 = g(x_{i-1})
\advance\dimen0 by\apnumB pt % dimen0 = g(x_{i-1}) + (x - x_{i-1} / (x_i-x_{i-1})
\advance\dimen0 by\dimen1 % dimen0 += additional linear correction
\edef\OUT{\expandafter\apNOPT\the\dimen0}% OUT = dimen0
}
\def\EXP#1{\relax\apINIT % OUT = EXP(#1) ...
\evalmdef\OUT{#1}\apEnum\OUT % OUT = #1
\localcounts \N \K ;%
\ifnum\apSIGN=0 \def\OUT{1}\apSIGN=1 \apRETURN \fi
\edef\digits{\the\apFRAC}\advance\apFRAC by4
\edef\signX{\the\apSIGN}%
\ifnum\apSIGN<0 \apSIGN=1 \apREMfirst\OUT \fi % remove "minus" sign
\def\testBig ##1##2##3\relax##4{\ifx##1.\apXfalse \else
\ifx##2.\ifnum##1<4 \apXfalse \else \apXtrue \fi \else \apXtrue
\fi \fi \ifapX}%
\expandafter\testBig \OUT.\relax
\iftrue \apEXPb \else \apEXPa \fi % OUT = e^OUT
\ifnum\signX<0 \K=-\apE \apDIV 1\OUT \apE=\K \fi % if (signX < 0) OUT = 1 / OUT
\apSIGN=1 % EXP is always positive
\apEND
}
\def\apEXPa{%
\def\testDot ##1##2\relax##3{\ifx##1.}%
\K=0 \N=0 % K = 0, N = 0
\loop \expandafter \testDot\OUT \relax % loop if (OUT >= 1)
\iftrue \else % OUT = OUT/2
\apDIV\OUT{2}% % K++
\advance\K by1 % repeat
\repeat % oriOUT = 2^K * OUT, OUT < 1
\advance\apFRAC by\K
\def\S{1}\def\Sn{1}\N=0 \let\X=\OUT % S = 1, Sn = 1, N = 0, X = OUT
\loop \advance\N by1 % loop N++
\do\Sn=\apDIV{\apMUL\Sn\X}{\the\N};% % Sn = Sn * X / N
\apTAYLOR\iftrue \repeat % S = S + Sn (... Taylor)
\N=0
\loop \ifnum\N < \K % loop if (N < K)
\apPOW\OUT{2}% % OUT = OUT^2
\advance\N by1 \repeat % N++
\apFRAC=\digits\relax \apROUND\OUT\apFRAC
}
\def\apTAYLOR#1{\ifnum\apSIGN=0 \let\OUT=\S \else \apPLUS\S\Sn \let\S=\OUT }
\def\apEXPb{%
\let\X=\OUT \apLNtenexec \apDIV\X\apLNten \let\D=\OUT
\apROUND\D{0}% % D = floor( X/ln(X) )
\ifnum\D<\apEX \advance\apFRAC by\D \relax \apLNtenexec \fi
\EXP{\X-\D*\apLNten}% % mantissa = EXP(X-D*LN(10))
\ifnum\D<\apEX \apROLL\OUT\D \apE=0 \else \apE=\D \relax \fi
\apFRAC=\digits \apROUND\OUT\apFRAC % OUT = mantissa * 10^D
}
\def\LN#1{\relax \apINIT % OUT = LN(#1) ...
\evalmdef\X{#1}% % X = #1
\localcounts \M \N \E;%
\E=\apE
\edef\digits{\the\apFRAC}\advance\apFRAC by4
\ifnum\apSIGN>0 \else \apERR{\string\LN: argument {\X} is out of range}\apRETURN\fi
\apDIG\OUT\relax \M=\apnumD % find M: X = mantissa * 10^M
\ifnum\M>-\E \def\sgnout{1}\else % if X in (0,1):
\def\sgnout{-1}% % sgnout = -1
\do\X=\apDIV 1\X;\E=-\E % X = 1/X
\apDIG\OUT\relax \M=\apnumD % find M: X = mantissa * 10^M
\fi % else sgnout = 1
\advance\M by-1 % M = M - 1
\ifnum\M=0 \else\apROLL\X{-\M}\fi % X = X * 10^(-M), now X in (1,10)
\advance\M by\E % M = M + E (sientific format of numbers)
\do\lnX=\apLNr\X;% % lnX = LN(X) ... roughly estimate
\do\A=\apDIV\X{\EXP\lnX};% % A = X / EXP(lnX) ... A =approx= 1
\apLNtaylor % OUT = LN(A)
\do\LNOUT=\apPLUS\OUT\lnX;% % LNOUT = OUT + LNrOUT
\ifnum\M>0 % if M > 0
\apLNtenexec % LNtenOUT = ln(10)
\apPLUS\LNOUT{\apMUL{\the\M}{\apLNten}}% OUT = LNOUT + M * LNten
\fi
\ifnum\apSIGN=0 \else \apSIGN=\sgnout \fi % if (OUT != 0) apSIGN = saved sign
\apROUND\OUT\digits % round result to desired precision
\ifnum\apSIGN<0 \xdef\OUT{-\OUT}\else \global\let\OUT=\OUT \fi
\apEND
}
\def\apLNtaylor{%
\apDIV{\apPLUS{\A}{-1}}{\apPLUS{\A}{1}}% % OUT = (A-1) / (A+1)
\ifnum\apSIGN=0 \def\OUT{0}\else % ln 1 = 0 else:
\let\Sn=\OUT \let\Kn=\OUT \let\S=\OUT % Sn = OUT, Kn = OUT, S = OUT
\apPOW\OUT{2}\apROUND\OUT\apFRAC \let\XX=\OUT % XX = OUT^2
\N=1 % N = 1
\loop \advance\N by2 % loop N = N + 2
\do\Kn=\apMUL\Kn\XX\apROUND\OUT\apFRAC;% Kn = Kn * XX
\do\Sn=\apDIV\Kn{\the\N};% % Sn = Kn / N
\apTAYLOR\iftrue \repeat % S = S + Sn (Taylor)
\apMUL\S{2}% % OUT = 2 * OUT
\fi
}
\def\apLNr#1{\dimen0=#1pt \apnumC=1
\apLNra {0}{.69}{1.098}{1.386}{1.609}{1.791}{1.9459}{2.079}{2.197}{\apLNrten}{}\relax
}
\def\apLNra #1#2{\advance\apnumC by1
\ifx\relax#2\relax \let\OUT=\apLNrten \let\apNext=\relax
\else
\ifdim\dimen0<\apnumC pt % linear interpolation:
\advance\dimen0 by-\apnumC pt \advance\dimen0 by1pt % dimen0 = x - x_{i-1}
\dimen1=#2pt \advance\dimen1 by-#1pt % dimen1 = f(x_i) - f(x_{i-1})
\dimen1=\expandafter\apNOPT\the\dimen0 \dimen1 % dimen1 = (x - x_{i-1}) * dimen1
\advance\dimen1 by#1pt % dimen1 = f(x_{i-1}) + dimen1
\edef\OUT{\expandafter\apNOPT\the\dimen1}% % OUT = dimen1
\def\apNext##1\relax{}%
\else \def\apNext{\apLNra{#2}}%
\fi\fi \apNext
}
\def\apLNrten{2.302585} % apLNrten = ln 10 (roughly)
\def\apLNtenexec{% % OUT = ln 10 ...
\expandafter\ifx\csname LNten:\the\apFRAC\endcsname \relax
\begingroup \apTOT=0
\do\A=\apDIV{10}{\EXP\apLNrten};% % A = 10 / exp(LNrten)
\apLNtaylor % OUT = ln A
\apPLUS\OUT\apLNrten % OUT = OUT + LNrten
\global\expandafter\let\csname LNten:\the\apFRAC\endcsname=\OUT
\endgroup
\fi
\expandafter\let\expandafter \apLNten \csname LNten:\the\apFRAC\endcsname
}
\def\apPIvalue{3.141592653589793238462643383279}
\def\apPIdigits{30}
\def\apPIexec{%
\expandafter\ifx\csname apPI:\the\apFRAC\endcsname \relax \apPIexecA \else
\expandafter\let\expandafter\apPI\csname apPI:\the\apFRAC\endcsname
\expandafter\let\expandafter\apPIhalf\csname apPIh:\the\apFRAC\endcsname
\fi
}
\def\apPIexecA{%
\ifnum\apPIdigits<\apFRAC \apPIexecB \fi
\let\apPI=\apPIvalue
\ifnum\apPIdigits>\apFRAC \apROUND\apPI\apFRAC \fi
\apnumP=\apTOT \apTOT=0 \apDIV\apPI2\let\apPIhalf=\OUT \apTOT=\apnumP
\global\expandafter\let\csname apPI:\the\apFRAC\endcsname=\apPI
\global\expandafter\let\csname apPIh:\the\apFRAC\endcsname=\apPIhalf
}
\def\apPIexecB{\apINIT
\localcounts \N \a \c;%
\apTOT=0 \advance\apFRAC by2
\def\apSQRTxo{800.199975006248}% initial value for Newton method for SQRT
\SQRT{640320}%
\let\sqrtval=\OUT
\N=0 \def\An{1}\def\Bn{1}\def\Cn{1}\def\S{13591409}%
\loop
\advance\N by 1
\a=\N \multiply\a by6 \advance\a by-1 \c=\a
\advance\a by-2 \multiply\c by\a % An = An * 8 * (6N-5) *
\advance\a by-2 \multiply\a by8 % * (6N-3) * (6N-1)
\apMUL\An{\apMUL{\the\a}{\the\c}}\let\An=\OUT
\c=\N \multiply\c by\N % Bn = Bn * n^3
\apMUL\Bn{\apMUL{\the\c}{\the\N}}\let\Bn=\OUT
\apMUL\Cn{-262537412640768000}\let\Cn=\OUT % Cn = Cn * K3
\apDIV{\apMUL\An{\apPLUS{13591409}{\apMUL{545140134}{\the\N}}}}{\apMUL\Bn\Cn}%
\let\Sn=\OUT % Sn = An * (K1 + K2 * N) / (Bn * Cn)
\apTAYLOR \iftrue \repeat
\advance\apFRAC by-2
\apDIV{\apMUL{\sqrtval}{53360}}\S
\global\let\apPIvalue=\OUT
\xdef\apPIdigits{\the\apFRAC}%
\apEND
}
\def\PI{\relax \apPIexec \let\OUT=\apPI}
\def\PIhalf{\relax \apPIexec \let\OUT=\apPIhalf}
\def\SIN{\relax \let\apSINCOSx=\apSINx \apSINCOSa}
\def\COS{\relax \let\apSINCOSx=\apCOSx \apSINCOSa}
\def\apSINCOSa#1{\apINIT
\advance\apFRAC by3
\evalmdef\X{#1}\apEnum\X
\def\signK{1}\apSINCOSo\apCOSx
\ifnum\apSIGN<0 \apREMfirst\X \def\sign{-}\else\def\sign{+}\fi
\ifx\apSINCOSx\apCOSx \def\sign{+}\fi
\edef\apFRACsave{\the\apFRAC}%
\apPIexec
\apFRAC=0 \apDIV\X\apPI % OUT = X div PI
\ifnum\apSIGN=0 \apSIGN=1 \else
\let\K=\OUT
\do\X=\apPLUS\X{-\apMUL\K\apPI};% X := X - K * PI
\apROLL\K{-1}\apROUND\K{0}%
\ifodd 0\XOUT\space \def\signK{-1}\else\def\signK{1}\fi
\fi
\apSINCOSo\apCOSx
\apFRAC=\apFRACsave \relax
\do\XmPIh=\apPLUS\X{-\apPIhalf};% XmPIh = | X - PI/2 |
\apSINCOSo\apSINx
\ifnum\apSIGN<0 \apREMfirst\XmPIh
\else % X in (PI/2, PI)
\do\X=\apPLUS\apPI{-\X};%
\ifx\apSINCOSx\apCOSx \apSIGN=-\signK \edef\signK{\the\apSIGN}\fi
\fi % X in (0, PI/2):
\apMINUS\X{.78}% % OUT = X - cca PI/4
\ifnum\apSIGN<0 \else % if X in (PI/4, PI/2) :
\let\X=\XmPIh % X = | X - PI/2 |; SIN <-> COS
\ifx\apSINCOSx\apSINx \let\apSINCOSx=\apCOSx \else \let\apSINCOSx=\apSINx \fi
\fi
\localcounts \N \NN;%
\do\XX=\apPOW\X{2}\ROUND\OUT\apFRAC;%
\apSINCOSx % X in (0, PI/4), initialize Taylor SIN X or COS X
\loop
\advance\N by1 \NN=\N
\advance\N by1 \multiply\NN by\N
\do\Sn=\apDIV{\apMUL\Sn\XX}{-\the\NN};% Sn = - Sn * X^2 / N*(N+1)
\apTAYLOR \iftrue\repeat
\apSIGN=\sign\signK
\ifnum\apTOT=0 \advance\apFRAC by-3 \else \apFRAC=\apTOT \fi
\ifnum\apFRAC<0 \apFRAC=-\apFRAC \fi
\apROUND\OUT\apFRAC
\def\X{0}\ifx\OUT\X \apSIGN=0 \fi
\ifnum\apSIGN<0 \edef\OUT{-\OUT}\fi
\apEND
}
\def\apSINx{\let\S=\X \N=1 \let\Sn=\X}
\def\apCOSx{\def\S{1}\N=0 \let\Sn=\S}
\def\apSINCOSo#1{\ifnum\apSIGN=0 \ifx#1\SCgo \apSIGN=\signK \let\OUT=\signK \fi \apRETURN\fi}
\def\TAN#1{\relax \apINIT
\advance\apFRAC by3
\evalmdef\X{#1}\apEnum\X
\advance\apFRAC by-3
\do\denom=\COS\X;%
\ifnum\apSIGN=0 \apERR{\string\TAN: argument {\X} is out of range}\apRETURN\fi
\SIN\X
\apDIV{\SIN\X}\denom
\apEND
}
\def\ATAN#1{\relax \apINIT
\advance\apFRAC by3
\evalmdef\X{#1}\apEnum\X
\ifnum\apSIGN=0 \def\OUT{0}\apRETURN\fi
\ifnum\apSIGN<0 \def\sign{-}\apREMfirst\X \else\def\sign{}\fi
\let\tmp=\X \apDIG\tmp\relax
\ifnum\apnumD>0 % if X > 1:
\apPIexec % OUT = apPIhalf - apATANox
\def\tmp{1}\ifx\tmp\X \apDIV\apPIhalf2\else \apATANox \apPLUS\apPIhalf{-\OUT}\fi
\else % else
\do\X=\apDIV{1}\X;% X := 1/X
\apATANox % OUT = apATANox
\fi
\ifnum\apTOT=0 \advance\apFRAC by-3 \else \apFRAC=\apTOT \fi
\ifnum\apFRAC<0 \apFRAC=-\apFRAC \fi
\apROUND\OUT\apFRAC
\ifx\sign\empty\apSIGN=1 \else \edef\OUT{-\OUT}\apSIGN=-1 \fi
\apEND
}
\def\apATANox{%
\localcounts \N;%
\do\XX=\apPLUS{1}{\apPOW\X{2}}\apROUND\OUT\apFRAC;% XX = 1 + X^2
\do\Sn=\apDIV\X\XX \apROUND\OUT\apFRAC;% % Sn = X / (1+X^2)
\N=1 \let\S=\Sn
\loop
\advance\N by1
\do\Sn=\apMUL{\the\N}\Sn;%
\advance\N by1
\do\Sn=\apDIV\Sn{\apMUL{\the\N}\XX};% Sn = Sn * N / ((N+1) * (1+X^2))
\apTAYLOR \iftrue \repeat
}
\def\ASIN#1{\relax \apINIT
\evalmdef\X{#1}\apEnum\X \edef\sign{\the\apSIGN}%
\apPLUS 1{-\apPOW\X2}% OUT = 1 - X^2
\ifnum\apSIGN<0 \apERR{\string\ASIN: argument {\X} is out of range}\apRETURN\fi
\do\sqrt=\SQRT\OUT;% sqrt = SRQT {1 - X^1}
\ifnum\apSIGN=0 \apPIexec
\ifnum\sign<0 \edef\OUT{-\apPIhalf}\apSIGN=-1 % ASIN(-1) = -PI/2
\else \let\OUT=\apPIhalf \apSIGN=1 \fi % ASIN(1) = PI/2
\apRETURN \fi
\ATAN{\X/\sqrt}% OUT = arctan ( X / SQRT {1 - X^2} )
\apEND
}
\def\ACOS#1{\relax \apPIexec \apPLUS\apPIhalf{-\ASIN{#1}}}
%%%%%%%%%%%% Printing expressions, sec 2.11 in apnum.pdf
\def\eprint#1#2{\bgroup \apnumA=0 \apnumE=1 \apEVALb#1\limits
\let\apEPe=\relax \apEPi #2\tmpb \apEPe \egroup
}
\def\apEPi{\let\apPLUS=\apEPplus \let\apMINUS=\apEPminus
\let\apMUL=\apEPmul \let\apDIV=\apEPdiv \let\apPOWx=\apEPpow \def\apPPn##1{##1}%
\let\EXP=\apEPexp \def\LN{\apEPf{ln}}\let\SQRT=\apEPsqrt
\def\SIN{\apEPf{sin}}\def\COS{\apEPf{cos}}\def\TAN{\apEPf{tan}}%
\def\ASIN{\apEPf{arcsin}}\def\ACOS{\apEPf{arccos}}\def\ATAN{\apEPf{arctan}}%
\let\PI=\pi \def\PIhalf{{\pi\over2}}%
\let\ABS=\apEPabs \let\FAC=\apEPfac \let\BINOM=\apEPbinom
\let\SGN=\apEPsgn \let\iDIV=\apEPidiv \let\iMOD=\apEPimod
\let\iFLOOR=\apEPifloor \let\iFRAC=\apEPifrac
\let\apEPk=\empty \let\apEPy=\empty \def\apEPx{.}%
\let\apEPi=\relax \apEPj
}
\def\apEPj{}
\def\apEPplus#1#2{\apEPp{#1}{?...}+\apEPp{#2}{!...}}
\def\apEPminus#1#2{\apEPp{#1}{?...}-\apEPp{#2}{!!..}}
\def\apEPmul#1#2{\def\tmpa{#1}\def\tmpb{-1}%
\ifx\tmpa\tmpb \if\apEPx!\left(-\apEPp{#2}{!!..}\right)\else
-\apEPp{#2}{!!..}\fi
\else \apEPp{#1}{?!..}\apMULop \apEPp{#2}{!!..}\fi
}
\def\apEPdiv#1#2{{\def\apEPx{.}#1}\over{\def\apEPx{.}#2}}
\def\apEPpow#1#2{%
\let\apEPy=\empty \apEPpowa{#1}\end{#2}%
\ifx\apEPy\empty \apEPp{#1}{!!!!}^{\def\apEPx{.}#2}\else#1\fi
}
\def\apEPpowa#1{\expandafter\apEPpowb#1;}
\def\apEPpowb#1#2;\end#3{\ifx#1\apEPf \def\apEPy{\let\apEPy=\empty\def\apEPx{.}#3}\fi}
\def\apEPf#1#2{\begingroup
\mathop{\rm#1}\nolimits
\ifx\apEPy\empty \else ^{\apEPy}\let\apEPy=\empty \fi
\def\apEPk{\mskip-\thinmuskip}%
\def\apEPx{.}%
\eprint{#2}{\expandafter\apEPb}\endgroup
}
\def\apEPb#1#2{\def\next{\apEPk\left(\def\apEPe{\right)}}%
\ifx\apPPn#1\expandafter\apEPd#2.\end{}{\let\next=\relax}.\fi
\ifx\apDIV#1\let\next=\relax \fi
\next\let\apEPk=\empty #1{#2}%
}
\def\apEPp#1#2{\apEPq#1\end\bgroup{\left(\def\apEPx{.}}#2#1\apEPq#1\end\egroup{\right)}#2}
\def\apEPq#1#2\end#3#4#5#6#7#8{
\ifx#5!\def\apEPx{!}\fi
\ifx#1\apEPplus \ifx#6!#4\else#3\fi\else
\ifx#1\apEPminus \ifx#6!#4\else#3\fi\else
\ifx#1\apEPmul \ifx#7!#4\else#3\fi\else
\ifx#1\apEPdiv \ifx#8!#4\else#3\fi\else
\ifx#1\apEPpow \ifx#8!#4\else#3\fi\else
\expandafter\apEPd#1.\end#3{}\apEPx\fi\fi\fi\fi\fi
}
\def\apEPd#1#2\end#3#4#5{\ifx-#1\if#5!\ifx#3\bgroup\left(\else\right)\fi\fi\else#4\fi}
\let\apMULop=\cdot
\def\apEPabs#1{\left|\eprint{#1}{}\right|}
\def\apEPfac#1{\eprint{#1}{\expandafter\apEPb}\,!}
\def\apEPbinom#1#2{{\eprint{#1}{}\choose\eprint{#2}{}}}
\def\apEPsqrt#1{\sqrt{\eprint{#1}{}}}
\def\apEPexp#1{{\rm e}^{\eprint{#1}{}}}
\def\apEPsgn#1{\mathop{\rm sign}\eprint{#1}{\expandafter\apEPb}}
\def\apEPdivmod#1#2#3{\left[\eprint{#2}{\expandafter\apEPb}%
\mathbin{\rm #1}\eprint{#3}{\expandafter\apEPb}\right]}
\def\apEPidiv{\apEPdivmod{div}}
\def\apEPimod{\apEPdivmod{mod}}
\def\apEPifloor#1{\left\lfloor\eprint{#1}{}\right\rfloor}
\def\apEPifrac#1{\left\{\eprint{#1}{}\right\}}
\def\corrnum#1{\edef#1{\expandafter\apEPc#1\end}}
\def\apEPc#1#2\end{\ifx#1-{-}\apEPc#2\end\else \ifx#1.0.#2\else #1#2\fi\fi}
%%%%%%%%%%%% Conclusion, sec. 2.12 in apnum.pdf
\let\PLUS=\apPLUS \let\MINUS=\apMINUS \let\MUL=\apMUL \let\DIV=\apDIV \let\POW=\apPOW
\let\SIGN=\apSIGN \let\ROUND=\apROUND \let\NORM=\apNORM \let\ROLL=\apROLL
\ifx\documentclass\undefined \else % please, don't remove this message
\message{WARNING: the author of apnum package recommends: Never use LaTeX.}\fi
\catcode`\@=\apnumZ
\endinput
1.0 <Nov 2014> - First version released
1.1 <Jan 2015>
- POW implemented more simple (by base 2 of exponent)
- \next renamed in order to avoid name conflict
1.2 <May 2015> - .5+.5=.1 bug fixed
1.3 <Dec 2015> -
- \apPPn corrected (empty \OUT bug fixed)
- \apEVAL: spaces ignored between parameters of function-like macros
- \apEVAL: in one group, \apEND introduced, \apOUTtmpb removed
- \apPLUS, etc. instead \PLUS introduced
- \apSTRIPfirst introduced
- \apEVALone, \apEVALtwo removed
- \addE renamed to \apEadd
- \ROLL, \NORM, \ROUND renamed to \apROLL, \apNORM, \apROUND
- \apREV removed
- \ABS, \iDIV, \iMOD, \iROUND, \iFRAC, \FAC rewriten
- \XOUT is empty, no "0000" after \apROLL\a0 (2.0000 bug fixed)
- \def#1{} corrected
- \localcount after \evaldef in order to avoid name conflict
- \PI added
- \apEnorm to \apEnum renamed
1.4 <Dec 2015>
- \ATAN, \ASIN, \ACOS added
- \SIN, \COS, \TAN added
- \apTOT=0 by default
1.4a \end -> \limits, internal change in \evaldef because LaTeX redefines \end
1.4b \apSINCOS: \ifx\OUT eq 0 added after \apROUND (bug fixed).
1.5 <Jan 2016>
- \iROUND replaced by \iFLOOR, \iFRAC corrected (for negative numbers)
- \eprint introduced
1.6 <Feb 2016>
- \evalmdef introduced
- \EXP for arg>=4 rewritten, \apEX register introduced
1.7 <Apr 2018>
- \eprint: bug removed (round brackets around negative constants)