summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/source
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2013-07-17 22:13:10 +0000
committerKarl Berry <karl@freefriends.org>2013-07-17 22:13:10 +0000
commiteec0a21648f82e9f409ed8151e1bc83f34bfec79 (patch)
treedd59250c0f9496a72010ad5b1097e0199888f06f /Master/texmf-dist/source
parentd9add100de8f268f2b408533a41cc5624ef1ee7a (diff)
minifp (17jul13)
git-svn-id: svn://tug.org/texlive/trunk@31221 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/source')
-rw-r--r--Master/texmf-dist/source/generic/minifp/minifp.dtx1388
1 files changed, 716 insertions, 672 deletions
diff --git a/Master/texmf-dist/source/generic/minifp/minifp.dtx b/Master/texmf-dist/source/generic/minifp/minifp.dtx
index 77f56099f70..c3621aef2ba 100644
--- a/Master/texmf-dist/source/generic/minifp/minifp.dtx
+++ b/Master/texmf-dist/source/generic/minifp/minifp.dtx
@@ -13,8 +13,8 @@
% minifp has maintenance status "author-maintained". The Current Maintainer
% is Daniel H. Luecking. The Base Interpreter is TeX (plain TeX or LaTeX).
%<*driver|sty>
-\def\MFPfiledate{2013/02/01}%
-\def\MFPfileversion{0.92}%
+\def\MFPfiledate{2013/05/28}%
+\def\MFPfileversion{0.95}%
%</driver|sty>
%
%<*driver>
@@ -49,8 +49,11 @@
\let\env\file
\def\sgn{\mathop{\mathrm{sgn}}\nolimits}
% \op is for abstract operations (e.g., \op{add}) as opposed to
-% the macro that performs it (e.g., \cs{Radd}).
+% the macro that performs it (e.g., \cs{Radd}). And \reg is for
+% a "register" (e.g., the 3 macros \MFP@x@Sgn, \MFP@x@Int and \MFP@x@Frc)
+% conceived of as a single entity.
\let\op\textit
+\def\reg#1{$#1$}
% The occasional bare \tt braces
\renewcommand\{{\char`\{}
\renewcommand\}{\char`\}}
@@ -60,7 +63,7 @@
\makeatletter
\newcommand\bsl{{\mytt\@backslashchar}}
-% Stupid lists!
+% better lists
\def\@listi{\leftmargin\leftmargini
\parsep \z@ \@plus\p@ \@minus\z@
\topsep 4\p@ \@plus\p@ \@minus2\p@
@@ -70,11 +73,10 @@
\renewcommand\labelitemii{\textasteriskcentered}
\renewcommand\labelitemiii{\textperiodcentered}
\leftmargini\parindent
-% Stupid index!
+% better index
\def\usage#1{\textrm{#1}}
\def\index@prologue{\section*{Index}\markboth{Index}{Index}%
- Numbers refer to the page(s) where the corresponding entry is described.
-}
+ Numbers refer to the page(s) where the corresponding entry is described.}
\def\IndexParms{%
\parindent \z@ \columnsep 15pt
\parskip 0pt plus 1pt
@@ -102,7 +104,7 @@
\end{document}
%</driver>
%\fi
-% \CheckSum{3339}
+% \CheckSum{3402}
% \CharacterTable
% {Upper-case \A\B\C\D\E\F\G\H\I\J\K\L\M\N\O\P\Q\R\S\T\U\V\W\X\Y\Z
% Lower-case \a\b\c\d\e\f\g\h\i\j\k\l\m\n\o\p\q\r\s\t\u\v\w\x\y\z
@@ -138,13 +140,13 @@
% to make simple calculations with real numbers. What \TeX{} provides is
% far too limited. In fact, its only native user-level support for real
% numbers is as factors for dimensions. For example one can ``multiply''
-% $3.1\times 0.2$ by \verb$\dimen0=0.2pt \dimen0=3.1\dimen0 $.
+% $3.1\times 0.2$ by the following: \verb$\dimen0=0.2pt \dimen0=3.1\dimen0$.
%
% Unfortunately \TeX{} stores dimensions as integer multiples of the
-% ``scaled point'' (\dim{sp}) with \dim{sp}${}=2^{-16}$\dim{pt}, and
+% ``scaled points'' (\dim{sp}) with \dim{sp}${}=2^{-16}$\dim{pt}, and
% therefore \dim{.2pt} is approximated by $\frac{13107}{65536}$, which is
% not exact. Then mutiplying by $3.1$ produces $\frac{40631}{65536}$. If
-% we ask for five digit accuracy, this produces $0.61998$\dim{pt} and not the
+% we ask \TeX{} to display this, it produces $0.61998$\dim{pt} and not the
% exact value $0.62$. This is sufficiently accurate for positioning
% elements on a page, but not for displaying automatically computed axis
% labels if five digit accuracy is needed.
@@ -205,30 +207,35 @@
% exceed the allowed eight digits) is always possible, but is much more
% likely with multiplication and division.
%
-% Multiplication is carried out internally to an exact 16-digit answer,
-% which is then rounded to an 8-digit result. Overflow (more than 8
-% digits in the integer part) is discarded. Division is internally
-% carried to nine digits after the decimal, which is then also rounded to
-% an 8-digit result.
+% Multiplication is carried out internally to an exact answer, with 16
+% digits on each side of the decimal point. The underflow digits (places 9
+% through 16 after the decimal point) are used to round to an 8-digit
+% result. Overflow digits (those to the left of the lowest 8 in the
+% integer part) are discarded, usually without warning. Division is
+% internally carried to nine digits after the decimal, which is then also
+% rounded to an 8-digit result. Overflow digits are ignored for division
+% also.
%
% We supply two kinds of operations in this package. There are stack-based
-% operations, in which the operands are popped from a stack and the
-% results pushed onto it, and argument-based, in which the operands (and a
+% operations, in which the operands are \op{popped} from a stack and the
+% results \op{pushed} onto it, and argument-based, in which the operands (and a
% macro to hold the result^^A
-% \footnote{Unlike most other packages for floating point
+% \footnote{Unlike most other packages for decimal
% arithmetic, \mfp{} puts the macro to hold the result
% last. This allows the calculation to be performed before the
-% macro is even read, and makes it somewhat easier for the
+% macro is even read, and this makes it somewhat easier for the
% stack- and argument-based versions to share code.}^^A
% ) are arguments of a macro. Both types load the arguments into internal
-% macros (think ``registers''), then call internal commands
+% macros (think of them as ``registers''), then call internal commands
% (think ``microcode'') which return the results in internal macros.
-% These results are then pushed onto the stack (stack-based operations) or
-% stored in a supplied macro argument (think ``variable'').
+% These results are then \op{pushed} onto the stack (stack-based
+% operations) or stored in a supplied macro argument (think ``variable'').
+% The difference lies entirely in where the operands come from (arguments
+% or stack) and where they go (macro or stack).
%
% The stack is implemented as an internal macro which is redefined with
-% each command. The binary operations act on the last two pushed objects
-% in the order they were pushed. For example, the sequence ``\op{push} 5,
+% each command. The binary operations act on the last two \op{pushed} objects
+% in the order they were \op{pushed}. For example, the sequence ``\op{push} 5,
% \op{push} 3, \op{subtract}'' performs $5-3$ by popping $3$ and $5$ into
% registers (thereby removing them from the stack), subtracting them
% and then pushing the result ($2$) onto the stack.
@@ -269,9 +276,14 @@
% \def\MFPextra{} \input minifp.sty \end{verbatim}
% The extras can also be loaded by means of the command
% \cs{MFPloadextra}, issued after \file{minifp.sty} is loaded.
+% As of version 0.95 \file{mfpextra} can be directly \cs{input}.
+% It will detect whether \file{minifp.sty} has been loaded and input it
+% if not. This will work only in plain \TeX{}.
%
% If the extra operations are not needed, some memory and time might be
-% saved by using \file{minifp.sty} alone.
+% saved by using \file{minifp.sty} alone. I have not seriously tried to
+% keep \file{mfpextra.tex} as small or fast as possible, but I do try
+% to improve the accuracy when I can.
%
% As previously mentioned, each of these operations come in two versions:
% a version that acts on operands and stores the result in a macro, and a
@@ -283,7 +295,7 @@
% because it is possible that stacks of other types will be implemented in
% the future.
%
-% For example, \verb$\MFPadd{1.3}{3.4}\X$ will add $1.20000000$ to
+% For example, \verb$\MFPadd{1.2}{3.4}\X$ will add $1.20000000$ to
% $3.40000000$ and then define \cs{X} to be the resulting
% \texttt{4.60000000}. These operand forms do not alter or even address
% the stack in any way. The stack-based version of the same operation
@@ -316,10 +328,10 @@
%
% In the following tables, an argument designated
% \meta{num} can be any decimal real number with at most 8
-% digits on each side of the decimal point, or they can be macros that
-% contain such a number. If the decimal dot is absent, the fractional part
-% will be taken to be $0$, if the integer part or the fractional part is
-% absent, it will be taken to be $0$. (One consequence of these rules is
+% digits on each side of the decimal point, or it can be a macro that
+% contains such a number. If the decimal dot is absent, the fractional part
+% will be taken to be zero, if the integer part or the fractional part is
+% absent, it will be taken to be zero. (One consequence of these rules is
% that all the following arguments produce the same internal
% representation of zero: \marg{0.0}, \marg{0.}, \marg{.0},
% \marg{0}, \marg{.}, and \marg{}\,.) Spaces may appear anywhere in the
@@ -381,7 +393,7 @@
% Stores $|$\meta{num}$|$ in \cs{macro}.\\
% \SpecialUsageIndex{\MFPdbl}^^A
% \cs{MFPdbl}\mmarg{num}\cs{macro}&
-% Stores 2\meta{num} in \cs{macro}.\\
+% Stores $2\times{}$\meta{num} in \cs{macro}.\\
% \SpecialUsageIndex{\MFPhalve}^^A
% \cs{MFPhalve}\mmarg{num}\cs{macro}&
% Stores \meta{num}/2, rounded to 8 places after the decimal point, in
@@ -389,12 +401,12 @@
% \SpecialUsageIndex{\MFPint}^^A
% \cs{MFPint}\mmarg{num}\cs{macro}&
% Replaces the part of \meta{num} after the decimal point with zeros
-% (keeps the sign unless the result is $0$) and stores the result in
+% (keeps the sign unless the result is zero) and stores the result in
% \cs{macro}.\\
% \SpecialUsageIndex{\MFPfrac}^^A
% \cs{MFPfrac}\mmarg{num}\cs{macro}&
-% Replaces the part of \meta{num} before the decimal point with $0$
-% (keeps the sign unless the result is $0$) and stores the result in
+% Replaces the part of \meta{num} before the decimal point with zero
+% (keeps the sign unless the result is zero) and stores the result in
% \cs{macro}.\\
% \SpecialUsageIndex{\MFPfloor}^^A
% \cs{MFPfloor}\mmarg{num}\cs{macro}&
@@ -432,7 +444,8 @@
% The command \cs{MFPzero} is useful for ``macro programs''. If you want
% to do something to a number depending on the outcome of a test, you may
% occasionally want to simply absorbed the number and output a default
-% result. (There are more efficient ways to simply store $0$ in a macro.)
+% result. This is more efficient than multiplying by zero (but less
+% efficient than simply defining the \cs{macro} to be zero.)
%
% Note that one could easily double, halve, square, increment,
% decrement or invert a \meta{num} using the binary versions of
@@ -445,17 +458,16 @@
% use the two argument versions, \cs{MFPmul}\mmarg{num}\marg{.5} is faster than
% \cs{MFPdiv}\mmarg{num}\marg{2}.
%
-% There is one command that takes no argument:
+% There is one command that takes no argument and returns no value:
%
% \medskip
% \centerline{%
% \begin{tabular}{lp{3.4in}}
-% \textit{Nullary Operations}&\\[3pt]
+% \textit{Do Nothing}&\\[3pt]
% \hline\hline
% \textbf{Command}&\textbf{operation}\\
% \hline
-% \SpecialUsageIndex{\MFPnoop}^^A
-% \cs{MFPnoop}& Does nothing.
+% \SpecialUsageIndex{\MFPnoop}\cs{MFPnoop}& Does nothing.
% \end{tabular}}
%
% \bigskip
@@ -502,7 +514,7 @@
% \medskip
% Issuing \verb$\MFPchk{\X}$ will check the sign of the number stored in
% the macro \cs{X}. Then \verb$\IFneg{A}{B}$ will produce `\verb$A$' if it
-% is negative and `\verb$B$' if it is $0$ or positive. Similarly,
+% is negative and `\verb$B$' if it is zero or positive. Similarly,
% \verb$\MFPcmp{\X}{1}$ will compare the number stored in \cs{X} to $1$.
% Afterward, \verb$\IFlt{A}{B}$ will produce `\verb$A$' if \cs{X} is less
% than $1$ and `\verb$B$' if \cs{X} is equal to or greater than $1$.
@@ -689,20 +701,20 @@
% Decreases by $1$. Slightly more efficient than the equivalent
% subtraction.\\
% \SpecialUsageIndex{\Rzero}\cs{Rzero}&
-% Replaces the number with $0$. Slightly more convenient than the
+% Replaces the number with zero. Slightly more convenient than the
% equivalent \cs{Rpop}\cs{X} followed by a \cs{Rpush}\marg{0}.\\
% \end{tabular}}
%
%\bigskip
%
%
-% There is one nullary operation, which does not read the stack nor
-% change it.
+% There is one operation, which does not read the stack nor change it
+% (nor do anything else).
%
% \medskip
% \centerline{%
% \begin{tabular}{lp{3.8in}}
-% \multicolumn2{c}{\textit{Nullary Operations}}\\
+% \multicolumn2{c}{\textit{Do Nothing}}\\
% \hline\hline
% \textbf{Command}&\textbf{operation}\\
% \hline
@@ -808,12 +820,12 @@
%
% \subsection{Errors}
%
-% If one tries to \op{pop}from an empty stack, an error message will be
+% If one tries to \op{pop} from an empty stack, an error message will be
% issued. Ignoring the error causes the macro to have the value stored
% in the macro \SpecialUsageIndex{\EndofStack}\verb$\EndofStack$.
% Its default is \texttt{0.00000000}.
%
-% If one tries to divide by $0$, an error message will be issued.
+% If one tries to divide by zero, an error message will be issued.
% Ignoring the error causes the result to be one of the following:
% \begin{itemize}
% \item Dividing $0$ by $0$ gives a result whose integer part is stored
@@ -859,10 +871,10 @@
\ifx \csname MFP@finish\endcsname\relax
\else \expandafter\endinput \fi
\expandafter\edef\csname MFP@finish\endcsname{%
- \catcode64=\the\catcode64 \space % @
- \catcode46=\the\catcode46 \space % .
- \catcode60=\the\catcode60 \space % <
- \catcode62=\the\catcode62 \space}% >
+ \catcode64=\the\catcode64 \space
+ \catcode46=\the\catcode46 \space
+ \catcode60=\the\catcode60 \space
+ \catcode62=\the\catcode62 \space}%
\ifx\ProvidesPackage\UndEfInEd
\newlinechar`\^^J%
\message{%
@@ -876,33 +888,27 @@
\DeclareOption{extra}{\def\MFPextra{}}%
\ProcessOptions\relax
\fi
-\catcode64=11 % @=letter (already is in LaTeX)
+\catcode64=11
\ifx\MFPextra\UndEfInEd
\def\MFP@loadextra{}%
\else
\def\MFP@loadextra{\input mfpextra\relax}%
\fi
-\def\MFPloadextra{%
- \edef\MFP@load@extra{%
- \catcode46=12 \catcode60=12 \catcode62=12 \catcode64=11
- \noexpand\input mfpextra\relax
- \catcode46=\the\catcode46\relax\catcode60=\the\catcode60\relax
- \catcode62=\the\catcode62\relax\catcode64=\the\catcode64\relax}%
- \MFP@load@extra}%
-\catcode46=12 % .
-\catcode60=12 % <
-\catcode62=12 % >
+\def\MFPloadextra{\input mfpextra\relax}%
+\catcode46=12
+\catcode60=12
+\catcode62=12
% \end{macrocode}
%
-% We check for \LaTeX{} (ignoring \LaTeX209); \cs{MFP@ifnoLaTeX}\dots\cs{mfp@end}
-% is skipped in LateX and executed otherwise.
+% We check for \LaTeX{} (ignoring \LaTeX209); \cs{MFP@ifnoLaTeX}\dots\cs{MFP@end}
+% is skipped in \LaTeX{} and executed otherwise.
% \begin{macrocode}
-\long\def\gobbleto@mfp@end#1\mfp@end{}%
-\ifx\mfp@end\UndEfInEd\def\mfp@end{\@empty}\fi
+\long\def\gobbleto@MFP@end#1\MFP@end{}%
+\def\MFP@end{\@empty}%
\ifx\documentclass\UndEfInEd
\def\MFP@ifnoLaTeX{}%
\else
- \let\MFP@ifnoLaTeX\gobbleto@mfp@end
+ \let\MFP@ifnoLaTeX\gobbleto@MFP@end
\fi
% \end{macrocode}
%
@@ -945,7 +951,7 @@
\long\def\@firstofone #1{#1}%
\long\def\@firstoftwo #1#2{#1}%
\long\def\@secondoftwo#1#2{#2}%
-\mfp@end
+\MFP@end
% \end{macrocode}
%
% We need to divide by both $10^4$ and $10^8$ several times. I could
@@ -958,15 +964,15 @@
%
% These are for manipulating digits. The \verb$\...ofmany$ commands
% require a sequence of arguments (brace groups or tokens) followed by
-% \verb$\mfp@end$. The minimum number of required parameters is surely
+% \verb$\MFP@end$. The minimum number of required parameters is surely
% obvious. For example, \cs{MFP@ninthofmany} must be used like\\
-% \indent\cs{MFP@ninthofmany}\meta{9 or more arguments}\cs{mfp@end}
+% \indent\cs{MFP@ninthofmany}\meta{9 or more arguments}\cs{MFP@end}\\
% All these are fully expandable.
% \begin{macrocode}
-\def\MFP@firstofmany#1#2\mfp@end{#1}%
-\def\MFP@fifthofmany#1#2#3#4#5#6\mfp@end{#5}%
-\def\MFP@ninthofmany#1#2#3#4#5#6#7#8{\MFP@firstofmany}%
-\def\MFP@firsteightofmany#1#2#3#4#5#6#7#8#9\mfp@end{#1#2#3#4#5#6#7#8}%
+\def\MFP@oneofmany#1#2\MFP@end{#1}%
+\def\MFP@fifthofmany#1#2#3#4#5#6\MFP@end{#5}%
+\def\MFP@ninthofmany#1#2#3#4#5#6#7#8{\MFP@oneofmany}%
+\def\MFP@eightofmany#1#2#3#4#5#6#7#8#9\MFP@end{#1#2#3#4#5#6#7#8}%
% \end{macrocode}
%
% \subsection{Processing numbers and the stack}
@@ -1011,7 +1017,8 @@
% \indent
% \cs{endgroup}\cs{def}\cs{MFP@z@Val}\marg{\meta{expansion-of-\cs{MFP@z@Val}}}\\
% which defines \cs{MFP@z@Val} outside the current group to equal its expansion
-% within the current group, provided it was started with \cs{begingroup}.
+% within the current group (provided the group was started with
+% \cs{begingroup}).
%
% We define a \cs{MFP@returned@values} to make all the conceivable produced
% values survive the group. The \cs{MFPcurr@Sgn} part is to permit testing
@@ -1020,30 +1027,27 @@
% I have been lax at making sure \cs{MFP@z@Ovr} is properly initiallized
% and properly checked whenever it could be relevant, and properly
% passed on. I think every internal command \cs{MFP@R}\textit{xxx}
-% should ensure it starts being $0$ and ends with a numerical value. I
-% notice that division might make it empty.
+% should ensure it starts being zero and ends with a numerical value. At
+% one time division could leave it undefined.
%
% \cs{MFP@subroutine} executes its argument (typically a single command) with
% a wrapper that initializes all the macros that might need initializing,
% and returns the necessary results.
% \begin{macrocode}
\def\MFP@endgroup@after#1{\edef\x{\endgroup#1}\x}%
-\def\MFP@endgroup@return{\MFP@endgroup@after\MFP@returned@values}%
-\def\MFP@def@after{\def\noexpand}%
+\def\MFP@afterdef{\def\noexpand}%
\def\MFP@returned@values{%
- \MFP@def@after\MFP@z@Val{\MFP@z@Sign\MFP@z@Int.\MFP@z@Frc}%
- \MFP@def@after\MFP@z@Ovr{\MFP@z@Ovr}%
- \MFP@def@after\MFP@z@Und{\MFP@z@Und}%
- \MFP@def@after\MFPcurr@Sgn{\MFP@z@Sgn}}%
+ \MFP@afterdef\MFP@z@Val{\MFP@z@Sign\MFP@z@Int.\MFP@z@Frc}%
+ \MFP@afterdef\MFP@z@Ovr{\MFP@z@Ovr}%
+ \MFP@afterdef\MFP@z@Und{\MFP@z@Und}%
+ \MFP@afterdef\MFPcurr@Sgn{\MFP@z@Sgn}}%
\def\MFP@subroutine#1{%
\begingroup
- \MFP@basic@init@z
+ \MFP@Rzero
+ \def\MFP@z@Ovr{0}%
+ \def\MFP@z@Und{0}%
#1%
- \MFP@endgroup@return}%
-\def\MFP@basic@init@z{%
- \MFP@Rzero
- \def\MFP@z@Ovr{0}%
- \def\MFP@z@Und{0}}%
+ \MFP@endgroup@after\MFP@returned@values}%
\def\MFP@Rzero{%
\def\MFP@z@Sgn{0}%
\def\MFP@z@Int{0}%
@@ -1052,7 +1056,7 @@
%
% \DescribeMacro{\EndofStack}
% We define here the error messages: popping from an empty stack and
-% dividing by $0$. In addition to the error messages, we provide some
+% dividing by zero. In addition to the error messages, we provide some
% default values that hopefully allow some operations to continue.
%
% We also have a warning or two.
@@ -1087,7 +1091,7 @@
% \DescribeMacro{\MaxRealInt}These are the largest possible integer and
% fractional parts of a real
% \DescribeMacro{\MaxRealFrac}number. They are returned for division by
-% $0$, for logarithm of $0$, and when overflow is detected in the
+% zero, for logarithm of zero, and when overflow is detected in the
% exponential function.
% \begin{macrocode}
\def\MaxRealInt {99999999}%
@@ -1096,11 +1100,11 @@
%
% \SpecialUsageIndex{\MaxRealInt}
% \SpecialUsageIndex{\MaxRealFrac}
-% These are the results returned when trying to divide by $0$. Two are
+% These are the results returned when trying to divide by zero. Two are
% \DescribeMacro{\xOverZeroInt}
% \DescribeMacro{\xOverZeroFrac}
-% used when dividing a nonzero number by $0$ and and two when trying to
-% divide $0$ by $0$.
+% used when dividing a nonzero number by zero and and two when trying to
+% divide zero by zero.
% \DescribeMacro{\ZeroOverZeroInt}
% \DescribeMacro{\ZeroOverZeroFrac}
% \begin{macrocode}
@@ -1130,7 +1134,7 @@
\def\MFPparse@real#1#2#3#4{%
\MFPnospace@def\MFPtemp@Val{#4}%
\MFPprocess@into@parts\MFPtemp@Val#1#2#3%
- \MFPpadto@eight#3}%
+ \MFP@padtoeight#3}%
\def\MFPparse@x{\MFPparse@real\MFP@x@Sgn\MFP@x@Int\MFP@x@Frc}%
\def\MFPparse@y{\MFPparse@real\MFP@y@Sgn\MFP@y@Int\MFP@y@Frc}%
% \end{macrocode}
@@ -1149,29 +1153,27 @@
% the sign, integer and fractional parts.
% \begin{macrocode}
\def\MFPprocess@into@parts#1#2#3#4{%
- \@xp\MFPsplit@dot#1..\mfp@end #3#4%
+ \@xp\MFPsplit@dot#1..\MFP@end #3#4%
% \end{macrocode}
%
-% This is the first place where having at most eight digits simplifies things.
-% At this point \arg3 could contain any number of consecutive signs
-% followed by any eight digits. It could be $0$, so to avoid losing the sign
-% we append a \texttt{1} (for up to nine digits). We temporarily define the
-% sign based on the result, but may need to drop it if both the integer
-% and fractional parts are $0$.
+% At this point \arg3 holds the part before the dot (or the whole thing
+% if there was no dot) and \arg4 holds the part after the dot, (or
+% nothing). Now is the first place where having at most eight digits
+% simplifies things. Note that \arg3 could contain any number of
+% consecutive signs followed by up to eight digits. It could be zero or
+% empty, so to avoid losing the sign we append a \texttt{1} (for up to
+% nine digits). We temporarily define the sign based on the result, but
+% may need to drop it if both the integer and fractional parts are zero.
%
-% Prepending a 0 to the fractional part permits it to be empty.
+% Prepending a zero to the fractional part pemits it to be empty.
% In the final \cs{edef}, \arg3 is made positive.
% \begin{macrocode}
- \ifnum#31<0
- \def#2{-1}%
- \else
- \def#2{1}%
+ \ifnum#31<0 \def#2{-1}%
+ \else \def#2{1}%
\fi
\ifnum #30=0
\def#3{0}%
- \ifnum 0#4=0
- \def#2{0}%
- \fi
+ \ifnum 0#4=0 \def#2{0}\fi
\fi
\edef#3{\number \ifnum #2<0 -\fi#3}}%
% \end{macrocode}
@@ -1179,15 +1181,15 @@
% This only copies the parts before and after the dot, \arg1 and \arg2,
% into macros \arg4 and \arg5.
% \begin{macrocode}
-\def\MFPsplit@dot#1.#2.#3\mfp@end#4#5{\edef#4{#1}\edef#5{#2}}%
+\def\MFPsplit@dot#1.#2.#3\MFP@end#4#5{\edef#4{#1}\edef#5{#2}}%
% \end{macrocode}
%
% This is used to pad the fractional part to eight places with zeros. If
% a number with more than eight digits survives to this point, it gets
% truncated.
% \begin{macrocode}
-\def\MFPpadto@eight#1{%
- \edef#1{\@xp\MFP@firsteightofmany#100000000\mfp@end}}%
+\def\MFP@padtoeight#1{%
+ \edef#1{\@xp\MFP@eightofmany#100000000\MFP@end}}%
% \end{macrocode}
%
% These take operands off the stack. We know already that there are no
@@ -1215,6 +1217,7 @@
% \end{macrocode}
%
% Sometimes only parts of the number needs changing (used in CHS, ABS).
+% This copies the integer and fractional parts of $x$ into $z$.
% \begin{macrocode}
\def\copyMFP@x{\edef\MFP@z@Int{\MFP@x@Int}\edef\MFP@z@Frc{\MFP@x@Frc}}%
% \end{macrocode}
@@ -1238,7 +1241,7 @@
% \end{macrocode}
%
% The macro \cs{Rpop} calls \cs{MFP@popit} followed by the contents of the
-% stack, the token \cs{mfp@end} and the macro to \op{pop} into. If the stack is
+% stack, the token \cs{MFP@end} and the macro to \op{pop} into. If the stack is
% not empty, \cs{doMFP@popit} will read the first group \arg1 into that macro
% \arg3, and then redefine the stack to be the rest of the argument \arg2.
% If the stack is empty, \cs{doMFP@EOS} will equate the macro to
@@ -1246,8 +1249,8 @@
% message.
% \begin{macrocode}
\def\MFP@popit{\if@EndofStack\doMFP@EOS\doMFP@popit}%
-\def\doMFP@EOS#1\mfp@end#2{\MFP@popempty@err\let#2\EndofStack}%
-\def\doMFP@popit#1#2\mfp@end#3{\edef\MFP@Rstack{#2}\edef#3{#1}}%
+\def\doMFP@EOS#1\MFP@end#2{\MFP@popempty@err\let#2\EndofStack}%
+\def\doMFP@popit#1#2\MFP@end#3{\edef\MFP@Rstack{#2}\edef#3{#1}}%
% \end{macrocode}
%
% \subsection{The user-level operations}
@@ -1365,7 +1368,7 @@
% \item[cmp] compare $x$ and $y$ (stack version does not change stack).
% \item[chk] examine the sign of $x$ (stack version does not change stack).
% \item[dup] stack only, duplicate the top element of the stack.
-% \item[push] stack only, put a value on top of the stack.
+% \item[push] stack only, put a value onto the top of the stack.
% \item[pop] stack only, remove the top element of the stack,
% store it in a variable.
% \item[exch] stack only, exchange top two elements of the stack.
@@ -1449,8 +1452,8 @@
% \begin{macrocode}
\let\Rnoop\relax
\def\Rcmp{%
- \MFPgetoperand@y\MFPgetoperand@x % get operands (last pushed is y)
- \MFP@Rcat\MFP@x@Val\MFP@Rcat\MFP@y@Val % put back: LOFI
+ \MFPgetoperand@y\MFPgetoperand@x
+ \MFP@Rcat\MFP@x@Val\MFP@Rcat\MFP@y@Val
\MFP@Rcmp}%
\def\Rchk{%
\MFPgetoperand@x
@@ -1461,7 +1464,7 @@
\edef\MFP@z@Val{\MFP@x@Sign\MFP@x@Int.\MFP@x@Frc}%
\edef\MFPcurr@Sgn{\MFP@x@Sgn}%
\MFPpush@result}%
- \def\Rpop{\@xp\MFP@popit\MFP@Rstack\mfp@end}%
+ \def\Rpop{\@xp\MFP@popit\MFP@Rstack\MFP@end}%
\def\Rexch{%
\Rpop\MFP@x@Val\Rpop\MFP@y@Val
\MFP@Rcattwo\MFP@y@Val\MFP@x@Val}%
@@ -1546,6 +1549,17 @@
\def\MFPmax{\MFP@op@Binary\MFP@Rmax}%
% \end{macrocode}
%
+% A \emph{nullary} operation is one that produces a result with no
+% operand. Thus, it could return a fixed constant, or it could perform
+% calculations that obtain input from the system (e.g., current time). At
+% the moment we don't define any.
+% \begin{macrocode}
+\def\MFP@stack@Nullary#1{%
+ \MFP@subroutine{#1}\MFPpush@result}%
+\def\MFP@op@Nullary#1{%
+ \MFP@subroutine{#1}\MFPstore@result}%
+% \end{macrocode}
+%
% These are the wrappers for unary operations. The operand versions have a
% second argument, the macro that stores the result. But this will be the
% argument of \cs{MFPstore@result}.
@@ -1604,7 +1618,7 @@
\MFP@tempb\MFP@x@Frc\relax
\ifodd\MFP@tempb
\def\MFP@z@Und{5}%
- \advance\MFP@tempb 1 % round up
+ \advance\MFP@tempb 1
\ifnum\MFP@ttteight=\MFP@tempb
\MFP@tempb0 \advance\MFP@tempa1
\fi
@@ -1623,13 +1637,15 @@
% \end{macrocode}
%
% The squaring operation just calls \cs{MFP@Rmul} after copying $x$ to
-% $y$.
+% $y$. Its gain in efficiency over a multiplication is that it can skip
+% preprocessing of the second (identical) operand.
% \begin{macrocode}
\def\MFP@Rsq{\MFP@Rcopy xy\MFP@Rmul}%
% \end{macrocode}
%
% The inversion operation just calls \cs{MFP@Rdiv} after copying $x$ to
-% $y$ and $1$ to $x$.
+% $y$ and $1$ to $x$. Its advantage over a divide is it skips the
+% preprocessing of $1$ as an operand.
% \begin{macrocode}
\def\MFP@Rinv{\MFP@Rcopy xy\MFP@Rload x110\MFP@Rdiv}%
% \end{macrocode}
@@ -1640,7 +1656,7 @@
\MFP@Rloadz {\ifnum\MFP@x@Int=0 0\else\MFP@x@Sgn\fi}\MFP@x@Int 0}%
% \end{macrocode}
%
-% Fractional part: replace integer part with a $0$.
+% Fractional part: replace integer part with a zero.
% \begin{macrocode}
\def\MFP@Rfrac{%
\MFP@Rloadz {\ifnum\MFP@x@Frc=0 0\else\MFP@x@Sgn\fi}0\MFP@x@Frc}%
@@ -1649,8 +1665,8 @@
% To increment and decrement by $1$, except in border cases, we need only
% address the integer part of a number. This doesn't seem so simple
% written out but, even so, it is more efficient than full-blown addition.
-% It would be very slightly more efficient to repeat the increment code in
-% decrementing, but it would be annoying to do so,
+% It would be very slightly more efficient if \cs{MFP@Rdecr} did not call
+% \cs{MFP@Rincr}, but instead was similarly coded.
% \begin{macrocode}
\def\MFP@Rincr{%
\ifnum\MFP@x@Sgn<0
@@ -1685,17 +1701,17 @@
% We use the same code to get floor or ceiling, the
% appropriate inequality character being its argument.
% \begin{macrocode}
-\def\MFP@Rfloororceil#1{%
+\def\MFP@Rfloorceil#1{%
\MFP@tempa\MFP@x@Int\relax
- \ifnum 0#1\MFP@x@Sgn
+ \ifnum \MFP@x@Sgn #10
\ifnum\MFP@x@Frc=0
\else
\advance\MFP@tempa1
\fi
\fi
\MFP@Rloadz{\ifnum\MFP@x@Int=0 0\else\MFP@x@Sgn\fi}\MFP@tempa0}%
-\def\MFP@Rfloor{\MFP@Rfloororceil>}%
-\def\MFP@Rceil {\MFP@Rfloororceil<}%
+\def\MFP@Rfloor{\MFP@Rfloorceil<}%
+\def\MFP@Rceil {\MFP@Rfloorceil>}%
% \end{macrocode}
%
% For multiplication, after the usual break into integer and fractional
@@ -1713,14 +1729,13 @@
\multiply\MFP@tempb by\MFP@tttfour
\advance\MFP@tempa-\MFP@tempb
\MFP@endgroup@after{%
- \MFP@def@after#2{#2}%
- \MFP@def@after#3{\number\MFP@tempa}%
+ \MFP@afterdef#2{#2}%
+ \MFP@afterdef#3{\number\MFP@tempa}%
}}%
%
-\def\MFP@x@split{%
+\def\MFP@@split{%
\MFP@split\MFP@x@Int\MFP@x@Int@ii\MFP@x@Int@i
- \MFP@split\MFP@x@Frc\MFP@x@Frc@i\MFP@x@Frc@ii}%
-\def\MFP@y@split{%
+ \MFP@split\MFP@x@Frc\MFP@x@Frc@i\MFP@x@Frc@ii
\MFP@split\MFP@y@Int\MFP@y@Int@ii\MFP@y@Int@i
\MFP@split\MFP@y@Frc\MFP@y@Frc@i\MFP@y@Frc@ii}%
% \end{macrocode}
@@ -1728,39 +1743,27 @@
% We will store the intermediate and final products in \cs{MFP@z@*}. Each one
% is ultimately reduced to four digits, like the parts of $x$ and $y$. As each
% base-$10000$ digit of $y$ is multiplied by a digit of $x$, we add the
-% result to the appropriate digit of the partial result $z$. Thus, we need
-% to zero out $z$ at the start (or treat the first iteration differently):
+% result to the appropriate digit of the partial result $z$.
%
% The underflow ends up in \cs{MFP@z@Frc@iv} and \cs{MFP@z@Frc@iii}.
% Overflow will be in \cs{MFP@z@Int@iii}. Unlike the rest, it can be up to
% eight digits because we do not need to carry results out of it.
-% \begin{macrocode}
-\def\MFPmore@init@z{%
- \def\MFP@z@Frc@iv {0}%
- \def\MFP@z@Frc@iii{0}%
- \def\MFP@z@Frc@ii {0}%
- \def\MFP@z@Frc@i {0}%
- \def\MFP@z@Int@i {0}%
- \def\MFP@z@Int@ii {0}%
- \def\MFP@z@Int@iii{0}}%
-% \end{macrocode}
%
-% This command prepends zeros so a number fills four slots. In the
-% ``make'' version, \arg1 is a macro holding the value and is redefined to
-% contain the result. A macro that calls these should ensure that \arg1 is
-% not empty and is less than 10,000.
+% This command prepends zeros so a number fills four slots. Here \arg1 is
+% a macro holding the value and it is redefined to contain the result. A
+% macro that calls this should ensure that \arg1 is not empty and is less
+% than 10,000.
% \begin{macrocode}
-\def\MFP@fourdigits#1{%
- \@xp\MFP@fifthofmany\number#1{}{0}{00}{000}\mfp@end\number#1}%
-\def\makeMFP@fourdigits#1{\edef#1{\MFP@fourdigits{#1}}}%
+\def\makeMFP@fourdigits#1{%
+ \edef#1{\@xp\MFP@fifthofmany\number#1{}{0}{00}{000}\MFP@end\number#1}}%
% \end{macrocode}
%
-% This is the same, but produce eight digits.
+% This is the same, but produces eight digits. Similarly \arg1 should be
+% nonempty and less than 100,000,000.
% \begin{macrocode}
-\def\MFP@eightdigits#1{%
- \@xp\MFP@ninthofmany\number#1%
- {}{0}{00}{000}{0000}{00000}{000000}{0000000}\mfp@end\number#1}%
-\def\makeMFP@eightdigits#1{\edef#1{\MFP@eightdigits{#1}}}%
+\def\makeMFP@eightdigits#1{%
+ \edef#1{\@xp\MFP@ninthofmany\number#1%
+ {}{0}{00}{000}{0000}{00000}{000000}{0000000}\MFP@end\number#1}}%
% \end{macrocode}
%
% The following macros implement carrying. The macros \cs{MFP@carrya} and
@@ -1779,22 +1782,19 @@
\begingroup
\MFP@carryi{#1}#2#3%
\MFP@endgroup@after{%
- \MFP@def@after#3{\number\MFP@tempa}%
- \MFP@def@after#2{\number\MFP@tempb}%
+ \MFP@afterdef#3{\number\MFP@tempa}%
+ \MFP@afterdef#2{\number\MFP@tempb}%
}}%
% \end{macrocode}
%
% This is the ``internal'' carry. \arg1, \arg2, and \arg3 are as in
-% \cs{MFP@carry}. Its advantage is that it can be used used where \arg2 and
-% \arg3 are not macros, leaving the result in \cs{MFP@tempa} and \cs{MFP@tempb}
-% with \cs{MFP@tempb} in the correct range, $[0,\mbox{\arg1})$. Its
-% disadvantage is it does not protect temporary registers. Warning:
-% never use it in the form \cs{MFP@carryi}\meta{num}\cs{MFP@tempa},
-% because this would copy \meta{num} to \cs{MFP@tempa}, losing the value
-% in the second argument before anything can be done. The other order is
-% okay, and \cs{MFP@tempb} can be used in either slot. Do not use it
-% without grouping if you want the values in the temp registers \texttt{a},
-% \texttt{b} or \texttt{c} preserved.
+% \cs{MFP@carry}. Its advantage is that it can be used used where \arg2
+% and \arg3 are not macros, leaving the result in \cs{MFP@tempa} and
+% \cs{MFP@tempb} with \cs{MFP@tempb} in the correct range,
+% $[0,\mbox{\arg1})$. Its disadvantage is it does not protect temporary
+% registers. Warning: do not use it with \arg2=\cs{MFP@tempa} and do not
+% use it without grouping if you want to preserve the values in these
+% temporary count registers.
% \begin{macrocode}
\def\MFP@carryi#1#2#3{%
\MFP@tempa=#3\relax
@@ -1814,7 +1814,7 @@
\MFP@tempa#1%
\advance\MFP@tempa#2\relax
\MFP@endgroup@after{%
- \MFP@def@after#3{\number\MFP@tempa}%
+ \MFP@afterdef#3{\number\MFP@tempa}%
}}%
% \end{macrocode}
%
@@ -1940,43 +1940,44 @@
% \end{macrocode}
%
% \cs{MFP@Rmul} first computes the (theoretical) sign of the product: if
-% $0$, return $0$, otherwise provisionally set the sign of the product and
-% call \cs{MFP@@Rmul}.
+% it is zero, return zero, otherwise provisionally set the sign of the product
+% and call \cs{MFP@@Rmul}.
% \begin{macrocode}
\def\MFP@Rmul{%
- \MFP@tempa\MFP@x@Sgn \multiply\MFP@tempa\MFP@y@Sgn\relax
- \ifnum 0=\MFP@tempa
- \MFP@Rzero
- \else
- \edef\MFP@z@Sgn{\number\MFP@tempa}%
- \@xp\MFP@@Rmul
- \fi}%
+ \ifnum\MFP@x@Sgn=0 \MFP@Rzero
+ \else\ifnum\MFP@y@Sgn=0 \MFP@Rzero
+ \else \edef\MFP@z@Sgn{\number\MFP@x@Sign\MFP@y@Sgn}%
+ \@XP\MFP@@Rmul
+ \fi\fi}%
% \end{macrocode}
%
-% \cs{MFP@@Rmul} splits the four expected macros into eight macros
-% considered to be four base-10000 digits for each of $x$ and $y$.
+% \cs{MFP@@Rmul} first initializes the macros that will hold the
+% base-10000 digits of $z$. Then it splits the four expected macros into
+% eight macros that hold the base-10000 digits for each of $x$ and $y$.
% Then each digit of $y$ is used to multiply the four digits of $x$ and the
-% results are added to corresponding digits of $z$, which have been
-% initialized to $0$ by \cs{MFPmore@init@z}.
+% results are added to corresponding digits of $z$.
% \begin{macrocode}
\def\MFP@@Rmul{%
- \MFPmore@init@z
- \MFP@x@split\MFP@y@split
+ \def\MFP@z@Frc@iv {0}\def\MFP@z@Frc@iii{0}%
+ \def\MFP@z@Frc@ii {0}\def\MFP@z@Frc@i {0}%
+ \def\MFP@z@Int@i {0}\def\MFP@z@Int@ii {0}%
+ \def\MFP@z@Int@iii{0}%
+ \MFP@@split
\MFP@multiplyfour \MFP@y@Frc@ii \MFP@z@Frc@i
- \MFP@z@Frc@ii \MFP@z@Frc@iii\MFP@z@Frc@iv
+ \MFP@z@Frc@ii \MFP@z@Frc@iii \MFP@z@Frc@iv
\MFP@multiplyfour \MFP@y@Frc@i \MFP@z@Int@i
- \MFP@z@Frc@i \MFP@z@Frc@ii \MFP@z@Frc@iii
+ \MFP@z@Frc@i \MFP@z@Frc@ii \MFP@z@Frc@iii
\MFP@multiplyfour \MFP@y@Int@i \MFP@z@Int@ii
- \MFP@z@Int@i \MFP@z@Frc@i \MFP@z@Frc@ii
+ \MFP@z@Int@i \MFP@z@Frc@i \MFP@z@Frc@ii
\MFP@multiplyfour \MFP@y@Int@ii \MFP@z@Int@iii
- \MFP@z@Int@ii \MFP@z@Int@i \MFP@z@Frc@i
+ \MFP@z@Int@ii \MFP@z@Int@i \MFP@z@Frc@i
% \end{macrocode}
% Now apply the carry routines on the underflow digits\dots
% \begin{macrocode}
\MFP@carrym\MFP@z@Frc@iv\MFP@z@Frc@iii
\MFP@carrym\MFP@z@Frc@iii\MFP@z@Frc@ii
% \end{macrocode}
-% \dots pause to round the lowest digit that will be kept\dots
+% \dots and pause to round the lowest digit that will be kept\dots
% \begin{macrocode}
\ifnum\MFP@z@Frc@iii<5000 \else
\MFP@tempb\MFP@z@Frc@ii
@@ -1992,13 +1993,13 @@
\MFP@carrym\MFP@z@Int@ii\MFP@z@Int@iii
% \end{macrocode}
% To end, we arrange for all macros to hold four digits (except
-% \cs{MFP@z@Int@ii} which doesn't need leading 0s, and \cs{MFP@z@Int@iii}
-% which also doesn't) and load them into the appropriate 8-digit macros.
-% The underflow digits are stored in \cs{MFP@z@Und} in case we ever need
-% to examine them, and the overflow in \cs{MFP@z@Ovr} in case we ever need
-% to implement an overflow error. Theoretically $z \ne 0$, but it is
-% possible that $z=0$ after rounding to eight places. If so, we must reset
-% \cs{MFP@z@Sgn}.
+% \cs{MFP@z@Int@ii} and \cs{MFP@z@Int@iii} which don't need leading 0s)
+% and load them into the appropriate 8-digit macros. The underflow digits
+% are stored in \cs{MFP@z@Und} in case we ever need to examine them (we
+% now do: in our unit conversion routine \cs{MFP@DPmul}), and the overflow
+% in \cs{MFP@z@Ovr} in case we ever want to implement an overflow error.
+% Theoretically $z \ne 0$, but it is possible that $z=0$ after reducing to
+% eight places. If so, we must reset \cs{MFP@z@Sgn}.
% \begin{macrocode}
\makeMFP@fourdigits\MFP@z@Frc@iv
\makeMFP@fourdigits\MFP@z@Frc@iii
@@ -2011,8 +2012,7 @@
\edef\MFP@z@Und{\MFP@z@Frc@iii\MFP@z@Frc@iv}%
\ifnum\MFP@z@Int>0
\else\ifnum\MFP@z@Frc>0
- \else
- \def\MFP@z@Sgn{0}%
+ \else \def\MFP@z@Sgn{0}%
\fi\fi}%
% \end{macrocode}
%
@@ -2031,12 +2031,12 @@
% only $16$ significant digits should be retained in any case.) If $d$ is
% $0$ and $n$ is $15$ we would need $-5$ digits. That means the first
% nonzero digit is in the 15th or 16th place after the dot and the
-% quotient is effectively $0$.
+% quotient is effectively zero.
%
% Here I explain why we normalize the parts in this way. If a numerator
% has the form $n_1.n_2$ and the denominator has the form $d_1.d_2$ then
-% TeX can easily obtain the integer part of $n_1/d_1$, because these are
-% within its range for integers. The resulting quotient (let's call it
+% \TeX{} can easily obtain the integer part of $n_1/d_1$, because these
+% are within its range for integers. The resulting quotient (let's call it
% $q_1$) is the largest integer satisfying $q_1d_1 \le n_1$. What we seek,
% however is the largest integer $q$ such that $q(d_1.d_2) \le n_1.n_2$.
% It can easily be shown that $q \le q_1$. It is true, but not so easily
@@ -2058,7 +2058,7 @@
% Since $d_2$ is no more than eight digits, $q_1 d_2$ is less than $9
% (10)^8$. Inequality (\ref{crucial}) is therefore satisfied if $n_1 - q_1
% d_1 \ge 9$. If that is not the case then the right side of
-% (\ref{crucial}) is computable within TeX's integer ranges and we can
+% (\ref{crucial}) is computable within \TeX's integer ranges and we can
% easily test the inequality. If the inequality holds, then $q = q_1$,
% otherwise $q = q_1 - 1$.
%
@@ -2069,19 +2069,18 @@
%
% Now I need to get it organized. \cs{MFP@Rdiv} will have \cs{MFP@x@*} and
% \cs{MFP@y@*} available. One step (could be first or last). Is to calculate
-% the sign. Let's do it first (because we need to check for $0$ anyway).
+% the sign. Let's do it first (because we need to check for zero anyway).
%
-% We invoke an error message upon division by $0$, but nevertheless return
+% We invoke an error message upon division by zero, but nevertheless return
% a value. By default it is $0$ for $0/0$ and the maximum possible real
-% for $x/0$ when $x$ is not $0$. If the numerator is $0$ and the
-% denominator not, we do nothing as $z$ was initialized to be $0$.
+% for $x/0$ when $x$ is not zero. If the numerator is zero and the
+% denominator not, we do nothing as $z$ was initialized to be zero.
%
-% If neither is $0$, we calculate the sign of the result and call
+% If neither is zero, we calculate the sign of the result and call
% \cs{MFP@@Rdiv} to divide the absolute values.
% \begin{macrocode}
\def\MFP@Rdiv{%
- \ifnum\MFP@y@Sgn=0
- \MFP@dividebyzero@err
+ \ifnum\MFP@y@Sgn=0 \MFP@dividebyzero@err
\ifnum\MFP@x@Sgn=0
\edef\MFP@z@Int{\ZeroOverZeroInt}%
\edef\MFP@z@Frc{\ZeroOverZeroFrac}%
@@ -2090,13 +2089,8 @@
\edef\MFP@z@Frc{\xOverZeroFrac}%
\fi
\edef\MFP@z@Sgn{\MFP@x@Sgn}%
- \else\ifnum\MFP@x@Sgn=0
- \MFP@Rzero
- \else
- \MFP@tempa\MFP@x@Sgn
- \multiply\MFP@tempa\MFP@y@Sgn
- \edef\MFP@z@Sgn{\number\MFP@tempa}%
- \MFP@@Rdiv
+ \else\ifnum\MFP@x@Sgn=0 \MFP@Rzero
+ \else \edef\MFP@z@Sgn{\number\MFP@x@Sign\MFP@y@Sgn}\MFP@@Rdiv
\fi\fi}%
% \end{macrocode}
%
@@ -2106,7 +2100,7 @@
% but knowing the shift will give us the correct quotient in the end.
%
% We first arrange that \cs{MFP@y@Int} is nonzero by making it \cs{MFP@y@Frc} if
-% it is $0$ (a shift of eight digits). Then the macro
+% it is zero (a shift of eight digits). Then the macro
% \cs{MFP@numdigits@toshift} computes $8$ minus the number of digits in
% \cs{MFP@y@Int}, which is how many positions left $y$ will be shifted.
% We then call \cs{MFP@doshift@y} on the concatenation of the digits in
@@ -2123,7 +2117,7 @@
\MFP@tempa=0
\fi
\advance\MFP@tempa\MFP@numdigits@toshift\MFP@y@Int\relax
- \@XP\MFP@doshift@y\@xp\MFP@y@Int\MFP@y@Frc0000000\mfp@end
+ \@XP\MFP@doshift@y\@xp\MFP@y@Int\MFP@y@Frc0000000\MFP@end
% \end{macrocode}
%
% We repeat all that on the numerator $x$, except shifting its digits
@@ -2142,7 +2136,7 @@
\advance\MFP@tempa -8
\fi
\advance\MFP@tempa-\MFP@numdigits@toshift\MFP@x@Int\relax
- \@XP\MFP@doshift@x\@xp\MFP@x@Int\MFP@x@Frc0000000\mfp@end
+ \@XP\MFP@doshift@x\@xp\MFP@x@Int\MFP@x@Frc0000000\MFP@end
% \end{macrocode}
%
% Since our result will have at most one digit in the integer part, a
@@ -2184,13 +2178,13 @@
% from the rest (\arg9) inside \cs{MFP@x@Frc}. The same with
% \cs{MFP@doshift@y}.
% \begin{macrocode}
-\def\MFP@numdigits@toshift#1{\@xp\MFP@ninthofmany#101234567\mfp@end}%
-\def\MFP@doshift@x#1#2#3#4#5#6#7#8#9\mfp@end{%
+\def\MFP@numdigits@toshift#1{\@xp\MFP@ninthofmany#101234567\MFP@end}%
+\def\MFP@doshift@x#1#2#3#4#5#6#7#8#9\MFP@end{%
\def\MFP@x@Int{#1#2#3#4#5#6#7#8}%
- \edef\MFP@x@Frc{\MFP@firsteightofmany#9\mfp@end}}%
-\def\MFP@doshift@y#1#2#3#4#5#6#7#8#9\mfp@end{%
+ \edef\MFP@x@Frc{\MFP@eightofmany#9\MFP@end}}%
+\def\MFP@doshift@y#1#2#3#4#5#6#7#8#9\MFP@end{%
\def \MFP@y@Int{#1#2#3#4#5#6#7#8}%
- \edef\MFP@y@Frc{\MFP@firsteightofmany#9\mfp@end}}%
+ \edef\MFP@y@Frc{\MFP@eightofmany#9\MFP@end}}%
% \end{macrocode}
%
% The loop counter is \cs{MFP@tempf}, \cs{MFP@tempa} is reserved for the
@@ -2252,7 +2246,7 @@
\fi
\advance\MFP@tempf -1
\ifnum\MFP@tempf>0
- \edef\MFP@x@Int{\MFP@x@Int0}% easy multiplications by 10
+ \edef\MFP@x@Int{\MFP@x@Int0}%
\edef\MFP@x@Frc{\MFP@x@Frc0}%
\MFP@carrya\MFP@x@Frc\MFP@x@Int
\@xp\MFP@Rdivloop
@@ -2289,7 +2283,7 @@
\advance \MFP@tempa -7
\ifnum\MFP@tempa>0
\def\MFP@z@Ovr{}%
- \@xp\MFPget@Ovrdigits\MFP@z@digits\mfp@end
+ \@xp\MFPget@Ovrdigits\MFP@z@digits\MFP@end
\else
\ifnum\MFP@tempa<-7
\edef\MFP@z@digits{00000000\MFP@z@digits}%
@@ -2307,7 +2301,7 @@
0000\else
00000%
\fi \MFP@z@digits}%
- \@xp\MFPget@Intdigits\MFP@z@digits\mfp@end
+ \@xp\MFPget@Intdigits\MFP@z@digits\MFP@end
\fi}%
% \end{macrocode}
%
@@ -2338,7 +2332,7 @@
\MFPget@Frcdigits}%
\def\MFPget@Frcdigits#1#2#3#4#5#6#7#8#9{%
\def\MFP@z@Frc{#1#2#3#4#5#6#7#8}%
- \def\MFP@z@Und{#9}\gobbleto@mfp@end}%
+ \def\MFP@z@Und{#9}\gobbleto@MFP@end}%
% \end{macrocode}
%
% The max amd min operations simply run the compare operation and use
@@ -2358,7 +2352,7 @@
% keeps the right number. For negative truncations we prepend zeros to the
% integer part so it too is exactly eight digits. These become the
% arguments of \cs{MFP@@iRtrunc}, which substitutes 0 for the last
-% \texttt{-\cs{MFP@tempa}} of them.
+% \texttt{-}\cs{MFP@tempa} of them.
%
% The macro to store the result in follows \arg2. It is read and
% defined by either \cs{MFP@Rtrunc} or \cs{MFP@iRtrunc}.
@@ -2373,15 +2367,15 @@
\@xp\MFP@Rtrunc
\fi}%
\def\MFP@Rtrunc#1{%
- \edef\MFP@x@Frc{\@xp\MFP@@Rtrunc\MFP@x@Frc\mfp@end}%
- \ifnum\MFP@x@Int=0 % possibly returns 0
+ \edef\MFP@x@Frc{\@xp\MFP@@Rtrunc\MFP@x@Frc\MFP@end}%
+ \ifnum\MFP@x@Int=0
\ifnum\MFP@x@Frc=0
\def\MFP@x@Sgn{0}%
\fi
\fi
\MFP@endgroup@after{%
- \MFP@def@after#1{\MFP@x@Sign\MFP@x@Int.\MFP@x@Frc}}}%
-\def\MFP@@Rtrunc#1#2#3#4#5#6#7#8#9\mfp@end{%
+ \MFP@afterdef#1{\MFP@x@Sign\MFP@x@Int.\MFP@x@Frc}}}%
+\def\MFP@@Rtrunc#1#2#3#4#5#6#7#8#9\MFP@end{%
\ifcase\MFP@tempa\or
#1\or
#1#2\or
@@ -2393,9 +2387,9 @@
#1#2#3#4#5#6#7#8\fi}%
\def\MFP@iRtrunc#1{%
\makeMFP@eightdigits\MFP@x@Int
- \edef\MFP@x@Val{\number\MFP@x@Sign\@xp\MFP@@iRtrunc\MFP@x@Int\mfp@end}%
- \MFP@endgroup@after{\MFP@def@after#1{\MFP@x@Val}}}%
-\def\MFP@@iRtrunc#1#2#3#4#5#6#7#8#9\mfp@end{%
+ \edef\MFP@x@Val{\number\MFP@x@Sign\@xp\MFP@@iRtrunc\MFP@x@Int\MFP@end}%
+ \MFP@endgroup@after{\MFP@afterdef#1{\MFP@x@Val}}}%
+\def\MFP@@iRtrunc#1#2#3#4#5#6#7#8#9\MFP@end{%
\ifcase-\MFP@tempa
#1#2#3#4#5#6#7#8\or
#1#2#3#4#5#6#70\or
@@ -2419,41 +2413,41 @@
\ifnum 0>\MFP@tempa
\edef\MFP@y@Tmp{%
\ifcase-\MFP@tempa\or
- 5\or % .5 x 10^1
- 50\or % .5 x 10^2
- 500\or % .5 x 10^3
- 5000\or % .5 x 10^4
- 50000\or % .5 x 10^5
- 500000\or % .5 x 10^6
- 5000000\else % .5 x 10^7
- 50000000\fi % .5 x 10^8
+ 5\or
+ 50\or
+ 500\or
+ 5000\or
+ 50000\or
+ 500000\or
+ 5000000\else
+ 50000000\fi
}%
\else
\edef\MFP@y@Tmp{%
\ifcase\MFP@tempa
- .5\or % .5 x 10^0
- .05\or % .5 x 10^{-1}
- .005\or % .5 x 10^{-2}
- .0005\or % .5 x 10^{-3}
- .00005\or % .5 x 10^{-4}
- .000005\or % .5 x 10^{-5}
- .0000005\or % .5 x 10^{-6}
- .00000005\else% .5 x 10^{-7}
- 0\fi %
+ .5\or
+ .05\or
+ .005\or
+ .0005\or
+ .00005\or
+ .000005\or
+ .0000005\or
+ .00000005\else
+ 0\fi
}%
\fi
\MFPchk{#2}\ifMFP@neg\edef\MFP@y@Tmp{-\MFP@y@Tmp}\fi
\MFPadd{#2}\MFP@y@Tmp\MFP@z@Tmp
- \MFP@endgroup@after{\MFP@def@after\MFP@z@Tmp{\MFP@z@Tmp}}%
+ \MFP@endgroup@after{\MFP@afterdef\MFP@z@Tmp{\MFP@z@Tmp}}%
\MFPtruncate{#1}\MFP@z@Tmp}%
% \end{macrocode}
%
% \DescribeMacro{\MFPstrip}
% Stripping zeros from the right end of the fractional part. The star form
-% differs only in the handling of a $0$ fractional part. So we check
-% whether it is $0$ and when it is, we either append `\texttt{.0}' or
+% differs only in the handling of a zero fractional part. So we check
+% whether it is zero and when it is, we either append `\texttt{.0}' or
% nothing. The rest of the code grabs a digit at a time and stops when the
-% rest are $0$.
+% rest are zero.
% \begin{macrocode}
\def\MFPstrip{%
\@ifstar{\MFP@strip{}}{\MFP@strip{.0}}}%
@@ -2462,15 +2456,15 @@
\ifnum \MFP@x@Frc=0
\edef#3{\MFP@x@Sign\MFP@x@Int#1}%
\else
- \edef#3{\MFP@x@Sign\MFP@x@Int.\@xp\MFP@@strip\MFP@x@Frc\mfp@end}%
+ \edef#3{\MFP@x@Sign\MFP@x@Int.\@xp\MFP@@strip\MFP@x@Frc\MFP@end}%
\fi}%
-\def\MFP@@strip#1#2\mfp@end{%
+\def\MFP@@strip#1#2\MFP@end{%
#1%
\ifnum 0#2>0
\@xp\MFP@@strip
\else
- \@xp\gobbleto@mfp@end
- \fi#2\mfp@end}%
+ \@xp\gobbleto@MFP@end
+ \fi#2\MFP@end}%
% \end{macrocode}
%
% \subsection{Miscellaneous}
@@ -2489,7 +2483,7 @@
\begingroup
\toks@\@xp{\MFPprogram@returns}%
\MFP@endgroup@after{%
- \MFP@def@after\MFPprogram@returns{\the\toks@ \MFP@def@after#1{#1}}%
+ \MFP@afterdef\MFPprogram@returns{\the\toks@ \MFP@afterdef#1{#1}}%
}}%
\def\MFP@ExportStack{\MFP@Export\MFP@Rstack}%
% \end{macrocode}
@@ -2514,7 +2508,7 @@
% strictly needed, allowing the parts to be specified as anything \TeX{}
% recognizes as a number and allowing any register name. This generality
% might reduce efficiency but it simplifies code. Because register
-% \texttt{z} is by far the most common one to load, we make more efficient
+% \reg{z} is by far the most common one to load, we make more efficient
% version of it.
% \begin{macrocode}
\def\MFP@Rload #1#2#3#4{%
@@ -2638,40 +2632,10 @@
%
% \bigskip
% The user could easily convert between radians and degrees using
-% multiplication and/or division. The commands \cs{Rdeg}, \cs{Rrad},
-% etc., aim to be a little more accurate.
-%
-% \subsection{Additional errors}
-%
-% These extra commands come with a few possible warnings and errors.
-%
-% \DescribeMacro{\LogOfZeroInt}
-% \DescribeMacro{\LogOfZeroFrac}
-% Trying to take the logarithm of $0$ will result in an error message.
-% If one allows \TeX{} to continue, the returned value will be negative,
-% with an integer part equal to the contents of \cs{LogOfZeroInt} and a
-% fractional part equal to the contents of \cs{LogOfZeroFrac}. The
-% defaults are both $99999999$.
-%
-% Trying to take the logarithm of a negative number will produce the
-% warning
-% \begin{verbatim}
-% MFP warning: Log of a negative number is complex.
-% Computing real part only. \end{verbatim}
-% The log of the absolute value is returned.
-%
-% Trying to take the square root of a negative number results in an
-% error. If you continue after the error message, the value $0$ is
-% returned.
-%
-% \SpecialUsageIndex{\MaxRealInt}
-% \SpecialUsageIndex{\MaxRealFrac}
-% Trying to take the exponential of a number larger than about $18.42$
-% will cause an error and the number returned has integer part
-% $99999999$ and fractional part $99999999$.
-%
-% Trying to take a negative power of $0$ returns the same as trying
-% to divide $1$ by $0$.
+% multiplication and/or division. One could similarly convert between
+% natural logarithms and base ten logarithms. The commands \cs{Rdeg},
+% \cs{Rrad}, \cs{Rlog} and \cs{Rln} (and their \cs{MFP...} counterparts)
+% aim for more accurate results.
%
% \subsection{Loading the extras}
%
@@ -2682,10 +2646,34 @@
% \DescribeMacro{\Rexp}\DescribeMacro{\Rsqrt}
% \DescribeMacro{\Rpow}
% We start \file{mfpextra} with the hook \cs{MFP@Rextra} that
-% \cs{starMFPprogram} will call to make available the extra operations
-% defined here.
+% \cs{startMFPprogram} will call to make available the extra operations
+% defined here. If \file{minifp.sty} has been loaded, this macro is
+% \cs{@empty}, otherwise it should be undefined. If it is undefined we
+% load \file{minifp.sty}. If it is then not \cs{@empty} we assume
+% \file{mfpextra.tex} was previously loaded and end input here.
% \begin{macrocode}
%<*extra>
+% check if mfpextra already loaded:
+\expandafter\ifx\csname MFP@xfinish\endcsname\relax
+\else \expandafter\endinput\fi
+\expandafter\edef\csname MFP@xfinish\endcsname{%
+ \catcode64=\the\catcode64 \space
+ \catcode46=\the\catcode46 \space
+ \catcode60=\the\catcode60 \space
+ \catcode62=\the\catcode62 \space}%
+\catcode64=11 % @
+\catcode46=12 % . (period)
+\catcode60=12 % <
+\catcode62=12 % >
+\ifx\MFP@Rextra\UndEfInEd \input minifp.sty \fi
+\ifx\MFP@Rextra\@empty
+\else
+ \immediate\write16{mfpextra.tex: already loaded.^^J}%
+ \MFP@xfinish
+ \expandafter\endinput
+\fi
+\immediate\write16{%
+ mfpextra.tex: extra operations for the MiniFP package.^^J}%
\def\MFP@Rextra{%
\def\Rcos {\MFP@stack@Unary\MFP@Rcos }%
\def\Rsin {\MFP@stack@Unary\MFP@Rsin }%
@@ -2720,7 +2708,36 @@
%
% \subsection{Error messages}
%
-% Messages for errors related to impossible powers, roots and logarithms.
+% These extra commands come with a few possible new warnings and errors.
+%
+% \DescribeMacro{\LogOfZeroInt}
+% \DescribeMacro{\LogOfZeroFrac}
+% Trying to take the logarithm of zero will result in an error message.
+% If one allows \TeX{} to continue, the returned value will be negative,
+% with an integer part whose absolute value is equal to the contents of
+% \cs{LogOfZeroInt} and a fractional part equal to the contents of
+% \cs{LogOfZeroFrac}. The defaults are both $99999999$.
+%
+% Trying to take the logarithm of a negative number will produce the
+% warning
+% \begin{verbatim}
+% MFP warning: Log of a negative number is complex.
+% Only the real part will be computed. \end{verbatim}
+% The log of the absolute value is returned.
+%
+% Trying to take the square root of a negative number has similar
+% behavior. It produces a warning and returns $0$.
+%
+% \SpecialUsageIndex{\MaxRealInt}
+% \SpecialUsageIndex{\MaxRealFrac}
+% Trying to take the exponential of a number larger than about $18.42$
+% will cause an error and the number returned has integer part
+% $99999999$ and fractional part $99999999$.
+%
+% Trying to take a negative power of $0$ produces an error and returns
+% the same value as trying to divide $1$ by $0$.
+%
+% Messages for errors related to impossible powers and logarithms.
% \begin{macrocode}
\def\MFP@logofzero@err{%
\MFP@errmsg{logarithm of zero}%
@@ -2731,13 +2748,9 @@
\def\LogOfZeroFrac{\MaxRealFrac}%
\def\MFP@expoverflow@err{%
\MFP@errmsg{Power too large}%
- {The power you tried to calcualate is too large for %
+ {The power you tried to calculate is too large for %
8 digits. If you continue, ^^Jthe value assigned will be %
\MaxRealInt.\MaxRealFrac.}}%
-\def\MFP@sqrtofneg@err{%
- \MFP@errmsg{square root of a negative number}%
- {You tried to take the square root of a negative value. What %
- were you thinking? If you continue, zero will be returned.}}%
\def\MFP@badpower@err{%
\MFP@errmsg{negative power of zero}%
{You tried to take a negative power of zero. What were you
@@ -2745,10 +2758,26 @@
\xOverZeroInt.\xOverZeroFrac.}}%
% \end{macrocode}
%
+% A debugging utility, \cs{MFPshowreg} displays the contents of a
+% register.
+% \begin{macrocode}
+\def\MFPshowreg #1{%
+\ifMFPdebug
+\begingroup
+ \edef\theregister{%
+ #1 = \expandafter \MFP@Sign
+ \csname MFP@#1@Sgn\endcsname %
+ \csname MFP@#1@Int\endcsname.%
+ \csname MFP@#1@Frc\endcsname}%
+ \show\theregister
+\endgroup
+\fi}%
+% \end{macrocode}
+%
% \subsection{Sine and Cosine}
%
-% For iterated code, the most common register to copy is \texttt{z} and
-% the most common place to copy it is to \texttt{x} or \texttt{y} so we
+% For iterated code, the most common register to copy is $z$ and
+% the most common place to copy it is to $x$ or $y$ so we
% make single commands to do those.
% \begin{macrocode}
\def\MFP@Rcopyz#1{\MFP@Rload {#1}\MFP@z@Sgn\MFP@z@Int\MFP@z@Frc}%
@@ -2764,34 +2793,72 @@
% \begin{verbatim}
% \MFPdeg\X\Y \MFPsin\Y\S \end{verbatim}
%
-% Our degree/radian conversions try to be more accurate than a simple
-% multiplication by $57.2957 7951$ or $0.0174 5329$. These conversion
-% factors are accurate to only eight digits, and the rounding error is
-% magnified by multiplication. Thus we will use 16 digits for these
-% constants. That is, we multiply first by $57.2957 7951$, then by the
-% next eight digits ($.30823208\times 10^{-8}$), performing the ``${}\times
-% 10^{-8}$'' by using the (rounded) integer part as the fractional part.
-%
-% The copying of \texttt{x} to \texttt{t} beforehand is so that we don't have to
-% remember which operations (in this case only \cs{MFP@Rmul} and
-% \cs{MFP@Radd}) leave register \texttt{x} unchanged. All operations
-% defined in \file{mfpextra.tex} overwrite the \texttt{x} register, as do
-% a few others.
-% \begin{macrocode}
-\def\MFP@Rdeg{%
- \MFP@Rcopy xt\MFP@Rload y1{57}{29577951}\MFP@Rmul
- \MFP@Rcopyz s\MFP@Rcopy tx\MFP@Rload y10{30823209}\MFP@Rmul
- \MFP@tempa\MFP@z@Int\relax
- \ifnum\MFP@z@Frc<50000000 \else \advance\MFP@tempa 1 \fi
- \MFP@Rload x{\ifnum\MFP@tempa>0 \MFP@z@Sgn\else0\fi}0\MFP@tempa
- \MFP@Rcopy sy\MFP@Radd}%
-\def\MFP@Rrad{%
- \MFP@Rcopy xt\MFP@Rload y10{0174 5329}\MFP@Rmul
- \MFP@Rcopyz s\MFP@Rcopy tx\MFP@Rload y10{25199433}\MFP@Rmul
- \MFP@tempa\MFP@z@Int\relax
- \ifnum\MFP@z@Frc<50000000 \else \advance\MFP@tempa 1 \fi
- \MFP@Rload x{\ifnum\MFP@tempa>0 \MFP@z@Sgn\else0\fi}0\MFP@tempa
- \MFP@Rcopy sy\MFP@Radd}%
+% For unit conversions such as radian to degree we try to be more accurate
+% than a multiplication by an eight-digit conversion factor allows.
+% If $x$ is large and the factor is off by $0.5\times 10^{-8}$, then the
+% result can be significantly off. But if we are able to give the
+% conversion factor 16 digits precision, then only the imprecision of $x$
+% will significantly affect the result.
+%
+% We express the conversion factor as an integer part and two eight-digit
+% fractional parts. We multiply $x$ by the integer and first fractional
+% part (\arg1 and \arg2) with a normal \cs{MFP@Rmul}, but we save the
+% underflow digits and undo the rounding that occured at the 8th digit.
+% Together these give us an essentially exact result. Then we multiply $x$
+% by the second fractional part (\arg3) and add the saved underflow to the
+% result. Finally, we round and add the result to the first product.
+% Argument \arg3, as well as the underflow digits, represent numbers less
+% than $10^{-8}$, so we effectively scale them up by $10^8$, round the
+% result to an integer and scale that back down.
+%
+% The registers $w$ and $v$ are used to save intermediate results.
+% The ``\texttt{DP}'' in \cs{MFP@DPmul} refers to the fact that we are
+% multiplying by a ``double precision'' real.
+% \begin{macrocode}
+\def\MFP@DPmul#1#2#3{%
+ \ifnum\MFP@x@Sgn=0
+ \MFP@Rzero
+ \else
+ \MFP@Rcopy xv%
+ \MFP@Rload y1{#1}{#2}\MFP@Rmul
+ \edef\MFP@w@Und{\MFP@z@Und}%
+ \ifnum\MFP@z@Frc@iii>4999
+ \MFP@tempa\MFP@z@Frc \advance\MFP@tempa-1
+ \edef\MFP@z@Frc{\number\MFP@tempa}%
+ \makeMFP@eightdigits\MFP@z@Frc
+ \fi
+ \MFP@Rcopyz w%
+ \MFP@Rcopy vx\MFP@Rload y10{#3}\MFP@Rmul
+ \MFP@Rcopyzx\MFP@Rload y\MFP@v@Sgn 0{\MFP@w@Und}\MFP@Radd
+ \MFP@tempa\MFP@z@Int\relax
+ \ifnum\MFP@z@Frc<50000000 \else \advance\MFP@tempa 1 \fi
+ \ifnum\MFP@tempa<\MFP@ttteight\relax
+ \MFP@Rload x{\ifnum\MFP@tempa>0 \MFP@z@Sgn\else0\fi}0\MFP@tempa
+ \else
+ \MFP@Rload x\MFP@z@Sgn10%
+ \fi
+ \MFP@Rcopy wy\MFP@Radd
+ \fi}%
+% \end{macrocode}
+%
+% Conversion factors:
+% \begin{itemize}
+% \item radians to degrees: $57.2957795130823209$
+% \item degrees to radians: $0.0174532925199433$
+% \item natural log to common log: $0.4342944819032518$
+% \item common log to natural log: $2.3025850929940457$
+% \end{itemize}
+%
+% Note that the comparatively large size of the first number means that
+% the $\pm0.5\cdot10^{-8}$ imprecision that $x$ implicitly carries will
+% be multiplied to approximately $\pm29.6\cdot 10^{-8}$ in the result.
+% The only way around this would be to operate with higher precision
+% internally. We do that in the code for computing angles.
+% \begin{macrocode}
+\def\MFP@Rdeg{\MFP@DPmul{57}{29577951}{30823209}}%
+\def\MFP@Rrad{\MFP@DPmul{0}{01745329}{25199433}}%
+\def\MFP@RbaseX{\MFP@DPmul{0}{43429448}{19032518}}%
+\def\MFP@RbaseE{\MFP@DPmul{2}{30258509}{29940457}}%
% \end{macrocode}
%
% There are very few angles that are expressible in eight digits whose sine
@@ -2813,71 +2880,59 @@
\MFP@Rcopyzx\MFP@Rsin}%
% \end{macrocode}
%
-% Return $0$ if $x$ is $0$. Then reduce $|x|$ by subtracting $360$ from
-% the integer part until it is less than $360$. Of course,
-% $\sin x = \sgn(x)\sin(|x|)$ so we only need to compute $\sin(|x|)$. For
-% $|x| < 180$, the sign will be that of $x$. For $|x| >= 180$, calculate
-% $\sin(|x| - 180)$ and the sign will be opposite that of $x$. Then again
-% return $0$ if $x$ is $0$. Finally, call \cs{MFP@@Rsin} to finish.
+% Reduce $|x|$ by subtracting $180$ from the integer part until it is less
+% than $180$. Of course, $\sin x = \sgn(x)\sin(|x|)$ so we only need to
+% compute $\sin(|x|)$. The sign will be that of $x$; each reduction by
+% $180$ changes the sign, but the reduction code keeps track of that. If
+% $x$ is 0 after the reduction, return zero.
% \begin{macrocode}
\def\MFP@Rsin{%
- \ifnum \MFP@x@Sgn=0 % return 0
- \MFP@Rzero
- \else
- \MFP@tempa\MFP@x@Int\relax
- \MFP@reduce@angle
- \ifnum \MFP@tempa<180
- \edef\MFP@sin@Sgn{\MFP@x@Sgn}%
- \else
- \edef\MFP@sin@Sgn{\number-\MFP@x@Sgn}%
- \advance \MFP@tempa-180
- \fi
- \MFP@tempb\MFP@x@Frc\relax
- \ifnum\MFP@tempa>0
- \MFP@@Rsin
- \else\ifnum\MFP@tempb>0
- \MFP@@Rsin
- \else
- \MFP@Rzero
- \fi\fi
- \fi}%
+ \MFP@tempa\MFP@x@Int
+ \MFP@tempb\MFP@x@Frc
+ \MFP@tempc\MFP@x@Sgn\relax
+ \MFP@reduce@angle
+ \ifnum\MFP@tempa>0 \MFP@@Rsin
+ \else\ifnum\MFP@tempb>0 \MFP@@Rsin
+ \else \MFP@Rzero
+ \fi\fi}%
% \end{macrocode}
%
-% This following reduces $|x|$ to the case $0 <= |x| < 360$. It assumes
-% the integer part is in count register \cs{MFP@tempa}.
+% This following reduces $|x|$ to the case $0 \le |x| < 180$. It assumes
+% the integer part is in count register \cs{MFP@tempa}, the sign in
+% \cs{MFP@tempc}.
% \begin{macrocode}
\def\MFP@reduce@angle{%
- \ifnum\MFP@tempa<360
+ \ifnum\MFP@tempa<180
\else
- \advance\MFP@tempa-360
+ \advance\MFP@tempa-180
+ \MFP@tempc-\MFP@tempc
\@xp\MFP@reduce@angle
\fi}%
% \end{macrocode}
%
% At this point, $|x|$ is represented by \cs{MFP@tempa} (integer part) and
-% \cs{MFP@tempb} (fractional part). Also, we already know the sign, stored
-% in \cs{MFP@sin@Sgn}. Moreover $0 < {}$\cs{MFP@tempa}${} < 180$. We now
+% \cs{MFP@tempb} (fractional part). Also, we already know the sign stored
+% in \cs{MFP@tempc}. Moreover $0 < {}$\cs{MFP@tempa}${} < 180$. We now
% reduce to $0 < |x| \le 90$ using $\sin(x) = \sin(180-|x|)$, and return
% $1$ if equal to $90$.
%
% The calculation of $180-x$ is optimized, taking advantage of the fact
% that both $x$ and the result are known to be positive. If the fractional
-% part is positive, we would normally borrow $1$ by reducing $180$, but
-% instead we increase the integer part of $x$ by one.
+% part is positive, we borrow $1$ by reducing $180$ to $179$.
% \begin{macrocode}
\def\MFP@@Rsin{%
\ifnum\MFP@tempa<90
\else
+ \MFP@tempa -\MFP@tempa
\ifnum\MFP@tempb>0
- \advance\MFP@tempa 1
\MFP@tempb -\MFP@tempb
\advance\MFP@tempb \MFP@ttteight\relax
+ \advance\MFP@tempa 179
+ \else \advance\MFP@tempa 180
\fi
- \MFP@tempa -\MFP@tempa
- \advance\MFP@tempa 180
\fi
\ifnum\MFP@tempa=90
- \MFP@Rloadz \MFP@sin@Sgn10%
+ \MFP@Rloadz \MFP@tempc10%
\else
% \end{macrocode}
%
@@ -2900,11 +2955,12 @@
% computations amount to concatenating the top six digits of
% \cs{MFP@tempb} to the digits of \cs{MFP@tempa}. This will produce the
% integer form of the fractional part of $x/100$ (the integer part of
-% $x/100$ is $0$).
+% $x/100$ is zero).
%
% Division by $100$ can turn a number into $0$. This is one place we can
-% lose accuracy in the last digit of the result. In compensation, the rest
-% of the calculations become extremely accurate.
+% lose accuracy (up to $\pm1$ in the last digit of the result). In
+% compensation, the rest of the calculations become very much more
+% accurate.
% \begin{macrocode}
\advance\MFP@tempb 50 \divide\MFP@tempb 100
\multiply\MFP@tempa 1000000 \advance\MFP@tempb\MFP@tempa
@@ -2915,43 +2971,34 @@
%
% We save some multiplications by working with $t=x^2$. As we don't need
% the original $x$ anymore, we simply replace it with the newly reduced
-% value. We also save this reduced $x$ in another register, \texttt{s}, as
+% value. We also save this reduced $x$ in another register, $s$, as
% we will need it again at the end, and our intermediate calculations do
-% not preserve the \texttt{x} register. Then we square $x$ and, if that
+% not preserve the $x$ register. Then we square $x$ and, if that
% square is $0$ we can skip all the power series and simply return $x$
-% converted to radians (that's the last multiplication). If $x^2$ is not
-% $0$, we save it in temporary register \texttt{t} and call our power
-% series. When this program is finished, all that remains is a final
-% multiplication by a conversion factor\dots
+% converted to radians. If $x^2$ is not zero, we save it in temporary
+% register $t$ and call our power series. When this program is
+% finished, all that remains is the final multiplication by a conversion
+% factor (\cs{MFP@DPmul}).
% \begin{macrocode}
- \MFP@Rload x10\MFP@tempb
- \MFP@Rcopy xs%
+ \MFP@Rload s\MFP@tempc0\MFP@tempb
+ \MFP@Rcopy sx%
\MFP@Rsq
\ifnum \MFP@z@Frc>0
\MFP@Rcopyz t\MFP@Rsin@prog
\else
\MFP@Rcopy sx%
\fi
- \MFP@Rload y11{74532925}\MFP@Rmul
-% \end{macrocode}
-% \dots except this fiddle about the sign. Theoretically, all cases
-% where $\sin x$ can be $0$ were previously weeded out. However, I am not
-% 100 percent certain that rounding in \cs{MFP@Rsin@prog} will never
-% lead to a value of $0$.
-% \begin{macrocode}
- \ifnum\MFP@z@Sgn=0 \else
- \let\MFP@z@Sgn\MFP@sin@Sgn
- \fi
+ \MFP@DPmul 1{74532925}{19943296}%
\fi
\fi}%
% \end{macrocode}
%
% \cs{MFP@Rsin@prog} is the power series computation. The power series
-% need only go to eight terms as the ninth would be less than $.5*10^{-8}$ and
-% so our 8-place computations would return $0$. Our 8-term series is:
+% need only go to the $x^{13}$ term as the next is less than $10^{-9}$ and
+% in our 8-place computations is indistingushable from $0$. Our series is:
% $$
% rx(1 - r^2t/3! + r^4t^2/5! - r^6t^3/7! + r^8t^4/9! - r^{10}t^5/11! +
-% r^{12}t^6/13! - r^{14}t^7/15!)
+% r^{12}t^6/13!)
% $$
% where $r$ is the factor that converts $x$ to radian measure
% (hectodegrees to radians). When $x$ is so small as to produce $t = 0$ we
@@ -2959,32 +3006,37 @@
%
% We minimize any multiplications of tiny numbers by computing this as
% $$
-% r(1 - gt(1 - ft(1 - et(1 - dt(1 - ct(1 - bt(1 - at))))))).
+% rx(1 - ft(1 - et(1 - dt(1 - ct(1 - bt(1 - at)))))).
% $$
-% Now $r = 1.74532925$ and $a$, $b$, etc., have formulas:
+% In this format, additional terms might actually make a difference,
+% because $at$ is not particularly small. However, the more computations
+% we have, the more errors accumulate. Therefore we take the fewest that
+% produce acceptable accuracy.
+%
+% Now $r = 1.7453292519943296$ and $a$, $b$, etc., have formulas:
% $$
% \vcenter{\centering
-% $\displaystyle a = r^2/15/14,\ b = r^2/13/12,\ c = r^2/11/10,\ d = r^2/9/8$,\\
-% $\displaystyle e = r^2/7/6,\ f = r^2/5/4,\ g = r^2/3/2$.\par
+% $\displaystyle a = r^2/13/12,\ b = r^2/11/10,\ c = r^2/9/8$,\\
+% $\displaystyle d = r^2/7/6,\ e = r^2/5/4,\ f = r^2/3/2$.\par
% }
% $$
% An alternative method would be to accumulate a sum, computing each term
% from the previous one (e.g., if $u = t^3/7!$ is the fourth term, the next
% one is $u*t*(1/(8*9))$). This is a bit more complicated to code and requires
% moving values around more. It would have the advantage that we can stop
-% whenever a term evaluates to $0$, making computation faster for small
-% values of $x$.
+% whenever a term evaluates to zero, making computation faster for small
+% values of $x$. I have not determined whether it would compromise
+% accuracy.
%
% We avoid divisions by precomputing the coefficients $a$, $b$, $c$, etc.
% Note that without the reduction in $x$, the value of $a$ for example
-% would be $0.00000145$, with only three significant figures of accuracy.
+% would be $0.00000195$, with only three significant figures of accuracy.
% Now we can have seven, and the accuracy is more-or-less determined by that
% of the reduced x.
% $$
% \vcenter{\centering
-% $\displaystyle a = 0.01450559,\ b = 0.01952675,\ c = 0.02769249$,\\
-% $\displaystyle d = 0.04230797,\ e = 0.07252796,\ f = 0.15230871$,\\
-% $\displaystyle g = 0.50769570$.\par
+% $\displaystyle a = 0.01952676,\ b = 0.02769249,\ c = 0.04230797,$,\\
+% $\displaystyle d = 0.07252796,\ e = 0.15230871,\ f = 0.50769570$.\par
% }
% $$
% It is important to note that the following operations step all over
@@ -2995,28 +3047,20 @@
% previous operation. Instead of simply subtracting, we optimize based
% on the fact that $z$ is known to be nonnegative and not larger than $1$.
%
-% The macro \cs{MFP@com@iter} `flips' the previous result then multiplies
+% The macro \cs{MFP@com@iter} `flipz' the previous result then multiplies
% by $t$ and the indicated coefficient. (The name of this macro stands for
% ``common iterated'' code; it is reused for some other power series.)
%
% For extra efficiency, the power series uses a ``small'' version of
% multiplication \cs{MFP@Rsmul}, used only when the factors are sure to
-% lie in $[0,1]$.
-%
-% Despite what I said above, our chosen method of computation has a
-% slightly improved accuracy (in numerical experiments) if we take it one
-% step further, but adding only half the last term. This splits the
-% difference between stopping at the 15th power or continuing to the 17th.
-% It has provably better worst-case accuracy, but on average, who knows?
-% We are right at the edge of our 8-digit accuracy anyway. The constant
-% \texttt{00559959} corresponds to half of $r^2/16/17$.
+% lie in $[0,1]$. This does not take into account the sign of $x$,
+% whence the ending \cs{edef}.
% \begin{macrocode}
\def\MFP@Rsin@prog{%
- \MFP@Rcopy tx%
- \MFP@Rload y10{00559959}\MFP@Rsmul\MFP@com@iter{01450559}%
- \MFP@com@iter{01952675}\MFP@com@iter{02769249}\MFP@com@iter{04230797}%
- \MFP@com@iter{07252796}\MFP@com@iter{15230871}\MFP@com@iter{50769570}%
- \MFP@flipz \MFP@Rcopyzx \MFP@Rcopy sy\MFP@Rsmul \MFP@Rcopyzx}%
+ \MFP@Rcopy tx\MFP@Rload y10{01952676}\MFP@Rsmul%
+ \MFP@com@iter{02769249}\MFP@com@iter{04230797}\MFP@com@iter{07252796}%
+ \MFP@com@iter{15230871}\MFP@com@iter{50769570}\MFP@flipz \MFP@Rcopyzx
+ \MFP@Rcopy sy\MFP@Rsmul\MFP@Rcopyzx\edef\MFP@x@Sgn{\MFP@s@Sgn}}%
\def\MFP@flipz{%
\ifnum\MFP@z@Sgn=0
\MFP@Rloadz 110%
@@ -3053,18 +3097,16 @@
% multiplication could round to $0$ and then, after subtraction, a $1$
% could occur. We handle those easy cases separately, so that in
% \cs{MFP@@Rsmul} we don't have to worry about the integer parts at all.
+%
+% Also, since these are completely internal, we don't even define the
+% overflow and underflow macros.
% \begin{macrocode}
\def\MFP@Rsmul{%
- \ifnum \MFP@x@Sgn=0
- \MFP@Rzero
- \else\ifnum \MFP@y@Sgn=0
- \MFP@Rzero
- \else\ifnum\MFP@x@Int>0
- \MFP@Rcopy yz%
- \else\ifnum\MFP@y@Int>0
- \MFP@Rcopy xz%
- \else
- \MFP@@Rsmul
+ \ifnum \MFP@x@Sgn=0 \MFP@Rzero
+ \else\ifnum \MFP@y@Sgn=0 \MFP@Rzero
+ \else\ifnum\MFP@x@Int>0 \MFP@Rcopy yz%
+ \else\ifnum\MFP@y@Int>0 \MFP@Rcopy xz%
+ \else \MFP@@Rsmul
\fi\fi\fi\fi}%
\def\MFP@@Rsmul{%
\MFP@split\MFP@x@Frc\MFP@x@Frc@i\MFP@x@Frc@ii
@@ -3084,15 +3126,11 @@
\advance\MFP@tempb1
\edef\MFP@z@Frc@ii{\number\MFP@tempb}\fi
\MFP@carrym\MFP@z@Frc@ii\MFP@z@Frc@i
- \makeMFP@fourdigits\MFP@z@Frc@iv
- \makeMFP@fourdigits\MFP@z@Frc@iii
\makeMFP@fourdigits\MFP@z@Frc@ii
\makeMFP@fourdigits\MFP@z@Frc@i
\def\MFP@z@Int{0}%
\edef\MFP@z@Frc{\MFP@z@Frc@i\MFP@z@Frc@ii}%
- \edef\MFP@z@Sgn{\ifnum\MFP@z@Frc=0 0\else 1\fi}%
- \edef\MFP@z@Und{\MFP@z@Frc@iii\MFP@z@Frc@iv}%
- \edef\MFP@z@Ovr{0}}%
+ \edef\MFP@z@Sgn{\ifnum\MFP@z@Frc=0 0\else 1\fi}}%
% \end{macrocode}
%
% \subsection{Polar angle}
@@ -3101,7 +3139,8 @@
% more general angle function. This is a binary operation that accepts
% the two coordinates of a point and computes its angle in polar
% coordinates. One then has, for example, $\arctan x =
-% \mathop{\rm angle}(1,x)$ and $\arccos x = \angle(x, \sqrt{1-x^2})$.
+% \mathop{\mathrm{angle}}(1,x)$ and $\arccos x = \mathop{\mathrm{angle}}
+% (x, \sqrt{1-x^2})$.
%
% We start, as usual, with a few reductions. When the $y$-part is $0$, we
% immediately return $0$ or $180$. If the $y$-part is negative, we compute
@@ -3109,11 +3148,12 @@
% compute the angle for $|x|$ and subtract it from $180$. Finally,
% reduced to both coordinates positive, if $y>x$ we compute the angle of
% $(y,x)$ and subtract that from $90$. Ultimately, we apply a power
-% series formula for $\angle(1,y/x)$ and get convergence when the
-% argument is less than $1$, but convergence is poor unless the argument
-% is less than $1/2$. When that is not the case, conceptually, we rotate
-% the picture clockwise by the arctangent of $1/2$, compute the angle of
-% the new point and then add a precomputed value of $\arctan(1/2)$.
+% series formula for $mathop{\mathrm{angle}}(1,y/x)$ and get convergence
+% when the argument is less than $1$, but convergence is poor unless the
+% argument is less than $1/2$. When that is not the case, conceptually, we
+% rotate the picture clockwise by the arctangent of $1/2$, compute the
+% angle of the new point and then add a precomputed value of
+% $\arctan(1/2)$.
% \begin{macrocode}
\def\MFP@Rangle{%
\ifcase\MFP@y@Sgn\relax
@@ -3139,18 +3179,16 @@
\def\MFP@x@Sgn{1}\MFP@@Rangle
\MFP@Rcopyzy\MFP@Rload x1{180}0\MFP@Rsub
\fi
- \let\MFP@z@Sgn\MFP@angle@Sgn
- }}%
+ \let\MFP@z@Sgn\MFP@angle@Sgn}}%
\def\MFP@@Rangle{%
\MFP@Rcmp
\ifMFP@neg
- \MFP@Rcopy xs\MFP@Rcopy yx\MFP@Rcopy sy%
+ \MFP@Rcopy xw\MFP@Rcopy yx\MFP@Rcopy wy%
\MFP@@@Rangle
\MFP@Rload x1{90}0\MFP@Rcopyzy\MFP@Rsub
\else
\MFP@@@Rangle
- \fi
-}%
+ \fi}%
% \end{macrocode}
%
% Precisely what we do when we are finally in the case $0<y<x$ is perform
@@ -3167,62 +3205,56 @@
% If we still have $y/x > 1/4$, we perform $(x'',y'') = (4x + y, 4y - x)$,
% which then satisfies $y''/x'' \le 1/4$. When either of these
% transformations is performed, we add the corresponding angle to the
-% ``angle-so-far'' in register \texttt{a}.
+% ``angle-so-far'' in register $a$.
%
% We could continue this iteration 32 times to get (theoretically) the
% angle in degrees to $\pm 10^{-8}$. That seems a bit long, plus the
% accumulation of errors over $32$ iterations could (in the worst case)
% produce less than $\pm10^{-7}$ accuracy.
%
-% To get the accuracy we need we work in ``scaled reals''. That is, we
-% get 10 decimal places of accuracy by letting two digits of the integer
-% part represent the first two digits after the decimal point, and the
-% eight digits of the fractional part represent digits 3 through 10 after the
-% point. The macro \cs{MFP@RmulC} (around line 19 of the definition of
-% \cs{MFP@@@Rangle}) is a quick multiplication by 100, converting the
-% argument of the arctangent command to a scaled real.
-%
-% Since we ultimately take the arctangent of $y/x$ we get best accuracy if
-% $y$ is first scaled, which means it must necessarily be smaller than
-% $10^6$ to start with. Our initial reductions can double its size, so we
-% first make sure it is smaller than $500\,000$, dividing by $200$. This
-% can actually reduce accuracy, but as it is only done when $y$ already
-% has $14$ significant figures, and since the quotient $100y/x$ has at
-% most $10$, the loss is not significant.
+% To get the accuracy we need, we work in ``scaled reals''. That is, we
+% get 10 effective decimal places of accuracy by letting an $x$ in the
+% range $0< x < 100$ stand for $0< x/100 < 1$.
+%
+% Our initial reductions can increase $x$ by a factor of about 13.
+% Moreover, we ultimately need to scale y by 100 when we convert to
+% scaled computations. Thus, if we make sure $x$ is less than
+% $1\,000\,000$, we will prevent overflow in both cases.
% \begin{macrocode}
\def\MFP@Rquad{\MFP@Rdbl\MFP@Rcopyzx\MFP@Rdbl}%
\def\MFP@@@Rangle{%
\MFP@Rcopy xs\MFP@Rcopy yt%
- \ifnum\MFP@y@Int<500000
+ \ifnum\MFP@x@Int<1000000
\else
- \MFP@Rload y1{200}0\MFP@Rdiv
- \MFP@Rcopyz s\MFP@Rcopy tx%
- \MFP@Rload y1{200}0\MFP@Rdiv
- \MFP@Rcopyz t%
+ \MFP@RdivC \MFP@Rcopyz s%
+ \MFP@Rcopy tx\MFP@RdivC \MFP@Rcopyz t%
\fi
- \MFP@Rcopy tx\MFP@Rdbl\MFP@Rcopyzx\MFP@Rcopy sy\MFP@Rcmp
- \ifMFP@pos
- \MFP@Rsub\MFP@Rcopyz u\MFP@Rcopy sx\MFP@Rdbl
- \MFP@Rcopyzx\MFP@Rcopy ty\MFP@Radd
- \MFP@Rcopyz s\MFP@Rcopy ut%
- \MFP@Rload a1{2656}{50511771}%
+ \ifnum\MFP@t@Sgn=0 \MFP@Rzero
\else
- \MFP@Rload a000%
- \fi
- \MFP@Rcopy tx\MFP@Rquad\MFP@Rcopyzx\MFP@Rcopy sy\MFP@Rcmp
- \ifMFP@pos
- \MFP@Rsub\MFP@Rcopyz u\MFP@Rcopy sx\MFP@Rquad
- \MFP@Rcopyzx\MFP@Rcopy ty\MFP@Radd
- \MFP@Rcopyz s\MFP@Rcopy ut%
- \MFP@Rcopy ax\MFP@Rload y1{1403}{62434679}%
- \MFP@Radd\MFP@Rcopy za%
- \fi
- \MFP@Rcopy tx\MFP@RmulC
- \MFP@Rcopyzx\MFP@Rcopy sy\MFP@Rdiv
- \MFP@Rcopyzx\MFP@Ratanc
- \MFP@Rcopyzx\MFP@Rdeg
- \MFP@Rcopyzx\MFP@Rcopy ay\MFP@Radd
- \MFP@Rcopyzx\MFP@RdivC}%
+ \MFP@Rcopy tx\MFP@Rdbl\MFP@Rcopyzx\MFP@Rcopy sy\MFP@Rcmp
+ \ifMFP@pos
+ \MFP@Rsub\MFP@Rcopyz u\MFP@Rcopy sx\MFP@Rdbl
+ \MFP@Rcopyzx\MFP@Rcopy ty\MFP@Radd
+ \MFP@Rcopyz s\MFP@Rcopy ut%
+ \MFP@Rload a1{2656}{50511771}%
+ \else
+ \MFP@Rload a000%
+ \fi
+ \MFP@Rcopy tx\MFP@Rquad\MFP@Rcopyzx\MFP@Rcopy sy\MFP@Rcmp
+ \ifMFP@pos
+ \MFP@Rsub\MFP@Rcopyz u\MFP@Rcopy sx\MFP@Rquad
+ \MFP@Rcopyzx\MFP@Rcopy ty\MFP@Radd
+ \MFP@Rcopyz s\MFP@Rcopy ut%
+ \MFP@Rcopy ax\MFP@Rload y1{1403}{62434679}%
+ \MFP@Radd\MFP@Rcopy za%
+ \fi
+ \MFP@Rcopy tx\MFP@RmulC
+ \MFP@Rcopyzx\MFP@Rcopy sy\MFP@Rdiv
+ \MFP@Rcopyzx\MFP@Ratanc
+ \MFP@Rcopyzx\MFP@Rdeg
+ \MFP@Rcopyzx\MFP@Rcopy ay\MFP@Radd
+ \MFP@Rcopyzx\MFP@RdivC
+ \fi}%
% \end{macrocode}
%
% Here are fast multiplication and division by 100. We need these because
@@ -3232,25 +3264,28 @@
% division by 100: $.5\times.5 = .25$ would be computed as $(50\times50) /
% 100 = 25$.
% \begin{macrocode}
+\def\MFP@twoofmany#1#2#3\MFP@end{#1#2}%
+\def\MFP@gobbletwo#1#2{}%
\def\MFP@RmulC{%
- \makeMFP@eightdigits\MFP@x@Frc
- \edef\MFP@Tmp{\number\MFP@x@Int.\MFP@x@Frc}%
- \@xp\MFP@@RmulC\MFP@Tmp\mfp@end}%
-\def\MFP@@RmulC#1.#2#3#4\mfp@end{%
- \MFP@Rloadz\MFP@x@Sgn{#1#2#3}{#400}}%
+ \edef\MFP@z@Int{\MFP@x@Int\@xp\MFP@twoofmany\MFP@x@Frc\MFP@end}%
+ \edef\MFP@z@Frc{\@xp\MFP@gobbletwo\MFP@x@Frc00}%
+ \edef\MFP@z@Sgn{\MFP@x@Sgn}}%
\def\MFP@RdivC{%
\makeMFP@eightdigits\MFP@x@Int
\makeMFP@eightdigits\MFP@x@Frc
- \@XP\MFP@@RdivC\@xp\MFP@x@Int\MFP@x@Frc\mfp@end}%
+ \@XP\MFP@@RdivC\@xp\MFP@x@Int\MFP@x@Frc\MFP@end}%
\def\MFP@@RdivC#1#2#3#4#5#6{%
\edef\MFP@z@Int{\number#1#2#3#4#5#6}%
\MFP@@@RdivC}%
-\def\MFP@@@RdivC#1#2#3#4#5#6#7#8#9\mfp@end{%
+\def\MFP@@@RdivC#1#2#3#4#5#6#7#8#9\MFP@end{%
\MFP@tempa#1#2#3#4#5#6#7#8\relax
\ifnum#9>49 \advance\MFP@tempa1 \fi
+ \edef\MFP@z@Frc{\number\MFP@tempa}%
+ \makeMFP@eightdigits\MFP@z@Frc
\edef\MFP@z@Sgn{\MFP@x@Sgn}%
- \ifnum\MFP@tempa=0 \ifnum\MFP@x@Int=0 \def\MFP@z@Sgn{0}\fi\fi
- \MFP@Rloadz\MFP@z@Sgn\MFP@z@Int\MFP@tempa}%
+ \ifnum\MFP@tempa=0
+ \ifnum\MFP@z@Int=0 \def\MFP@z@Sgn{0}\fi
+ \fi}%
% \end{macrocode}
%
% Finally, we compute the arctan of a scaled real producing a result
@@ -3264,12 +3299,12 @@
% the sum
% $$
% x\left(1 - \frac{u}{3} + \frac{u^2}{5} - \frac{u^3}{7} + \cdots
-% - \frac{u^7}{2\cdot 15}\right),
+% - \frac{u^7}{15}\right),
% $$
% where $u = x^2$.
%
% We start with the common iterated code. It assumes a scaled value in x
-% to be multiplied by the saved (scaled) value of $x^2$ (in register u)
+% to be multiplied by the saved (scaled) value of $x^2$ (in register $u$)
% and by a coefficient (supplied in separate integer and fractional
% parts). It ends with the new value in x.
% \begin{macrocode}
@@ -3280,22 +3315,25 @@
\MFP@Rcopyzy\MFP@Rload x1{100}{00000000}%
\MFP@Rsub\MFP@Rcopyzx}%
\def\MFP@Ratanc{%
- \MFP@Rcopy xs%
- \MFP@Rcopy xy\MFP@scaledmul\MFP@Rcopyz u%
- \MFP@Rcopyzx\MFP@Rload y1{86}{66666667}%
- \MFP@scaledmul
- \MFP@Rcopyzy\MFP@Rload x1{100}{00000000}%
- \MFP@Rsub\MFP@Rcopyzx
- \MFP@atan@iter{84}{61538462}\MFP@atan@iter{81}{81818182}%
- \MFP@atan@iter{77}{77777778}\MFP@atan@iter{71}{42857143}%
- \MFP@atan@iter{60}{00000000}\MFP@atan@iter{33}{33333333}%
- \MFP@Rcopy sy\MFP@scaledmul}%
+ \MFP@Rcopy xs\MFP@Rcopy xy\MFP@scaledmul
+ \ifnum \MFP@z@Sgn=0
+ \MFP@Rcopy sz%
+ \else
+ \MFP@Rcopyz u\MFP@Rcopyzx
+ \MFP@Rload y1{86}{66666667}\MFP@scaledmul
+ \MFP@Rcopyzy\MFP@Rload x1{100}{00000000}\MFP@Rsub\MFP@Rcopyzx
+ \MFP@atan@iter{84}{61538462}\MFP@atan@iter{81}{81818182}%
+ \MFP@atan@iter{77}{77777778}\MFP@atan@iter{71}{42857143}%
+ \MFP@atan@iter{60}{00000000}\MFP@atan@iter{33}{33333333}%
+ \MFP@Rcopy sy\MFP@scaledmul
+ \fi}%
% \end{macrocode}
%
% \subsection{Logarithms}
%
-% Now for logarithms. We are going to compute a base 10 logarithm. This
-% allows the first step of the calculation to be essentially trivial: to
+% Now for logarithms. We are going to compute both common logarithms
+% (base $10$) and natural logarithms (base $e$). The first step of the
+% calculation is be essentially trivial and works with base 10: to
% get the integer part of the log for numbers with positive integer part,
% count the digits in the integer part and subtract $1$. For numbers less
% than one, count the number of zeros at the beginning of the fractional
@@ -3306,18 +3344,27 @@
% an acceptable number of of terms. Then we proceed as in the code for
% sine.
%
-% We keep the value-so-far in register \texttt{s} and the modified
-% $x$-value in register \texttt{t}. The log of \texttt{t} need only be
-% added to register \texttt{s} for the final result.
+% The power series produces a logarithm in base $e$ so we ultimately get
+% the answer in two parts, with the parts calculated for different bases.
+% The last step is to multiply the second part by a conversion factor and
+% add the first to it. For natural log, convert the first and add the
+% second. Which one is to be returned is passed as a boolean
+%
+% We keep the value-so-far in register $s$ and the modified
+% $x$-value in register $t$.
% \begin{macrocode}
-\def\MFP@Rlog{%
+\newif\ifMFP@natural
+\def\MFP@Rlog{\MFP@naturalfalse\MFP@Rlog@}%
+\def\MFP@Rln{\MFP@naturaltrue\MFP@Rlog@}%
+\def\MFP@Rlog@{%
\ifnum\MFP@x@Sgn=0
\MFP@logofzero@err
\MFP@Rloadz{-1}\LogOfZeroInt\LogOfZeroFrac
\else
\ifnum \MFP@x@Sgn<0
- \MFP@warn{The logarithm of a negative number is complex.%
+ \MFP@warn{The logarithm of a negative number is complex.
\MFP@msgbreak Only the real part will be computed}%
+ \def\MFP@x@Sgn{1}%
\fi
\MFP@Rload s000%
% \end{macrocode}
@@ -3326,115 +3373,113 @@
% number of places that will be shifted in \cs{MFP@tempa}. We use
% \cs{number} to strip the leading zeros and (essentially) we count
% the number of digits that remain. Then we shift left, putting the first
-% digit into the integer part of \texttt{s} and the rest into the
-% fractional part. Despite its name, \cs{MFP@numzeros} actually returns
-% one more than the number of zeros, which is the number of places to
-% shift, and the integer part of the logarithm.
+% digit into the integer part of \reg{s} and the rest into the
+% fractional part.
% \begin{macrocode}
\ifnum \MFP@x@Int=0
\edef\MFP@x@Tmp{\number\MFP@x@Frc}%
- \MFP@tempa=\MFP@numzeros\MFP@x@Tmp\relax
+ \MFP@tempa=\MFP@numshiftL\MFP@x@Tmp\relax
\def\MFP@s@Sgn{-1}%
- \edef\MFP@t@Int{\@xp\MFP@firstofmany\MFP@x@Tmp\mfp@end}%
+ \edef\MFP@t@Int{\@xp\MFP@oneofmany\MFP@x@Tmp\MFP@end}%
\edef\MFP@t@Frc{\@xp\@gobble\MFP@x@Tmp0}%
- \MFPpadto@eight\MFP@t@Frc
+ \MFP@padtoeight\MFP@t@Frc
\else
% \end{macrocode}
% When the integer part is not $0$, we get the number of digits to
-% shift again in \cs{MFP@tempa}. We actually want one less than the
-% number of digits, so that is what \cs{MFP@numdigits} actually produces.
+% shift again in \cs{MFP@tempa}. It will be one less than the number of
+% integer digits.
% \begin{macrocode}
- \MFP@tempa\MFP@numdigits\MFP@x@Int
- \edef\MFP@x@Tmp{\MFP@x@Int\MFP@x@Frc}% gather all digits
- \edef\MFP@s@Sgn{1}%
- \edef\MFP@t@Int{\@xp\MFP@firstofmany\MFP@x@Tmp\mfp@end}%
+ \MFP@tempa\MFP@numshiftR\MFP@x@Int
+ \edef\MFP@x@Tmp{\MFP@x@Int\MFP@x@Frc}%
+ \ifnum\MFP@tempa>0 \def\MFP@s@Sgn{1}\fi
+ \edef\MFP@t@Int{\@xp\MFP@oneofmany\MFP@x@Tmp\MFP@end}%
\edef\MFP@x@Tmp{\@xp\@gobble\MFP@x@Tmp}%
- \edef\MFP@t@Frc{\@xp\MFP@firsteightofmany\MFP@x@Tmp\mfp@end}%
+ \edef\MFP@t@Frc{\@xp\MFP@eightofmany\MFP@x@Tmp\MFP@end}%
\fi
% \end{macrocode}
%
-% The integer part of $\log x$ is now known, so save it in value-so-far.
-% Also, set the sign of the reduced argument to positive. Then call
-% \cs{MFP@log@reduce}, which reduces $x$ to at most $10^{1/16} \approx
-% 1.155\,$. Finally, if the reduced $x$ is $1$, return the value so far,
-% otherwise call the power series program.
+% Now the integer part of $\log_{10} x$ is known. We save it in $s$
+% Also, set the sign of the reduced argument (positive). Then call
+% \cs{MFP@Rlog@reduce}, which reduces $x$ to less than $1.161\,$ while
+% possibly increasing $s$. For the natural log, we convert the value in
+% $s$.
+%
+% If the reduced $x$ is $1$, return the value in $s$, otherwise call the
+% power series program (discarding the integer part of $t$, which should
+% be a $1$). Finally, convert the returned result if necessary and add
+% register $s$ to it.
% \begin{macrocode}
\edef\MFP@s@Int{\number\MFP@tempa}%
\def\MFP@t@Sgn{1}%
\MFP@Rlog@reduce
+ \ifMFP@natural \MFP@Rcopy sx\MFP@RbaseE \MFP@Rcopy zs\fi
\ifnum\MFP@t@Frc=0
\MFP@Rcopy sz%
\else
- \MFP@Rlog@prog
+ \def\MFP@t@Int{0}\MFP@Rlog@prog
+ \ifMFP@natural\else \MFP@Rcopyzx \MFP@RbaseX \fi
+ \MFP@Rcopy sy\MFP@Rcopyzx\MFP@Radd
\fi
\fi}%
-\def\showreg #1{% for debugging
-\begingroup
- \edef\reg{%
- #1 = (\csname MFP@#1@Sgn\endcsname)%
- \csname MFP@#1@Int\endcsname.%
- \csname MFP@#1@Frc\endcsname}%
- \show\reg
-\endgroup}%
% \end{macrocode}
%
-% We count the number of digits by lining them up, followed by the
-% possible numbers, and picking out the ninth argument. We count leading
-% zeros in the fractional part by removing leading zeros, lining up the
-% remaining digits, followed by the possible numbers, and again picking
-% the ninth. In \cs{MFP@Rlog@reduce} we divide by the square root of 10 if
-% the number is larger than that (adding $.5$ to value-so-far). We repeat
-% with the 4th, 8th and 16th roots. For speed, instead of dividing, we
-% actually multiply by the precalculated reciprocal, passed as the third
-% argument.
+% We determine the size of a right shift by lining up the digits in
+% the integer part, followed by the possible numbers, and picking out the
+% ninth argument. Similarly, to get a left shift we line up the digits
+% of the fractional part (minus the leading zeros) followed by the
+% possible numbers, and again picking the ninth.
+% \begin{macrocode}
+\def\MFP@numshiftR#1{\@xp\MFP@ninthofmany#176543210\MFP@end}%
+\def\MFP@numshiftL#1{\@xp\MFP@ninthofmany#112345678\MFP@end}%
+% \end{macrocode}
+%
+% In \cs{MFP@Rlog@reduce} we divide by the square root of 10 if the number
+% is significantly larger than that (adding $.5$ to value-so-far). We
+% repeat with the 4th, 8th and 16th roots. It seems that this could be
+% where errors can accumulate, so the divisions are done with double
+% precision multiplication and $x$ is scaled by 100. Our check whether
+% $x > \sqrt{10}$ is rather rough: comparing the first three digits only,
+% but even in the worst case, the final $x$ is less than $1.1605$, so at
+% most $0.161$ is fed to the power series.
% \begin{macrocode}
-\def\MFP@numdigits#1{\@xp\MFP@ninthofmany#176543210\mfp@end}%
-\def\MFP@numzeros#1{%
- \@xp\MFP@ninthofmany#112345678\mfp@end}%
\def\MFP@Rlog@reduce{%
- \MFP@reduceonce 3{16227766}{31622777}{50000000}%
- \MFP@reduceonce 1{77827941}{56234133}{25000000}%
- \MFP@reduceonce 1{33352143}{74989421}{12500000}%
- \MFP@reduceonce 1{15478198}{86596433}{06250000}}%
+ \MFP@Rcopy tx\MFP@RmulC\MFP@Rcopyz t%
+ \MFP@reduceonce {316}{31622776}{60168379}{50000000}%
+ \MFP@reduceonce {177}{56234132}{51903491}{25000000}%
+ \MFP@reduceonce {133}{74989420}{93324558}{12500000}%
+ \MFP@reduceonce {115}{86596432}{33600654}{06250000}%
+ \MFP@Rcopy tx\MFP@RdivC\MFP@Rcopyz t}%
\def\MFP@reduceonce#1#2#3#4{%
- \MFP@Rcopy tx\MFP@Rload y1{#1}{#2}\MFP@Rcmp
- \ifMFP@neg\else
- \MFP@Rload y10{#3}\MFP@Rmul
- \MFP@Rcopyz t\MFP@Rcopy sx\MFP@Rload y10{#4}\MFP@Radd
+ \ifnum\MFP@t@Int>#1\relax
+ \MFP@Rcopy tx%
+ \MFP@DPmul 0{#2}{#3}\MFP@Rcopyz t%
+ \MFP@Rcopy sx\MFP@Rload y10{#4}\MFP@Radd
\MFP@Rcopyz s%
\fi}%
% \end{macrocode}
%
-% Now we have a value for $t$ of the form $1 + u$ with $0\le u < 0.155$
-% (approx.) We will use the formula
+% Now we have a value for $t$ of the form $1 + u$ with $0\le u < 0.161$.
+% We will use the formula
% $$
-% \log (1 + u) = \frac{1}{\ln 10} \sum_{n=0}^\infty (-1)^n \frac{u^{n+1}}{n+1}.
+% \ln (1 + u) = \sum_{n=1}^\infty (-1)^{n-1} \frac{u^n}{n}.
% $$
-% We only need to carry it far enough to assure that the next term would
-% be $0$ in our finite resolution arithmetic, that is $.155^{k}/k/\ln10 < .5\times
-% 10^{-8}$. This is satisfied by $k=9$, so we only need eight terms.
+% We only need to carry it far enough to assure that the remaining terms
+% would be zero in our finite resolution arithmetic, that is
+% $(.161)^k/k < .5\times 10^{-8}$. This is satisfied by $k=10$.
+% So we carry the sum to 9 places.
%
% Again, we compute this by
% $$
-% ru(1 - au(1 - bu(1-cu(1-du(1-eu(1-fu(1-gu(1-hu))))))))
+% u(1-au(1-bu(1-cu(1-du(1-eu(1-fu(1-gu(1-hu))))))))
% $$
-% where $a= 1/2$, $b = 2/3$,\dots,$h=7/8/2$, and $r = 1/ln 10 = \log e$.
-% This allows us to reuse \cs{MFP@com@iter}.
-%
-% The for the natural log $\ln x$ we multiply the common logarithm by
-% $\ln 10$. It seems we are redundantly canceling the earlier division
-% by $\ln 10$, but remember we are also multiplying the integer part,
-% which did not have such a division in its code.
+% where $a= 1/2$, $b = 2/3$,\dots, $g=7/8$, and $h=8/9$
+% This arrangement allows us to reuse \cs{MFP@com@iter}.
% \begin{macrocode}
\def\MFP@Rlog@prog{%
- \MFP@Rcopy tx\MFP@Rdecr
- \MFP@Rcopyz t%
- \MFP@Rcopyzx\MFP@Rload y10{43750000}\MFP@Rsmul
- \MFP@com@iter{85714286}\MFP@com@iter{83333333}\MFP@com@iter{80000000}%
- \MFP@com@iter{75000000}\MFP@com@iter{66666667}\MFP@com@iter{50000000}%
- \MFP@com@iter{43429448}\MFP@Rcopyzx\MFP@Rcopy sy\MFP@Radd}%
-\def\MFP@Rln{%
- \MFP@Rlog\MFP@Rcopyzx\MFP@Rload y12{30258509}\MFP@Rmul}%
+ \MFP@Rcopy tx\MFP@Rload y10{88888889}\MFP@Rsmul
+ \MFP@com@iter{87500000}\MFP@com@iter{85714286}\MFP@com@iter{83333333}%
+ \MFP@com@iter{80000000}\MFP@com@iter{75000000}\MFP@com@iter{66666667}%
+ \MFP@com@iter{50000000}\MFP@flipz\MFP@Rcopyzx\MFP@Rcopy ty\MFP@Rsmul}%
% \end{macrocode}
%
% \subsection{Powers}
@@ -3472,10 +3517,10 @@
% \begin{macrocode}
\def\MFP@Rexp@pos{%
\MFP@Rload y1{18}{42068074}\MFP@Rcmp
- \ifMFP@pos % overflow
+ \ifMFP@pos
\MFP@expoverflow@err
\MFP@Rloadz 1\MaxRealInt\MaxRealFrac
- \else % handle integer part
+ \else
\MFP@tempa\MFP@x@Int
\edef\MFP@powerof@e{%
1\ifcase\MFP@tempa
@@ -3503,7 +3548,7 @@
\ifnum\MFP@x@Frc=0
\else
\MFP@Rcopyz s%
- \MFP@tempa=\@xp\MFP@firstofmany\MFP@x@Frc\mfp@end
+ \MFP@tempa=\@xp\MFP@oneofmany\MFP@x@Frc\MFP@end
\edef\MFP@powerof@e{%
y1\ifcase\MFP@tempa
10\or
@@ -3520,7 +3565,7 @@
\edef\MFP@t@Frc{0\@xp\@gobble\MFP@x@Frc}%
\MFP@Rcopy sx\@xp\MFP@Rload\MFP@powerof@e\MFP@Rmul
\ifnum\MFP@t@Frc=0
- \else % handle the rest
+ \else
\MFP@Rcopyz s\MFP@Rload t10\MFP@t@Frc
\MFP@Rexp@pos@prog
\MFP@Rcopy sx\MFP@Rcopyzy\MFP@Rmul
@@ -3537,16 +3582,17 @@
% $$
% 1 + (x + x/2(x + x/3(x + x/4(x + x/5(x + x/6)))))
% $$
-% We start by loading $x$ (now in register \texttt{t}) into register
-% \texttt{z}, then repeatedly run \cs{MFP@Rexp@iter} feeding it the
+% We start by loading $x$ (now in register $t$) into register
+% $z$, then repeatedly run \cs{MFP@Rexp@iter} feeding it the
% successive values of $1/n$. This iterator first multiplies the most
-% recent result (the \texttt{z} register) by $1/n$, then that by $x$ and
+% recent result (the $z$ register) by $1/n$, then that by $x$ and
% then adds $x$ to that. The final step is to add $1$.
% \begin{macrocode}
\def\MFP@Rexp@pos@prog{%
- \MFP@Rcopy tz\MFP@Rexp@iter{16666667}\MFP@Rexp@iter{20000000}%
- \MFP@Rexp@iter{25000000}\MFP@Rexp@iter{33333333}%
- \MFP@Rexp@iter{50000000}\MFP@Rcopyzx\MFP@Rincr}%
+ \MFP@Rcopy tz\MFP@Rexp@iter{14285714}\MFP@Rexp@iter{16666667}%
+ \MFP@Rexp@iter{20000000}\MFP@Rexp@iter{25000000}%
+ \MFP@Rexp@iter{33333333}\MFP@Rexp@iter{50000000}\MFP@Rcopyzx
+ \MFP@Rincr}%
\def\MFP@Rexp@iter#1{%
\MFP@Rcopyzx\MFP@Rload y10{#1}\MFP@Rsmul
\MFP@Rcopyzx\MFP@Rcopy ty\MFP@Rsmul
@@ -3556,16 +3602,14 @@
% This is because an absolute error in $x$ converts to a relative error
% in $e^x$, That is, knowing $x$ only to $10^{-8}$ means $e^x$ is off by
% (about) $e^x\cdot 10^{-8}$. Roughly speaking, this means only about $8$
-% places of $e^x$ are accurate, so if the integer part of $e^x$ has six
-% places then only two places after the decimal are significant. Even if
-% $x$ is exact, we can only represent $e$ itself to eight decimals and the
-% repeated multiplications accumulate errors in such a way that one still
-% cannot get more than eight significant figures.
+% places of $e^x$ are accurate, so if (for example) the integer part of
+% $e^x$ has six places then only two places after the decimal are
+% significant.
%
% \bigskip
% The first issue with negative exponents is that it doesn't take much to
% produce a value of $e^{-x}$ that rounds to $0$. Any $x > 19.11382792$. So
-% we start by comparing to that value and simply return $0$ if $x$ is
+% we start by comparing to that value and simply return zero if $x$ is
% larger.
%
% We perform exactly the same reductions as for positive exponents,
@@ -3577,7 +3621,7 @@
\MFP@Rcmp
\ifMFP@pos
\MFP@Rloadz 000%
- \else % handle integer part
+ \else
\MFP@tempa\MFP@x@Int
\edef\MFP@powerof@e{%
\ifcase\MFP@tempa
@@ -3604,9 +3648,9 @@
000\fi}%
\@xp\MFP@Rloadz\MFP@powerof@e
\ifnum\MFP@x@Frc=0
- \else % handle first decimal digit
+ \else
\MFP@Rcopyz s%
- \MFP@tempa=\@xp\MFP@firstofmany\MFP@x@Frc\mfp@end
+ \MFP@tempa=\@xp\MFP@oneofmany\MFP@x@Frc\MFP@end
\edef\MFP@powerof@e{%
y1\ifcase\MFP@tempa
10\or
@@ -3623,10 +3667,10 @@
\edef\MFP@t@Frc{0\@xp\@gobble\MFP@x@Frc}%
\MFP@Rcopy sx\@xp\MFP@Rload\MFP@powerof@e\MFP@Rmul
\ifnum\MFP@t@Frc=0
- \else % handle the rest
+ \else
\MFP@Rcopyz s\MFP@Rload t10\MFP@t@Frc
\MFP@Rexp@neg@prog
- \MFP@Rcopy sx\MFP@Rcopyzy\MFP@Rmul
+ \MFP@Rcopyzx\MFP@Rcopy sy\MFP@Rmul
\fi
\fi
\fi}%
@@ -3643,10 +3687,11 @@
% rather than call \cs{MFP@com@iter} with a useless multiplication by $1$.
% \begin{macrocode}
\def\MFP@Rexp@neg@prog{%
- \MFP@Rcopy tx\MFP@Rload y10{16666667}\MFP@Rsmul
- \MFP@com@iter{20000000}\MFP@com@iter{25000000}%
- \MFP@com@iter{33333333}\MFP@com@iter{50000000}%
- \MFP@flipz\MFP@Rcopyzx\MFP@Rcopy ty\MFP@Rsmul\MFP@flipz}%
+ \MFP@Rcopy tx\MFP@Rload y10{14285712}\MFP@Rsmul
+ \MFP@com@iter{16666667}\MFP@com@iter{20000000}%
+ \MFP@com@iter{25000000}\MFP@com@iter{33333333}%
+ \MFP@com@iter{50000000}\MFP@flipz\MFP@Rcopyzx
+ \MFP@Rcopy ty\MFP@Rsmul\MFP@flipz}%
% \end{macrocode}
%
% The most efficient way to take an integer power of a number $x$ is to
@@ -3675,11 +3720,14 @@
% The above scheme requires at most $\lfloor\log_2 n\rfloor$ squarings
% and at most $\lceil \log_2 n \rceil$ multiplications for $x^n$, while
% directly multiplying $x\cdot x \cdots x$ would require $n-1$
-% multiplications.
+% multiplications. I have tested with an exponent equal to $8000$ and it
+% takes only about $25$ times as long as a single multiplication (rather
+% than $7999$ times).
%
-% For negative powers we can either find the positive power and take its
-% reciprocal or take the reciprocal of $x$ and find its positive power. We
-% do the first so that overflow can be detected in \cs{MFP@@Rpow}.
+% For negative powers we can either find the positive power of $x$ and
+% take its reciprocal or take the reciprocal of $x$ and find its positive
+% power. We do the second so that overflow can be detected in
+% \cs{MFP@@Rpow}.
% \begin{macrocode}
\def\MFP@Rpow{%
\ifnum\MFP@y@Frc>0
@@ -3697,34 +3745,27 @@
\MFP@badpower@err
\MFP@Rloadz 1\xOverZeroInt\xOverZeroFrac
\fi
- \else % integer power of nonzero number
- % get the sign in case an overflow interrupts the calculations.
+ \else
\ifnum\MFP@x@Sgn>0
\def\MFP@power@Sgn{1}%
\else
\edef\MFP@power@Sgn{\ifodd\MFP@loopctr -\fi 1}%
-\ifMFPdebug
- \show\MFP@power@Sgn
-\fi
\fi
\ifnum\MFP@y@Sgn<0 \MFP@Rinv \MFP@Rcopyzx\fi
\ifnum\MFP@loopctr=1
\MFP@Rloadz \MFP@power@Sgn\MFP@x@Int\MFP@x@Frc
\else
\MFP@@Rpow
-\ifMFPdebug
- \showreg z%
-\fi
\fi
\fi
- \fi}% %
+ \fi}%
% \end{macrocode}
%
% This implements the algorithm discussed above. We save $x$ in register
-% \texttt{q}, initialize the starting value of $1$ in \texttt{p} and then
+% $q$, initialize the starting value of $1$ in \reg{p} and then
% run the loop. If the binary digit just read is a 1 (i.e., \cs{ifodd} is
% true), it multiplies $p$ and $q$. It also saves the last product (copies
-% \texttt{z} to \texttt{p}). This need not be done on the last iteration,
+% \reg{z} to \reg{p}). This need not be done on the last iteration,
% but must not be moved out of the \cs{ifodd} conditional because
% intervening computations modify $z$. If there are more iterations to do
% (i.e., the \cs{ifnum} is true), this squares $q$ and reduces the
@@ -3735,8 +3776,8 @@
% break the loop and return $\pm\infty$.
% \begin{macrocode}
\def\MFP@@Rpow{%
- \MFP@Rcopy xq% initialize register to be squared
- \MFP@Rload p110% initialize register to hold partial products
+ \MFP@Rcopy xq%
+ \MFP@Rload p110%
\MFP@Rpow@loop}%
\def\MFP@Rpow@loop{%
\ifodd\MFP@loopctr
@@ -3758,13 +3799,12 @@
\def\MFP@handle@expoverflow{%
\MFP@expoverflow@err
\MFP@loopctr=0
- \MFP@Rloadz\MFP@power@Sgn\MaxRealInt\MaxRealFrac
-}%
+ \MFP@Rloadz\MFP@power@Sgn\MaxRealInt\MaxRealFrac}%
% \end{macrocode}
%
% \subsection{The square root}
%
-% One can combine logarithms and exponentials to can get any power: to get
+% One can combine logarithms and exponentials to get any power: to get
% $x^y$, compute $e^{y\ln x}$. This has the disadvantage that it doesn't
% work if $x$ is negative. Most powers of negative numbers are not
% defined, but certainly integer powers are. Thus we have defined
@@ -3778,11 +3818,15 @@
% $\sqrt{9} = 3$. In fact, if a square root can be expressed exactly
% within our 8-digit precision, our code will find it.
%
-% For the square root we return $0$ if $x$ is not positive. If the integer
+% For the square root we return zero if $x$ is not positive. If the integer
% part of $x$ is $0$, we copy the fractional part to the integer part
% (that is, we multiply by $10^{8}$, remembering to multiply by $10^{-4}$
-% later). This makes the square root of such numbers slightly more
-% accurate. We then compute the square root using an algorithm that will
+% later). This makes the square root of such numbers rather more
+% accurate. (To get around some other rare but annoying inaccuracies, we
+% go through a similar process when the integer part of $x$ is at most $4$
+% digits, multiplying by $10^4$ before and by $10^{-2}$ after.)
+%
+% We then compute the square root using an algorithm that will
% be exact whenever possible. We perform one additional processing step.
% To explain it, note that our algorithm actually produces the largest
% number $s$ with four digits right of the decimal place that satisfies $s^2
@@ -3801,46 +3845,46 @@
% With this value, $s + \bar\epsilon$ misses the exact square root by at
% most $\epsilon^2/(2s) < .5\cdot 10^{-8}$, because $s \ge 1$.
% The final result $s + \bar\epsilon$ is equivalent to computing the
-% average $s$ and $x/s$. This, possibly divided by $10^4$ is the
+% average $s$ and $x/s$. This, possibly divided by $10^4$ or $10^2$ is the
% returned value.
%
-% I originally tried power series methods, but they failed to produce
-% exact answers when they existed (unless they were inconveniently carried
-% to nine decimals and then rounded to eight). Then I tried the ``exact when
-% possible'' algorithm to get $s$, but correcting it as follows: find
-% $\sqrt{x/s^2}$ by power series and multiply by $s$. But this turned out
-% to be remarkably inaccurate, being paradoxically worst when $s$ is
-% already close, but not equal, to $\sqrt{x}$. Finally, I hit on the above
-% simple and efficient correction, which turns out to be also the most
-% accurate. By tests it produces a result correct in all but the last
-% (eighth) decimal place, and within $1$ of the correct value in that
-% place. Of course, it still produces exact results when that is possible
-% within our digit limits.
+% By tests, with rare exceptions, our computations produces a result
+% correct in all eight decimal places. In the rare exception, the last
+% place is within $1$ of the correct value.
% \begin{macrocode}
\def\MFP@Rsqrt{%
\ifcase\MFP@x@Sgn\relax
\MFP@Rzero
\or
\ifnum\MFP@x@Int=0
- \let\MFP@sqrt@reduce=Y%
+ \def\MFP@sqrt@reduce{2}%
\edef\MFP@x@Int{\number\MFP@x@Frc}%
\edef\MFP@x@Frc{00000000}%
+ \else\ifnum\MFP@x@Int<10000
+ \def\MFP@sqrt@reduce{1}%
+ \edef\MFP@x@Int{\MFP@x@Int\@xp\MFP@fourofmany\MFP@x@Frc\MFP@end}%
+ \edef\MFP@x@Frc{\@xp\MFP@gobblefour\MFP@x@Frc0000}%
\else
- \let\MFP@sqrt@reduce=N%
- \fi
+ \def\MFP@sqrt@reduce{0}%
+ \fi\fi
\MFP@Rcopy xt%
\MFP@Isqrt
\MFP@Rcopyz s\MFP@Rcopyzy
\MFP@Rcopy tx\MFP@Rdiv
\MFP@Rcopy sx\MFP@Rcopyzy\MFP@Radd
\MFP@Rcopyzx\MFP@Rhalve
- \ifx Y\MFP@sqrt@reduce
- \MFP@Rcopyzx\MFP@Rload y10{10000}\MFP@Rmul
+ \ifcase \MFP@sqrt@reduce\relax
+ \or
+ \MFP@Rcopyzx\MFP@Rload y10{01000000}\MFP@Rmul
+ \or
+ \MFP@Rcopyzx\MFP@Rload y10{00010000}\MFP@Rmul
\fi
\else
- \MFP@sqrtofneg@err
+ \MFP@warn{Square root of a negative number. Zero will be returned.}%
\MFP@Rzero
\fi}%
+\def\MFP@fourofmany#1#2#3#4#5\MFP@end{#1#2#3#4}%
+\def\MFP@gobblefour#1#2#3#4{}%
% \end{macrocode}
%
% There is a rather straightforward pencil and paper algorithm that
@@ -3907,7 +3951,7 @@
\def\MFP@Isqrt{%
\MFP@ItoQ\MFP@x@Int\MFP@x@Frc
\MFP@tempa=0 \MFP@tempb=0 \MFP@tempc=0
- \expandafter\MFP@Isqrt@loop\MFP@ItoQ@Tmp\mfp@end
+ \expandafter\MFP@Isqrt@loop\MFP@ItoQ@Tmp\MFP@end
\MFP@tempa=\MFP@tempc
\divide\MFP@tempc\MFP@tttfour
\edef\MFP@z@Int{\number\MFP@tempc}%
@@ -3925,7 +3969,7 @@
% it written out, I am surprise at how concise and elegant it is!
% \begin{macrocode}
\def\MFP@Isqrt@loop#1{%
- \ifx\mfp@end #1%
+ \ifx\MFP@end #1%
\else
\multiply\MFP@tempa 2 \multiply\MFP@tempb 4 \multiply\MFP@tempc 2
\advance \MFP@tempb#1\relax
@@ -3935,18 +3979,18 @@
\advance\MFP@tempa 1
\fi
\expandafter\MFP@Isqrt@loop
- \fi
-}%
+ \fi}%
+\MFP@xfinish
%</extra>
% \end{macrocode}
-% For my own benefit: this finds the next binary digit and updates the
-% square root (in \cs{MFP@tempc}) by appending that digit. The new digit
-% is also appended to the end of \cs{MFP@tempa}. This is subtracted from
-% \cs{MFP@tempb}, but only if the last digit is a 1. Then the next
-% quadrenary digit is appended to \cs{MFP@tempb}. Finally, the last binary
-% digit found is added (not appended) to \cs{MFP@tempa}. The ``appending''
-% of a digit means a multiplication by $2$ (or $4$) and the addition of the
-% digit. We perform such additions only if the digit is a 1, and we
-% determine if the digit is 1 or 0 by the \cs{ifnum} test.
+% For my own benefit: the above code finds the next binary digit and
+% updates the square root (in \cs{MFP@tempc}) by appending that digit. The
+% new digit is also appended to the end of \cs{MFP@tempa}. This is
+% subtracted from \cs{MFP@tempb}, but only if the last digit is a 1. Then
+% the next quadrenary digit is appended to \cs{MFP@tempb}. Finally, the
+% last binary digit found is added (not appended) to \cs{MFP@tempa}. The
+% ``appending'' of a digit means a multiplication by $2$ (or $4$) and the
+% addition of the digit. We perform such additions only if the digit is a
+% 1, and we determine if the digit is 1 or 0 by the \cs{ifnum} test.
%\Finale
%