diff options
Diffstat (limited to 'macros/latex/contrib/l3kernel/l3fp-extended.dtx')
-rw-r--r-- | macros/latex/contrib/l3kernel/l3fp-extended.dtx | 1275 |
1 files changed, 0 insertions, 1275 deletions
diff --git a/macros/latex/contrib/l3kernel/l3fp-extended.dtx b/macros/latex/contrib/l3kernel/l3fp-extended.dtx deleted file mode 100644 index bcd16a4b68..0000000000 --- a/macros/latex/contrib/l3kernel/l3fp-extended.dtx +++ /dev/null @@ -1,1275 +0,0 @@ -% \iffalse meta-comment -% -%% File: l3fp-extended.dtx -% -% Copyright (C) 2011-2024 The LaTeX Project -% -% It may be distributed and/or modified under the conditions of the -% LaTeX Project Public License (LPPL), either version 1.3c of this -% license or (at your option) any later version. The latest version -% of this license is in the file -% -% https://www.latex-project.org/lppl.txt -% -% This file is part of the "l3kernel bundle" (The Work in LPPL) -% and all files in that bundle must be distributed together. -% -% ----------------------------------------------------------------------- -% -% The development version of the bundle can be found at -% -% https://github.com/latex3/latex3 -% -% for those people who are interested. -% -%<*driver> -\documentclass[full,kernel]{l3doc} -\begin{document} - \DocInput{\jobname.dtx} -\end{document} -%</driver> -% \fi -% -% \title{^^A -% The \pkg{l3fp-extended} module\\ -% Manipulating numbers with extended precision, for internal use^^A -% } -% \author{^^A -% The \LaTeX{} Project\thanks -% {^^A -% E-mail: -% \href{mailto:latex-team@latex-project.org} -% {latex-team@latex-project.org}^^A -% }^^A -% } -% \date{Released 2024-03-14} -% -% \maketitle -% -% \begin{documentation} -% -% \end{documentation} -% -% \begin{implementation} -% -% \section{\pkg{l3fp-extended} implementation} -% -% \begin{macrocode} -%<*package> -% \end{macrocode} -% -% \begin{macrocode} -%<@@=fp> -% \end{macrocode} -% -% \subsection{Description of fixed point numbers} -% -% This module provides a few functions to manipulate positive floating -% point numbers with extended precision ($24$ digits), but mostly -% provides functions for fixed-point numbers with this precision ($24$ -% digits). Those are used in the computation of -% Taylor series for the logarithm, exponential, and trigonometric -% functions. Since we eventually only care about the $16$ first digits -% of the final result, some of the calculations are not performed with -% the full $24$-digit precision. In other words, the last two blocks of -% each fixed point number may be wrong as long as the error is small -% enough to be rounded away when converting back to a floating point -% number. The fixed point numbers are expressed as -% \begin{quote} -% \Arg{a_1} \Arg{a_2} \Arg{a_3} \Arg{a_4} \Arg{a_5} \Arg{a_6} |;| -% \end{quote} -% where each \meta{a_i} is exactly $4$ digits (ranging from |0000| to -% |9999|), except \meta{a_1}, which may be any \enquote{not-too-large} -% non-negative integer, with or without leading zeros. Here, -% \enquote{not-too-large} depends on the specific function (see the -% corresponding comments for details). Checking for overflow is the -% responsibility of the code calling those functions. The fixed point -% number $a$ corresponding to the representation above is $a = -% \sum_{i=1}^{6} \meta{a_i} \cdot 10^{-4i}$. -% -% Most functions we define here have the form -% \begin{syntax} -% \cs{@@_fixed_\meta{calculation}:wwn} \meta{operand_1} |;| \meta{operand_2} |;| \Arg{continuation} -% \end{syntax} -% They perform the \meta{calculation} on the two \meta{operands}, then -% feed the result ($6$ brace groups followed by a semicolon) to the -% \meta{continuation}, responsible for the next step of the calculation. -% Some functions only accept an \texttt{N}-type \meta{continuation}. -% This allows constructions such as -% \begin{quote} -% \cs{@@_fixed_add:wwn} \meta{X_1} |;| \meta{X_2} |;| \\ -% \cs{@@_fixed_mul:wwn} \meta{X_3} |;| \\ -% \cs{@@_fixed_add:wwn} \meta{X_4} |;| \\ -% \end{quote} -% to compute $(X_1+X_2)\cdot X_3 + X_4$. This turns out to be very -% appropriate for computing continued fractions and Taylor series. -% -% At the end of the calculation, the result is turned back to a floating -% point number using \cs{@@_fixed_to_float_o:wN}. This function has to -% change the exponent of the floating point number: it must be used -% after starting an integer expression for the overall exponent of the -% result. -% -% \subsection{Helpers for numbers with extended precision} -% -% \begin{variable}{\c_@@_one_fixed_tl} -% The fixed-point number~$1$, used in \pkg{l3fp-expo}. -% \begin{macrocode} -\tl_const:Nn \c_@@_one_fixed_tl - { {10000} {0000} {0000} {0000} {0000} {0000} ; } -% \end{macrocode} -% \end{variable} -% -% \begin{macro}[EXP]{\@@_fixed_continue:wn} -% This function simply calls the next function. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_continue:wn #1; #2 { #2 #1; } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_fixed_add_one:wN} -% \begin{syntax} -% \cs{@@_fixed_add_one:wN} \meta{a} |;| \meta{continuation} -% \end{syntax} -% This function adds $1$ to the fixed point \meta{a}, by changing -% $a_1$ to $10000+a_1$, then calls the \meta{continuation}. This -% requires $a_1 + 10000 < 2^{31}$. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_add_one:wN #1#2; #3 - { - \exp_after:wN #3 \exp_after:wN - { \int_value:w \@@_int_eval:w \c_@@_myriad_int + #1 } #2 ; - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_fixed_div_myriad:wn} -% Divide a fixed point number by $10000$. This is a little bit more -% subtle than just removing the last group and adding a leading group -% of zeros: the first group~|#1| may have any number of digits, and we -% must split~|#1| into the new first group and a second group of -% exactly $4$~digits. The choice of shifts allows~|#1| to be in the -% range $[0, 5\cdot 10^{8}-1]$. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_div_myriad:wn #1#2#3#4#5#6; - { - \exp_after:wN \@@_fixed_mul_after:wwn - \int_value:w \@@_int_eval:w \c_@@_leading_shift_int - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_trailing_shift_int - + #1 ; {#2}{#3}{#4}{#5}; - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_fixed_mul_after:wwn} -% The fixed point operations which involve multiplication end by -% calling this auxiliary. It braces the last block of digits, and -% places the \meta{continuation} |#3| in front. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_mul_after:wwn #1; #2; #3 { #3 {#1} #2; } -% \end{macrocode} -% \end{macro} -% -% \subsection{Multiplying a fixed point number by a short one} -% -% \begin{macro}[EXP]{\@@_fixed_mul_short:wwn} -% \begin{syntax}\parskip=0pt\obeylines -% \cs{@@_fixed_mul_short:wwn} -% | |\Arg{a_1} \Arg{a_2} \Arg{a_3} \Arg{a_4} \Arg{a_5} \Arg{a_6} |;| -% | |\Arg{b_0} \Arg{b_1} \Arg{b_2} |;| \Arg{continuation} -% \end{syntax} -% Computes the product $c=ab$ of $a=\sum_i \meta{a_i} 10^{-4i}$ and -% $b=\sum_i \meta{b_i} 10^{-4i}$, rounds it to the closest multiple of -% $10^{-24}$, and leaves \meta{continuation} \Arg{c_1} \ldots{} -% \Arg{c_6} |;| in the input stream, where each of the \meta{c_i} are -% blocks of $4$~digits, except \meta{c_1}, which is any \TeX{} -% integer. Note that indices for \meta{b} start at~$0$: for instance -% a second operand of |{0001}{0000}{0000}| leaves the first operand -% unchanged (rather than dividing it by $10^{4}$, as -% \cs{@@_fixed_mul:wwn} would). -% \begin{macrocode} -\cs_new:Npn \@@_fixed_mul_short:wwn #1#2#3#4#5#6; #7#8#9; - { - \exp_after:wN \@@_fixed_mul_after:wwn - \int_value:w \@@_int_eval:w \c_@@_leading_shift_int - + #1*#7 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #1*#8 + #2*#7 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #1*#9 + #2*#8 + #3*#7 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #2*#9 + #3*#8 + #4*#7 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #3*#9 + #4*#8 + #5*#7 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_trailing_shift_int - + #4*#9 + #5*#8 + #6*#7 - + ( #5*#9 + #6*#8 + #6*#9 / \c_@@_myriad_int ) - / \c_@@_myriad_int ; ; - } -% \end{macrocode} -% \end{macro} -% -% \subsection{Dividing a fixed point number by a small integer} -% -% \begin{macro}[EXP]{\@@_fixed_div_int:wwN} -% \begin{macro}[EXP] -% { -% \@@_fixed_div_int:wnN, \@@_fixed_div_int_auxi:wnn, -% \@@_fixed_div_int_auxii:wnn, \@@_fixed_div_int_pack:Nw, -% \@@_fixed_div_int_after:Nw -% } -% \begin{syntax} -% \cs{@@_fixed_div_int:wwN} \meta{a} |;| \meta{n} |;| \meta{continuation} -% \end{syntax} -% Divides the fixed point number \meta{a} by the (small) integer -% $0<\meta{n}<10^4$ and feeds the result to the \meta{continuation}. -% There is no bound on $a_1$. -% -% The arguments of the \texttt{i} auxiliary are 1: one of the $a_{i}$, -% 2: $n$, 3: the \texttt{ii} or the \texttt{iii} auxiliary. It -% computes a (somewhat tight) lower bound $Q_{i}$ for the ratio -% $a_{i}/n$. -% -% The \texttt{ii} auxiliary receives $Q_{i}$, $n$, and $a_{i}$ as -% arguments. It adds $Q_{i}$ to a surrounding integer expression, and -% starts a new one with the initial value $9999$, which ensures that -% the result of this expression has $5$ digits. The auxiliary -% also computes $a_{i}-n\cdot Q_{i}$, placing the result in front of -% the $4$ digits of $a_{i+1}$. The resulting $a'_{i+1} = 10^{4} -% (a_{i} - n \cdot Q_{i}) + a_{i+1}$ serves as the first argument for -% a new call to the \texttt{i} auxiliary. -% -% When the \texttt{iii} auxiliary is called, the situation looks like -% this: -% \begin{quote} -% \cs{@@_fixed_div_int_after:Nw} \meta{continuation} \\ -% $-1 + Q_{1}$ \\ -% \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{2}$ \\ -% \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{3}$ \\ -% \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{4}$ \\ -% \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{5}$ \\ -% \cs{@@_fixed_div_int_pack:Nw} $9999$ \\ -% \cs{@@_fixed_div_int_auxii:wnn} $Q_{6}$ |;| \Arg{n} \Arg{a_{6}} -% \end{quote} -% where expansion is happening from the last line up. The -% \texttt{iii} auxiliary adds $Q_{6} + 2 \simeq a_{6}/n + 1$ to the -% last $9999$, giving the integer closest to $10000 + a_{6}/n$. -% -% Each \texttt{pack} auxiliary receives $5$ digits followed by a -% semicolon. The first digit is added as a carry to the integer -% expression above, and the $4$ other digits are braced. Each call to -% the \texttt{pack} auxiliary thus produces one brace group. The last -% brace group is produced by the \texttt{after} auxiliary, which -% places the \meta{continuation} as appropriate. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_div_int:wwN #1#2#3#4#5#6 ; #7 ; #8 - { - \exp_after:wN \@@_fixed_div_int_after:Nw - \exp_after:wN #8 - \int_value:w \@@_int_eval:w - 1 - \@@_fixed_div_int:wnN - #1; {#7} \@@_fixed_div_int_auxi:wnn - #2; {#7} \@@_fixed_div_int_auxi:wnn - #3; {#7} \@@_fixed_div_int_auxi:wnn - #4; {#7} \@@_fixed_div_int_auxi:wnn - #5; {#7} \@@_fixed_div_int_auxi:wnn - #6; {#7} \@@_fixed_div_int_auxii:wnn ; - } -\cs_new:Npn \@@_fixed_div_int:wnN #1; #2 #3 - { - \exp_after:wN #3 - \int_value:w \@@_int_eval:w #1 / #2 - 1 ; - {#2} - {#1} - } -\cs_new:Npn \@@_fixed_div_int_auxi:wnn #1; #2 #3 - { - + #1 - \exp_after:wN \@@_fixed_div_int_pack:Nw - \int_value:w \@@_int_eval:w 9999 - \exp_after:wN \@@_fixed_div_int:wnN - \int_value:w \@@_int_eval:w #3 - #1*#2 \@@_int_eval_end: - } -\cs_new:Npn \@@_fixed_div_int_auxii:wnn #1; #2 #3 { + #1 + 2 ; } -\cs_new:Npn \@@_fixed_div_int_pack:Nw #1 #2; { + #1; {#2} } -\cs_new:Npn \@@_fixed_div_int_after:Nw #1 #2; { #1 {#2} } -% \end{macrocode} -% \end{macro} -% \end{macro} -% -% \subsection{Adding and subtracting fixed points} -% -% \begin{macro}[EXP]{\@@_fixed_add:wwn, \@@_fixed_sub:wwn} -% \begin{macro}[EXP] -% { -% \@@_fixed_add:Nnnnnwnn, -% \@@_fixed_add:nnNnnnwn, -% \@@_fixed_add_pack:NNNNNwn, -% \@@_fixed_add_after:NNNNNwn -% } -% \begin{syntax} -% \cs{@@_fixed_add:wwn} \meta{a} |;| \meta{b} |;| \Arg{continuation} -% \end{syntax} -% Computes $a+b$ (resp.\ $a-b$) and feeds the result to the -% \meta{continuation}. This function requires $0\leq a_{1},b_{1}\leq -% 114748$, its result must be positive (this happens automatically for -% addition) and its first group must have at most~$5$ digits: $(a\pm -% b)_{1}<100000$. The two functions only differ by -% a sign, hence use a common auxiliary. It would be nice to grab the -% $12$ brace groups in one go; only $9$ parameters are allowed. Start -% by grabbing the sign, $a_{1}, \ldots, a_{4}$, the rest of $a$, -% and $b_{1}$ and $b_{2}$. The second auxiliary receives the rest of -% $a$, the sign multiplying $b$, the rest of $b$, and the -% \meta{continuation} as arguments. After going down through the -% various level, we go back up, packing digits and bringing the -% \meta{continuation} (|#8|, then |#7|) from the end of the argument -% list to its start. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_add:wwn { \@@_fixed_add:Nnnnnwnn + } -\cs_new:Npn \@@_fixed_sub:wwn { \@@_fixed_add:Nnnnnwnn - } -\cs_new:Npn \@@_fixed_add:Nnnnnwnn #1 #2#3#4#5 #6; #7#8 - { - \exp_after:wN \@@_fixed_add_after:NNNNNwn - \int_value:w \@@_int_eval:w 9 9999 9998 + #2#3 #1 #7#8 - \exp_after:wN \@@_fixed_add_pack:NNNNNwn - \int_value:w \@@_int_eval:w 1 9999 9998 + #4#5 - \@@_fixed_add:nnNnnnwn #6 #1 - } -\cs_new:Npn \@@_fixed_add:nnNnnnwn #1#2 #3 #4#5 #6#7 ; #8 - { - #3 #4#5 - \exp_after:wN \@@_fixed_add_pack:NNNNNwn - \int_value:w \@@_int_eval:w 2 0000 0000 #3 #6#7 + #1#2 ; {#8} ; - } -\cs_new:Npn \@@_fixed_add_pack:NNNNNwn #1 #2#3#4#5 #6; #7 - { + #1 ; {#7} {#2#3#4#5} {#6} } -\cs_new:Npn \@@_fixed_add_after:NNNNNwn 1 #1 #2#3#4#5 #6; #7 - { #7 {#1#2#3#4#5} {#6} } -% \end{macrocode} -% \end{macro} -% \end{macro} -% -% \subsection{Multiplying fixed points} -% -% ^^A todo: may a_1 or b_1 be = 10000? Used in ediv_epsi later. -% \begin{macro}[EXP]{\@@_fixed_mul:wwn} -% \begin{macro}[EXP]{\@@_fixed_mul:nnnnnnnw} -% \begin{syntax} -% \cs{@@_fixed_mul:wwn} \meta{a} |;| \meta{b} |;| \Arg{continuation} -% \end{syntax} -% Computes $a\times b$ and feeds the result to \meta{continuation}. -% This function requires $0\leq a_{1}, b_{1} < 10000$. Once more, we -% need to play around the limit of $9$ arguments for \TeX{} macros. -% Note that we don't need to obtain an exact rounding, contrarily to -% the |*| operator, so things could be harder. We wish to perform -% carries in -% \begin{align*} -% a \times b = -% & a_{1} \cdot b_{1} \cdot 10^{-8} \\ -% & + (a_{1} \cdot b_{2} + a_{2} \cdot b_{1}) \cdot 10^{-12} \\ -% & + (a_{1} \cdot b_{3} + a_{2} \cdot b_{2} -% + a_{3} \cdot b_{1}) \cdot 10^{-16} \\ -% & + (a_{1} \cdot b_{4} + a_{2} \cdot b_{3} -% + a_{3} \cdot b_{2} + a_{4} \cdot b_{1}) \cdot 10^{-20} \\ -% & + \Bigl(a_{2} \cdot b_{4} + a_{3} \cdot b_{3} + a_{4} \cdot b_{2} -% \\ & \qquad -% + \frac{a_{3} \cdot b_{4} + a_{4} \cdot b_{3} -% + a_{1} \cdot b_{6} + a_{2} \cdot b_{5} -% + a_{5} \cdot b_{2} + a_{6} \cdot b_{1}}{10^{4}} -% \\ & \qquad -% + a_{1} \cdot b_{5} + a_{5} \cdot b_{1}\Bigr) \cdot 10^{-24} -% + O(10^{-24}), -% \end{align*} -% where the $O(10^{-24})$ stands for terms which are at most $5\cdot -% 10^{-24}$; ignoring those leads to an error of at most -% $5$~\texttt{ulp}. Note how the first $15$~terms only depend on -% $a_{1},\ldots{},a_{4}$ and $b_{1},\ldots,b_{4}$, while the last -% $6$~terms only depend on $a_{1},a_{2},a_{5},a_{6}$, and the -% corresponding parts of~$b$. Hence, the first function grabs -% $a_{1},\ldots,a_{4}$, the rest of $a$, and $b_{1},\ldots,b_{4}$, and -% writes the $15$ first terms of the expression, including a left -% parenthesis for the fraction. The \texttt{i} auxiliary receives -% $a_{5}$, $a_{6}$, $b_{1}$, $b_{2}$, $a_{1}$, $a_{2}$, $b_{5}$, -% $b_{6}$ and finally the \meta{continuation} as arguments. It writes -% the end of the expression, including the right parenthesis and the -% denominator of the fraction. The \meta{continuation} -% is finally placed in front of the $6$ brace groups by -% \cs{@@_fixed_mul_after:wwn}. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_mul:wwn #1#2#3#4 #5; #6#7#8#9 - { - \exp_after:wN \@@_fixed_mul_after:wwn - \int_value:w \@@_int_eval:w \c_@@_leading_shift_int - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #1*#6 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #1*#7 + #2*#6 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #1*#8 + #2*#7 + #3*#6 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_middle_shift_int - + #1*#9 + #2*#8 + #3*#7 + #4*#6 - \exp_after:wN \@@_pack:NNNNNw - \int_value:w \@@_int_eval:w \c_@@_trailing_shift_int - + #2*#9 + #3*#8 + #4*#7 - + ( #3*#9 + #4*#8 - + \@@_fixed_mul:nnnnnnnw #5 {#6}{#7} {#1}{#2} - } -\cs_new:Npn \@@_fixed_mul:nnnnnnnw #1#2 #3#4 #5#6 #7#8 ; - { - #1*#4 + #2*#3 + #5*#8 + #6*#7 ) / \c_@@_myriad_int - + #1*#3 + #5*#7 ; ; - } -% \end{macrocode} -% \end{macro} -% \end{macro} -% -% \subsection{Combining product and sum of fixed points} -% -% \begin{macro}[EXP] -% { -% \@@_fixed_mul_add:wwwn, -% \@@_fixed_mul_sub_back:wwwn, -% \@@_fixed_mul_one_minus_mul:wwn, -% } -% \begin{syntax} -% \cs{@@_fixed_mul_add:wwwn} \meta{a} |;| \meta{b} |;| \meta{c} |;| \Arg{continuation} -% \cs{@@_fixed_mul_sub_back:wwwn} \meta{a} |;| \meta{b} |;| \meta{c} |;| \Arg{continuation} -% \cs{@@_fixed_one_minus_mul:wwn} \meta{a} |;| \meta{b} |;| \Arg{continuation} -% \end{syntax} -% Sometimes called |FMA| (fused multiply-add), these functions -% compute $a\times b + c$, $c - a\times b$, and $1 - a\times b$ and -% feed the result to the \meta{continuation}. Those functions require -% $0\leq a_{1}, b_{1}, c_{1} \leq 10000$. Since those functions are -% at the heart of the computation of Taylor expansions, we -% over-optimize them a bit, and in particular we do not factor out the -% common parts of the three functions. -% -% For definiteness, consider the task of computing $a\times b + c$. -% We perform carries in -% \begin{align*} -% a \times b + c = -% & (a_{1} \cdot b_{1} + c_{1} c_{2})\cdot 10^{-8} \\ -% & + (a_{1} \cdot b_{2} + a_{2} \cdot b_{1}) \cdot 10^{-12} \\ -% & + (a_{1} \cdot b_{3} + a_{2} \cdot b_{2} + a_{3} \cdot b_{1} -% + c_{3} c_{4}) \cdot 10^{-16} \\ -% & + (a_{1} \cdot b_{4} + a_{2} \cdot b_{3} + a_{3} \cdot b_{2} -% + a_{4} \cdot b_{1}) \cdot 10^{-20} \\ -% & + \Big(a_{2} \cdot b_{4} + a_{3} \cdot b_{3} + a_{4} \cdot b_{2} -% \\ & \qquad -% + \frac{a_{3} \cdot b_{4} + a_{4} \cdot b_{3} -% + a_{1} \cdot b_{6} + a_{2} \cdot b_{5} -% + a_{5} \cdot b_{2} + a_{6} \cdot b_{1}}{10^{4}} -% \\ & \qquad -% + a_{1} \cdot b_{5} + a_{5} \cdot b_{1} -% + c_{5} c_{6} \Big) \cdot 10^{-24} -% + O(10^{-24}), -% \end{align*} -% where $c_{1} c_{2}$, $c_{3} c_{4}$, $c_{5} c_{6}$ denote the -% $8$-digit number obtained by juxtaposing the two blocks of digits of -% $c$, and $\cdot$ denotes multiplication. The task is obviously -% tough because we have $18$ brace groups in front of us. -% -% Each of the three function starts the first two levels (the first, -% corresponding to $10^{-4}$, is empty), with $c_{1} c_{2}$ in the -% first level, calls the \texttt{i} auxiliary with arguments described -% later, and adds a trailing ${} + c_{5}c_{6}$ |;| -% \Arg{continuation}~|;|. The ${} + c_{5}c_{6}$ piece, which is -% omitted for \cs{@@_fixed_one_minus_mul:wwn}, is taken in the -% integer expression for the $10^{-24}$ level. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_mul_add:wwwn #1; #2; #3#4#5#6#7#8; - { - \exp_after:wN \@@_fixed_mul_after:wwn - \int_value:w \@@_int_eval:w \c_@@_big_leading_shift_int - \exp_after:wN \@@_pack_big:NNNNNNw - \int_value:w \@@_int_eval:w \c_@@_big_middle_shift_int + #3 #4 - \@@_fixed_mul_add:Nwnnnwnnn + - + #5 #6 ; #2 ; #1 ; #2 ; + - + #7 #8 ; ; - } -\cs_new:Npn \@@_fixed_mul_sub_back:wwwn #1; #2; #3#4#5#6#7#8; - { - \exp_after:wN \@@_fixed_mul_after:wwn - \int_value:w \@@_int_eval:w \c_@@_big_leading_shift_int - \exp_after:wN \@@_pack_big:NNNNNNw - \int_value:w \@@_int_eval:w \c_@@_big_middle_shift_int + #3 #4 - \@@_fixed_mul_add:Nwnnnwnnn - - + #5 #6 ; #2 ; #1 ; #2 ; - - + #7 #8 ; ; - } -\cs_new:Npn \@@_fixed_one_minus_mul:wwn #1; #2; - { - \exp_after:wN \@@_fixed_mul_after:wwn - \int_value:w \@@_int_eval:w \c_@@_big_leading_shift_int - \exp_after:wN \@@_pack_big:NNNNNNw - \int_value:w \@@_int_eval:w \c_@@_big_middle_shift_int + - 1 0000 0000 - \@@_fixed_mul_add:Nwnnnwnnn - - ; #2 ; #1 ; #2 ; - - ; ; - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_fixed_mul_add:Nwnnnwnnn} -% \begin{syntax} -% \cs{@@_fixed_mul_add:Nwnnnwnnn} \meta{op} |+| \meta{c_3} \meta{c_4} |;| -% ~~\meta{b} |;| \meta{a} |;| \meta{b} |;| \meta{op} -% ~~|+| \meta{c_5} \meta{c_6} |;| -% \end{syntax} -% Here, \meta{op} is either |+| or |-|. Arguments |#3|, |#4|, |#5| -% are \meta{b_1}, \meta{b_2}, \meta{b_3}; arguments |#7|, |#8|, |#9| -% are \meta{a_1}, \meta{a_2}, \meta{a_3}. We can build three levels: -% $a_{1} \cdot b_{1}$ for $10^{-8}$, $(a_{1} \cdot b_{2} + a_{2} \cdot -% b_{1})$ for $10^{-12}$, and $(a_{1} \cdot b_{3} + a_{2} \cdot b_{2} -% + a_{3} \cdot b_{1} + c_{3} c_{4})$ for $10^{-16}$. The $a$--$b$ -% products use the sign |#1|. Note that |#2| is empty for -% \cs{@@_fixed_one_minus_mul:wwn}. We call the \texttt{ii} auxiliary -% for levels $10^{-20}$ and $10^{-24}$, keeping the pieces of \meta{a} -% we've read, but not \meta{b}, since there is another copy later in -% the input stream. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_mul_add:Nwnnnwnnn #1 #2; #3#4#5#6; #7#8#9 - { - #1 #7*#3 - \exp_after:wN \@@_pack_big:NNNNNNw - \int_value:w \@@_int_eval:w \c_@@_big_middle_shift_int - #1 #7*#4 #1 #8*#3 - \exp_after:wN \@@_pack_big:NNNNNNw - \int_value:w \@@_int_eval:w \c_@@_big_middle_shift_int - #1 #7*#5 #1 #8*#4 #1 #9*#3 #2 - \exp_after:wN \@@_pack_big:NNNNNNw - \int_value:w \@@_int_eval:w \c_@@_big_middle_shift_int - #1 \@@_fixed_mul_add:nnnnwnnnn {#7}{#8}{#9} - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_fixed_mul_add:nnnnwnnnn} -% \begin{syntax} -% \cs{@@_fixed_mul_add:nnnnwnnnn} \meta{a} |;| \meta{b} |;| \meta{op} -% ~~|+| \meta{c_5} \meta{c_6} |;| -% \end{syntax} -% Level $10^{-20}$ is $(a_{1} \cdot b_{4} + a_{2} \cdot b_{3} + a_{3} -% \cdot b_{2} + a_{4} \cdot b_{1})$, multiplied by the sign, which was -% inserted by the \texttt{i} auxiliary. Then we prepare level -% $10^{-24}$. We don't have access to all parts of \meta{a} and -% \meta{b} needed to make all products. Instead, we prepare the -% partial expressions -% \begin{align*} -% & b_{1} + a_{4} \cdot b_{2} + a_{3} \cdot b_{3} + a_{2} \cdot b_{4} + a_{1} \\ -% & b_{2} + a_{4} \cdot b_{3} + a_{3} \cdot b_{4} + a_{2} . -% \end{align*} -% Obviously, those expressions make no mathematical sense: we -% complete them with $a_{5} \cdot {}$ and ${} \cdot b_{5}$, and with -% $a_{6} \cdot b_{1} + a_{5} \cdot {}$ and ${} \cdot b_{5} + a_{1} -% \cdot b_{6}$, and of course with the trailing ${} + c_{5} c_{6}$. -% To do all this, we keep $a_{1}$, $a_{5}$, $a_{6}$, and the -% corresponding pieces of \meta{b}. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_mul_add:nnnnwnnnn #1#2#3#4#5; #6#7#8#9 - { - ( #1*#9 + #2*#8 + #3*#7 + #4*#6 ) - \exp_after:wN \@@_pack_big:NNNNNNw - \int_value:w \@@_int_eval:w \c_@@_big_trailing_shift_int - \@@_fixed_mul_add:nnnnwnnwN - { #6 + #4*#7 + #3*#8 + #2*#9 + #1 } - { #7 + #4*#8 + #3*#9 + #2 } - {#1} #5; - {#6} - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_fixed_mul_add:nnnnwnnwN} -% \begin{syntax} -% \cs{@@_fixed_mul_add:nnnnwnnwN} \Arg{partial_1} \Arg{partial_2} -% ~~\Arg{a_1} \Arg{a_5} \Arg{a_6} |;| \Arg{b_1} \Arg{b_5} \Arg{b_6} |;| -% ~~\meta{op} |+| \meta{c_5} \meta{c_6} |;| -% \end{syntax} -% Complete the \meta{partial_1} and \meta{partial_2} expressions as -% explained for the \texttt{ii} auxiliary. The second one is divided -% by $10000$: this is the carry from level $10^{-28}$. The trailing -% ${} + c_{5} c_{6}$ is taken into the expression for level -% $10^{-24}$. Note that the total of level $10^{-24}$ is in the -% interval $[-5\cdot 10^{8}, 6\cdot 10^{8}$ (give or take a couple of -% $10000$), hence adding it to the shift gives a $10$-digit number, as -% expected by the packing auxiliaries. See \pkg{l3fp-aux} for the -% definition of the shifts and packing auxiliaries. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_mul_add:nnnnwnnwN #1#2 #3#4#5; #6#7#8; #9 - { - #9 (#4* #1 *#7) - #9 (#5*#6+#4* #2 *#7+#3*#8) / \c_@@_myriad_int - } -% \end{macrocode} -% \end{macro} -% -% \subsection{Extended-precision floating point numbers} -% -% In this section we manipulate floating point numbers with roughly $24$ -% significant figures (\enquote{extended-precision} numbers, in short, -% \enquote{ep}), which take the form of an integer exponent, followed by a -% comma, then six groups of digits, ending with a semicolon. The first -% group of digit may be any non-negative integer, while other groups of -% digits have $4$~digits. In other words, an extended-precision number -% is an exponent ending in a comma, then a fixed point number. The -% corresponding value is $0.\meta{digits}\cdot 10^{\meta{exponent}}$. -% This convention differs from floating points. -% -% \begin{macro}[EXP]{\@@_ep_to_fixed:wwn} -% \begin{macro}[EXP] -% {\@@_ep_to_fixed_auxi:www, \@@_ep_to_fixed_auxii:nnnnnnnwn} -% Converts an extended-precision number with an exponent at most~$4$ -% and a first block less than $10^{8}$ to a fixed point number whose -% first block has $12$~digits, hopefully starting with many zeros. -% \begin{macrocode} -\cs_new:Npn \@@_ep_to_fixed:wwn #1,#2 - { - \exp_after:wN \@@_ep_to_fixed_auxi:www - \int_value:w \@@_int_eval:w 1 0000 0000 + #2 \exp_after:wN ; - \exp:w \exp_end_continue_f:w - \prg_replicate:nn { 4 - \int_max:nn {#1} { -32 } } { 0 } ; - } -\cs_new:Npn \@@_ep_to_fixed_auxi:www 1#1; #2; #3#4#5#6#7; - { - \@@_pack_eight:wNNNNNNNN - \@@_pack_twice_four:wNNNNNNNN - \@@_pack_twice_four:wNNNNNNNN - \@@_pack_twice_four:wNNNNNNNN - \@@_ep_to_fixed_auxii:nnnnnnnwn ; - #2 #1#3#4#5#6#7 0000 ! - } -\cs_new:Npn \@@_ep_to_fixed_auxii:nnnnnnnwn #1#2#3#4#5#6#7; #8! #9 - { #9 {#1#2}{#3}{#4}{#5}{#6}{#7}; } -% \end{macrocode} -% \end{macro} -% \end{macro} -% -% ^^A todo: make it work when the arg is zero. -% \begin{macro}[EXP]{\@@_ep_to_ep:wwN} -% \begin{macro}[rEXP]{\@@_ep_to_ep_loop:N, \@@_ep_to_ep_end:www} -% \begin{macro}[EXP]{\@@_ep_to_ep_zero:ww} -% Normalize an extended-precision number. More precisely, leading -% zeros are removed from the mantissa of the argument, decreasing its -% exponent as appropriate. Then the digits are packed into $6$~groups -% of~$4$ (discarding any remaining digit, not rounding). Finally, the -% continuation~|#8| is placed before the resulting exponent--mantissa -% pair. The input exponent may in fact be given as an integer -% expression. The \texttt{loop} auxiliary grabs a digit: if it -% is~$0$, decrement the exponent and continue looping, and otherwise -% call the \texttt{end} auxiliary, which places all digits in the -% right order (the digit that was not~$0$, and any remaining digits), -% followed by some~$0$, then packs them up neatly in $3\times2=6$ -% blocks of four. At the end of the day, remove with \cs{@@_use_i:ww} -% any digit that did not make it in the final mantissa (typically only -% zeros, unless the original first block has more than~$4$ digits). -% \begin{macrocode} -\cs_new:Npn \@@_ep_to_ep:wwN #1,#2#3#4#5#6#7; #8 - { - \exp_after:wN #8 - \int_value:w \@@_int_eval:w #1 + 4 - \exp_after:wN \use_i:nn - \exp_after:wN \@@_ep_to_ep_loop:N - \int_value:w \@@_int_eval:w 1 0000 0000 + #2 \@@_int_eval_end: - #3#4#5#6#7 ; ; ! - } -\cs_new:Npn \@@_ep_to_ep_loop:N #1 - { - \if_meaning:w 0 #1 - - 1 - \else: - \@@_ep_to_ep_end:www #1 - \fi: - \@@_ep_to_ep_loop:N - } -\cs_new:Npn \@@_ep_to_ep_end:www - #1 \fi: \@@_ep_to_ep_loop:N #2; #3! - { - \fi: - \if_meaning:w ; #1 - - 2 * \c_@@_max_exponent_int - \@@_ep_to_ep_zero:ww - \fi: - \@@_pack_twice_four:wNNNNNNNN - \@@_pack_twice_four:wNNNNNNNN - \@@_pack_twice_four:wNNNNNNNN - \@@_use_i:ww , ; - #1 #2 0000 0000 0000 0000 0000 0000 ; - } -\cs_new:Npn \@@_ep_to_ep_zero:ww \fi: #1; #2; #3; - { \fi: , {1000}{0000}{0000}{0000}{0000}{0000} ; } -% \end{macrocode} -% \end{macro} -% \end{macro} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_ep_compare:wwww} -% \begin{macro}[EXP]{\@@_ep_compare_aux:wwww} -% In \pkg{l3fp-trig} we need to compare two extended-precision -% numbers. This is based on the same function for positive floating -% point numbers, with an extra test if comparing only $16$ decimals is -% not enough to distinguish the numbers. Note that this function only -% works if the numbers are normalized so that their first block is -% in~$[1000,9999]$. -% \begin{macrocode} -\cs_new:Npn \@@_ep_compare:wwww #1,#2#3#4#5#6#7; - { \@@_ep_compare_aux:wwww {#1}{#2}{#3}{#4}{#5}; #6#7; } -\cs_new:Npn \@@_ep_compare_aux:wwww #1;#2;#3,#4#5#6#7#8#9; - { - \if_case:w - \@@_compare_npos:nwnw #1; {#3}{#4}{#5}{#6}{#7}; \exp_stop_f: - \if_int_compare:w #2 = #8#9 \exp_stop_f: - 0 - \else: - \if_int_compare:w #2 < #8#9 - \fi: 1 - \fi: - \or: 1 - \else: -1 - \fi: - } -% \end{macrocode} -% \end{macro} -% \end{macro} -% -% ^^A todo: doc that neither operand may be zero (or fix ep_to_ep above) -% \begin{macro}[EXP]{\@@_ep_mul:wwwwn, \@@_ep_mul_raw:wwwwN} -% Multiply two extended-precision numbers: first normalize them to -% avoid losing too much precision, then multiply the mantissas |#2| -% and~|#4| as fixed point numbers, and sum the exponents |#1| -% and~|#3|. The result's first block is in $[100,9999]$. -% \begin{macrocode} -\cs_new:Npn \@@_ep_mul:wwwwn #1,#2; #3,#4; - { - \@@_ep_to_ep:wwN #3,#4; - \@@_fixed_continue:wn - { - \@@_ep_to_ep:wwN #1,#2; - \@@_ep_mul_raw:wwwwN - } - \@@_fixed_continue:wn - } -\cs_new:Npn \@@_ep_mul_raw:wwwwN #1,#2; #3,#4; #5 - { - \@@_fixed_mul:wwn #2; #4; - { \exp_after:wN #5 \int_value:w \@@_int_eval:w #1 + #3 , } - } -% \end{macrocode} -% \end{macro} -% -% \subsection{Dividing extended-precision numbers} -% -% \newcommand{\eTeXfrac}[2]{\left[\frac{#1}{#2}\right]} -% -% Divisions of extended-precision numbers are difficult to perform with -% exact rounding: the technique used in \pkg{l3fp-basics} for $16$-digit -% floating point numbers does not generalize easily to $24$-digit -% numbers. Thankfully, there is no need for exact rounding. -% -% Let us call \meta{n} the numerator and \meta{d} the denominator. -% After a simple normalization step, we can assume that -% $\meta{n}\in[0.1,1)$ and $\meta{d}\in[0.1,1)$, and compute -% $\meta{n}/(10\meta{d})\in(0.01,1)$. In terms of the $6$~blocks of -% digits $\meta{n_1}\cdots\meta{n_6}$ and the $6$~blocks -% $\meta{d_1}\cdots\meta{d_6}$, the condition translates to -% $\meta{n_1},\meta{d_1}\in[1000,9999]$. -% -% We first find an integer estimate $a \simeq 10^{8} / \meta{d}$ by -% computing -% \begin{align*} -% \alpha &= \eTeXfrac{10^{9}}{\meta{d_1}+1} \\ -% \beta &= \eTeXfrac{10^{9}}{\meta{d_1}} \\ -% a &= 10^{3} \alpha + (\beta-\alpha) \cdot -% \left(10^{3}-\eTeXfrac{\meta{d_2}}{10}\right) - 1250, -% \end{align*} -% where $\eTeXfrac{\bullet}{\bullet}$ denotes \eTeX{}'s rounding -% division, which rounds ties away from zero. The idea is to -% interpolate between $10^{3}\alpha$ and $10^{3}\beta$ with a parameter -% $\meta{d_2}/10^{4}$, so that when $\meta{d_2}=0$ one gets $a = -% 10^{3}\beta-1250 \simeq 10^{12} / \meta{d_1} \simeq 10^{8} / -% \meta{d}$, while when $\meta{d_2}=9999$ one gets $a = -% 10^{3}\alpha-1250 \simeq 10^{12} / (\meta{d_1} + 1) \simeq 10^{8} / -% \meta{d}$. The shift by $1250$ helps to ensure that $a$ is an -% underestimate of the correct value. We shall prove that -% \[ -% 1 - 1.755\cdot 10^{-5} < \frac{\meta{d}a}{10^{8}} < 1 . -% \] -% We can then compute the inverse of $\meta{d}a/10^{8} = 1 - \epsilon$ -% using the relation $1/(1-\epsilon) \simeq (1+\epsilon)(1+\epsilon^{2}) -% + \epsilon^{4}$, which is correct up to a relative error of -% $\epsilon^5 < 1.6\cdot 10^{-24}$. This allows us to find the desired -% ratio as -% \[ -% \frac{\meta{n}}{\meta{d}} -% = \frac{\meta{n}a}{10^{8}} -% \bigl( (1+\epsilon)(1+\epsilon^{2}) + \epsilon^{4}\bigr) . -% \] -% -% Let us prove the upper bound first (multiplied by $10^{15}$). Note -% that $10^{7} \meta{d} < 10^{3} \meta{d_1} + 10^{-1} (\meta{d_2} + 1)$, -% and that \eTeX{}'s division $\eTeXfrac{\meta{d_2}}{10}$ underestimates -% $10^{-1}(\meta{d_2} + 1)$ by $0.5$ at most, as can be checked -% for each possible last digit of \meta{d_2}. Then, -% \begin{align} -% 10^{7} \meta{d}a -% & < -% \left(10^{3}\meta{d_1} -% + \eTeXfrac{\meta{d_2}}{10} + \frac{1}{2}\right) -% \left(\left(10^{3}-\eTeXfrac{\meta{d_2}}{10}\right) \beta -% + \eTeXfrac{\meta{d_2}}{10} \alpha - 1250\right) -% \\ -% & < -% \left(10^{3}\meta{d_1} -% + \eTeXfrac{\meta{d_2}}{10} + \frac{1}{2}\right) -% \\ & \qquad -% \left( -% \left(10^{3}-\eTeXfrac{\meta{d_2}}{10}\right) -% \left(\frac{10^{9}}{\meta{d_1}} + \frac{1}{2} \right) -% + \eTeXfrac{\meta{d_2}}{10} -% \left(\frac{10^{9}}{\meta{d_1}+1} + \frac{1}{2} \right) -% - 1250 -% \right) -% \\ -% & < -% \left(10^{3} \meta{d_1} -% + \eTeXfrac{\meta{d_2}}{10} + \frac{1}{2}\right) -% \left(\frac{10^{12}}{\meta{d_1}} -% - \eTeXfrac{\meta{d_2}}{10} -% \frac{10^{9}}{\meta{d_1}(\meta{d_1}+1)} -% - 750\right) -% \end{align} -% We recognize a quadratic polynomial in $[\meta{d_2}/10]$ with a -% negative leading coefficient: this polynomial is bounded above, -% according to $([\meta{d_2}/10]+a)(b-c[\meta{d_2}/10]) \leq -% (b+ca)^2/(4c)$. Hence, -% \[ -% 10^{7} \meta{d}a -% < \frac{10^{15}}{\meta{d_1}(\meta{d_1}+1)} \left( -% \meta{d_1} + \frac{1}{2} + \frac{1}{4} 10^{-3} -% - \frac{3}{8} \cdot 10^{-9} \meta{d_1}(\meta{d_1}+1) \right)^2 -% \] -% Since \meta{d_1} takes integer values within $[1000,9999]$, it is a -% simple programming exercise to check that the squared expression is -% always less than $\meta{d_1}(\meta{d_1}+1)$, hence $10^{7} \meta{d} a -% < 10^{15}$. The upper bound is proven. We also find that -% $\frac{3}{8}$ can be replaced by slightly smaller numbers, but nothing -% less than $0.374563\ldots$, and going back through the derivation of -% the upper bound, we find that $1250$ is as small a shift as we can -% obtain without breaking the bound. -% -% Now, the lower bound. The same computation as for the upper bound -% implies -% \[ -% 10^{7} \meta{d}a -% > \left(10^{3} \meta{d_1} + \eTeXfrac{\meta{d_2}}{10} -% - \frac{1}{2}\right) -% \left(\frac{10^{12}}{\meta{d_1}} -% - \eTeXfrac{\meta{d_2}}{10} \frac{10^{9}}{\meta{d_1}(\meta{d_1}+1)} -% - 1750\right) -% \] -% This time, we want to find the minimum of this quadratic polynomial. -% Since the leading coefficient is still negative, the minimum is -% reached for one of the extreme values $[y/10]=0$ or $[y/10]=100$, and -% we easily check the bound for those values. -% -% We have proven that the algorithm gives us a precise enough -% answer. Incidentally, the upper bound that we derived tells us that -% $a < 10^{8}/\meta{d} \leq 10^{9}$, hence we can compute $a$ safely as -% a \TeX{} integer, and even add $10^{9}$ to it to ease grabbing of all -% the digits. The lower bound implies $10^{8} - 1755 < a$, which we do -% not care about. -% -% ^^A todo: provide ep_inv, not ep_div? -% ^^A todo: make extra sure that the result's first block cannot be 99 -% ^^A todo: doc that neither operand may be zero (or fix ep_to_ep) -% \begin{macro}[EXP]{\@@_ep_div:wwwwn} -% Compute the ratio of two extended-precision numbers. The result is -% an extended-precision number whose first block lies in the range -% $[100,9999]$, and is placed after the \meta{continuation} once we -% are done. First normalize the inputs so that both first block lie -% in $[1000,9999]$, then call \cs{@@_ep_div_esti:wwwwn} -% \meta{denominator} \meta{numerator}, responsible for estimating the -% inverse of the denominator. -% \begin{macrocode} -\cs_new:Npn \@@_ep_div:wwwwn #1,#2; #3,#4; - { - \@@_ep_to_ep:wwN #1,#2; - \@@_fixed_continue:wn - { - \@@_ep_to_ep:wwN #3,#4; - \@@_ep_div_esti:wwwwn - } - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP] -% { -% \@@_ep_div_esti:wwwwn, -% \@@_ep_div_estii:wwnnwwn, -% \@@_ep_div_estiii:NNNNNwwwn -% } -% The \texttt{esti} function evaluates $\alpha=10^{9} / (\meta{d_1} + -% 1)$, which is used twice in the expression for $a$, and combines the -% exponents |#1| and~|#4| (with a shift by~$1$ because we later compute -% $\meta{n}/(10\meta{d})$. Then the \texttt{estii} function evaluates -% $10^{9} + a$, and puts the exponent~|#2| after the -% continuation~|#7|: from there on we can forget exponents and focus -% on the mantissa. The \texttt{estiii} function multiplies the -% denominator~|#7| by $10^{-8}a$ (obtained as $a$ split into the -% single digit~|#1| and two blocks of $4$~digits, |#2#3#4#5| -% and~|#6|). The result $10^{-8}a\meta{d}=(1-\epsilon)$, and a -% partially packed $10^{-9}a$ (as a block of four digits, and five -% individual digits, not packed by lack of available macro parameters -% here) are passed to \cs{@@_ep_div_epsi:wnNNNNn}, which computes -% $10^{-9}a/(1-\epsilon)$, that is, $1/(10\meta{d})$ and we finally -% multiply this by the numerator~|#8|. -% \begin{macrocode} -\cs_new:Npn \@@_ep_div_esti:wwwwn #1,#2#3; #4, - { - \exp_after:wN \@@_ep_div_estii:wwnnwwn - \int_value:w \@@_int_eval:w 10 0000 0000 / ( #2 + 1 ) - \exp_after:wN ; - \int_value:w \@@_int_eval:w #4 - #1 + 1 , - {#2} #3; - } -\cs_new:Npn \@@_ep_div_estii:wwnnwwn #1; #2,#3#4#5; #6; #7 - { - \exp_after:wN \@@_ep_div_estiii:NNNNNwwwn - \int_value:w \@@_int_eval:w 10 0000 0000 - 1750 - + #1 000 + (10 0000 0000 / #3 - #1) * (1000 - #4 / 10) ; - {#3}{#4}#5; #6; { #7 #2, } - } -\cs_new:Npn \@@_ep_div_estiii:NNNNNwwwn 1#1#2#3#4#5#6; #7; - { - \@@_fixed_mul_short:wwn #7; {#1}{#2#3#4#5}{#6}; - \@@_ep_div_epsi:wnNNNNNn {#1#2#3#4}#5#6 - \@@_fixed_mul:wwn - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP] -% { -% \@@_ep_div_epsi:wnNNNNNn, -% \@@_ep_div_eps_pack:NNNNNw, -% \@@_ep_div_epsii:wwnNNNNNn, -% } -% The bounds shown above imply that the \texttt{epsi} function's first -% operand is $(1-\epsilon)$ with $\epsilon\in[0,1.755\cdot 10^{-5}]$. -% The \texttt{epsi} function computes $\epsilon$ as $1-(1-\epsilon)$. -% Since $\epsilon<10^{-4}$, its first block vanishes and there is no -% need to explicitly use~|#1| (which is $9999$). Then \texttt{epsii} -% evaluates $10^{-9}a/(1-\epsilon)$ as -% $(1+\epsilon^2)(1+\epsilon)(10^{-9}a \epsilon) + 10^{-9}a$. -% Importantly, we compute $10^{-9}a \epsilon$ before multiplying it -% with the rest, rather than multiplying by $\epsilon$ and then -% $10^{-9}a$, as this second option loses more precision. Also, the -% combination of \texttt{short_mul} and \texttt{div_myriad} is both -% faster and more precise than a simple \texttt{mul}. -% \begin{macrocode} -\cs_new:Npn \@@_ep_div_epsi:wnNNNNNn #1#2#3#4#5#6; - { - \exp_after:wN \@@_ep_div_epsii:wwnNNNNNn - \int_value:w \@@_int_eval:w 1 9998 - #2 - \exp_after:wN \@@_ep_div_eps_pack:NNNNNw - \int_value:w \@@_int_eval:w 1 9999 9998 - #3#4 - \exp_after:wN \@@_ep_div_eps_pack:NNNNNw - \int_value:w \@@_int_eval:w 2 0000 0000 - #5#6 ; ; - } -\cs_new:Npn \@@_ep_div_eps_pack:NNNNNw #1#2#3#4#5#6; - { + #1 ; {#2#3#4#5} {#6} } -\cs_new:Npn \@@_ep_div_epsii:wwnNNNNNn 1#1; #2; #3#4#5#6#7#8 - { - \@@_fixed_mul:wwn {0000}{#1}#2; {0000}{#1}#2; - \@@_fixed_add_one:wN - \@@_fixed_mul:wwn {10000} {#1} #2 ; - { - \@@_fixed_mul_short:wwn {0000}{#1}#2; {#3}{#4#5#6#7}{#8000}; - \@@_fixed_div_myriad:wn - \@@_fixed_mul:wwn - } - \@@_fixed_add:wwn {#3}{#4#5#6#7}{#8000}{0000}{0000}{0000}; - } -% \end{macrocode} -% \end{macro} -% -% \subsection{Inverse square root of extended precision numbers} -% -% The idea here is similar to division. Normalize the input, -% multiplying by powers of $100$ until we have $x\in[0.01,1)$. Then -% find an integer approximation $r \in [101, 1003]$ of -% $10^{2}/\sqrt{x}$, as the fixed point of iterations of the Newton -% method: essentially $r \mapsto (r + 10^{8} / (x_{1} r)) / 2$, starting -% from a guess that optimizes the number of steps before convergence. -% In fact, just as there is a slight shift when computing divisions to -% ensure that some inequalities hold, we replace $10^{8}$ by a -% slightly larger number which ensures that $r^2 x \geq 10^{4}$. -% This also causes $r \in [101, 1003]$. Another correction to the above -% is that the input is actually normalized to $[0.1,1)$, and we use -% either $10^{8}$ or $10^{9}$ in the Newton method, depending on the -% parity of the exponent. Skipping those technical hurdles, once we -% have the approximation~$r$, we set $y = 10^{-4} r^{2} x$ (or rather, -% the correct power of~$10$ to get $y\simeq 1$) and compute $y^{-1/2}$ -% through another application of Newton's method. This time, the -% starting value is $z=1$, each step maps $z \mapsto z(1.5-0.5yz^2)$, -% and we perform a fixed number of steps. Our final result combines~$r$ -% with $y^{-1/2}$ as $x^{-1/2} = 10^{-2} r y^{-1/2}$. -% -% ^^A todo: doc that the operand may not be zero (or fix ep_to_ep above) -% \begin{macro}[EXP]{\@@_ep_isqrt:wwn} -% \begin{macro}[EXP] -% {\@@_ep_isqrt_aux:wwn, \@@_ep_isqrt_auxii:wwnnnwn} -% First normalize the input, then check the parity of the -% exponent~|#1|. If it is even, the result's exponent will be -% $-|#1|/2$, otherwise it will be $(|#1|-1)/2$ (except in the case -% where the input was an exact power of $100$). The \texttt{auxii} -% function receives as~|#1| the result's exponent just computed, as -% |#2| the starting value for the iteration giving~$r$ (the -% values~$168$ and~$535$ lead to the least number of iterations before -% convergence, on average), as |#3| and~|#4| one empty argument and -% one~|0|, depending on the parity of the original exponent, as |#5| -% and~|#6| the normalized mantissa ($|#5|\in[1000,9999]$), and as |#7| -% the continuation. It sets up the iteration giving~$r$: the -% \texttt{esti} function thus receives the initial two guesses |#2| -% and~$0$, an approximation~|#5| of~$10^{4}x$ (its first block of -% digits), and the empty/zero arguments |#3| and~|#4|, followed by the -% mantissa and an altered continuation where we have stored the -% result's exponent. -% \begin{macrocode} -\cs_new:Npn \@@_ep_isqrt:wwn #1,#2; - { - \@@_ep_to_ep:wwN #1,#2; - \@@_ep_isqrt_auxi:wwn - } -\cs_new:Npn \@@_ep_isqrt_auxi:wwn #1, - { - \exp_after:wN \@@_ep_isqrt_auxii:wwnnnwn - \int_value:w \@@_int_eval:w - \int_if_odd:nTF {#1} - { (1 - #1) / 2 , 535 , { 0 } { } } - { 1 - #1 / 2 , 168 , { } { 0 } } - } -\cs_new:Npn \@@_ep_isqrt_auxii:wwnnnwn #1, #2, #3#4 #5#6; #7 - { - \@@_ep_isqrt_esti:wwwnnwn #2, 0, #5, {#3} {#4} - {#5} #6 ; { #7 #1 , } - } -% \end{macrocode} -% \end{macro} -% \end{macro} -% -% \begin{macro}[EXP] -% { -% \@@_ep_isqrt_esti:wwwnnwn, -% \@@_ep_isqrt_estii:wwwnnwn, -% \@@_ep_isqrt_estiii:NNNNNwwwn -% } -% If the last two approximations gave the same result, we are done: -% call the \texttt{estii} function to clean up. Otherwise, evaluate -% $(\meta{prev} + 1.005 \cdot 10^{\text{$8$ or $9$}} / (\meta{prev} -% \cdot x)) / 2$, as the next approximation: omitting the $1.005$ -% factor, this would be Newton's method. We can check by brute force -% that if |#4| is empty (the original exponent was even), the process -% computes an integer slightly larger than $100 / \sqrt{x}$, while if -% |#4| is~$0$ (the original exponent was odd), the result is an -% integer slightly larger than $100 / \sqrt{x/10}$. Once we are done, -% we evaluate $100 r^2 / 2$ or $10 r^2 / 2$ (when the exponent is even -% or odd, respectively) and feed that to \texttt{estiii}. This third -% auxiliary finds $y_{\text{even}} / 2 = 10^{-4} r^2 x / 2$ or -% $y_{\text{odd}} / 2 = 10^{-5} r^2 x / 2$ (again, depending on -% earlier parity). A simple program shows that $y\in [1, 1.0201]$. -% The number $y/2$ is fed to \cs{@@_ep_isqrt_epsi:wN}, which computes -% $1/\sqrt{y}$, and we finally multiply the result by~$r$. -% \begin{macrocode} -\cs_new:Npn \@@_ep_isqrt_esti:wwwnnwn #1, #2, #3, #4 - { - \if_int_compare:w #1 = #2 \exp_stop_f: - \exp_after:wN \@@_ep_isqrt_estii:wwwnnwn - \fi: - \exp_after:wN \@@_ep_isqrt_esti:wwwnnwn - \int_value:w \@@_int_eval:w - (#1 + 1 0050 0000 #4 / (#1 * #3)) / 2 , - #1, #3, {#4} - } -\cs_new:Npn \@@_ep_isqrt_estii:wwwnnwn #1, #2, #3, #4#5 - { - \exp_after:wN \@@_ep_isqrt_estiii:NNNNNwwwn - \int_value:w \@@_int_eval:w 1000 0000 + #2 * #2 #5 * 5 - \exp_after:wN , \int_value:w \@@_int_eval:w 10000 + #2 ; - } -\cs_new:Npn \@@_ep_isqrt_estiii:NNNNNwwwn 1#1#2#3#4#5#6, 1#7#8; #9; - { - \@@_fixed_mul_short:wwn #9; {#1} {#2#3#4#5} {#600} ; - \@@_ep_isqrt_epsi:wN - \@@_fixed_mul_short:wwn {#7} {#80} {0000} ; - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[EXP]{\@@_ep_isqrt_epsi:wN, \@@_ep_isqrt_epsii:wwN} -% Here, we receive a fixed point number $y/2$ with $y\in[1,1.0201]$. -% Starting from $z = 1$ we iterate $z \mapsto z(3/2 - z^2 y/2)$. In -% fact, we start from the first iteration $z=3/2-y/2$ to avoid useless -% multiplications. The \texttt{epsii} auxiliary receives $z$ as~|#1| -% and $y$ as~|#2|. -% \begin{macrocode} -\cs_new:Npn \@@_ep_isqrt_epsi:wN #1; - { - \@@_fixed_sub:wwn {15000}{0000}{0000}{0000}{0000}{0000}; #1; - \@@_ep_isqrt_epsii:wwN #1; - \@@_ep_isqrt_epsii:wwN #1; - \@@_ep_isqrt_epsii:wwN #1; - } -\cs_new:Npn \@@_ep_isqrt_epsii:wwN #1; #2; - { - \@@_fixed_mul:wwn #1; #1; - \@@_fixed_mul_sub_back:wwwn #2; - {15000}{0000}{0000}{0000}{0000}{0000}; - \@@_fixed_mul:wwn #1; - } -% \end{macrocode} -% \end{macro} -% -% \subsection{Converting from fixed point to floating point} -% ^^A todo: doc -% -% After computing Taylor series, we wish to convert the result from -% extended precision (with or without an exponent) to the public -% floating point format. The functions here should be called within an -% integer expression for the overall exponent of the floating point. -% -% \begin{macro}[rEXP]{\@@_ep_to_float_o:wwN, \@@_ep_inv_to_float_o:wwN} -% An extended-precision number is simply a comma-delimited exponent -% followed by a fixed point number. Leave the exponent in the current -% integer expression then convert the fixed point number. -% \begin{macrocode} -\cs_new:Npn \@@_ep_to_float_o:wwN #1, - { + \@@_int_eval:w #1 \@@_fixed_to_float_o:wN } -\cs_new:Npn \@@_ep_inv_to_float_o:wwN #1,#2; - { - \@@_ep_div:wwwwn 1,{1000}{0000}{0000}{0000}{0000}{0000}; #1,#2; - \@@_ep_to_float_o:wwN - } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[rEXP]{\@@_fixed_inv_to_float_o:wN} -% Another function which reduces to converting an extended precision -% number to a float. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_inv_to_float_o:wN - { \@@_ep_inv_to_float_o:wwN 0, } -% \end{macrocode} -% \end{macro} -% -% \begin{macro}[rEXP]{\@@_fixed_to_float_rad_o:wN} -% Converts the fixed point number~|#1| from degrees to radians then to -% a floating point number. This could perhaps remain in -% \pkg{l3fp-trig}. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_to_float_rad_o:wN #1; - { - \@@_fixed_mul:wwn #1; {5729}{5779}{5130}{8232}{0876}{7981}; - { \@@_ep_to_float_o:wwN 2, } - } -% \end{macrocode} -% \end{macro} -% -% ^^A todo: make exponents end in ',' consistently throughout l3fp -% \begin{macro}[rEXP] -% {\@@_fixed_to_float_o:wN, \@@_fixed_to_float_o:Nw} -% \begin{syntax} -% \ldots{} \cs{@@_int_eval:w} \meta{exponent} \cs{@@_fixed_to_float_o:wN} \Arg{a_1} \Arg{a_2} \Arg{a_3} \Arg{a_4} \Arg{a_5} \Arg{a_6} |;| \meta{sign} -% \end{syntax} -% yields -% \begin{quote} -% \meta{exponent'} |;| \Arg{a'_1} \Arg{a'_2} \Arg{a'_3} \Arg{a'_4} |;| -% \end{quote} -% And the \texttt{to_fixed} version gives six brace groups instead of -% $4$, ensuring that $1000\leq\meta{a'_1}\leq 9999$. At this stage, we -% know that \meta{a_1} is positive (otherwise, it is sign of an error -% before), and we assume that it is less than $10^8$.\footnote{Bruno: -% I must double check this assumption.} -% -%^^A todo: round properly when rounding to infinity: I need the sign. -% \begin{macrocode} -\cs_new:Npn \@@_fixed_to_float_o:Nw #1#2; - { \@@_fixed_to_float_o:wN #2; #1 } -\cs_new:Npn \@@_fixed_to_float_o:wN #1#2#3#4#5#6; #7 - { % for the 8-digit-at-the-start thing - + \@@_int_eval:w \c_@@_block_int - \exp_after:wN \exp_after:wN - \exp_after:wN \@@_fixed_to_loop:N - \exp_after:wN \use_none:n - \int_value:w \@@_int_eval:w - 1 0000 0000 + #1 \exp_after:wN \@@_use_none_stop_f:n - \int_value:w 1#2 \exp_after:wN \@@_use_none_stop_f:n - \int_value:w 1#3#4 \exp_after:wN \@@_use_none_stop_f:n - \int_value:w 1#5#6 - \exp_after:wN ; - \exp_after:wN ; - } -\cs_new:Npn \@@_fixed_to_loop:N #1 - { - \if_meaning:w 0 #1 - - 1 - \exp_after:wN \@@_fixed_to_loop:N - \else: - \exp_after:wN \@@_fixed_to_loop_end:w - \exp_after:wN #1 - \fi: - } -\cs_new:Npn \@@_fixed_to_loop_end:w #1 #2 ; - { - \if_meaning:w ; #1 - \exp_after:wN \@@_fixed_to_float_zero:w - \else: - \exp_after:wN \@@_pack_twice_four:wNNNNNNNN - \exp_after:wN \@@_pack_twice_four:wNNNNNNNN - \exp_after:wN \@@_fixed_to_float_pack:ww - \exp_after:wN ; - \fi: - #1 #2 0000 0000 0000 0000 ; - } -\cs_new:Npn \@@_fixed_to_float_zero:w ; 0000 0000 0000 0000 ; - { - - 2 * \c_@@_max_exponent_int ; - {0000} {0000} {0000} {0000} ; - } -\cs_new:Npn \@@_fixed_to_float_pack:ww #1 ; #2#3 ; ; - { - \if_int_compare:w #2 > 4 \exp_stop_f: - \exp_after:wN \@@_fixed_to_float_round_up:wnnnnw - \fi: - ; #1 ; - } -\cs_new:Npn \@@_fixed_to_float_round_up:wnnnnw ; #1#2#3#4 ; - { - \exp_after:wN \@@_basics_pack_high:NNNNNw - \int_value:w \@@_int_eval:w 1 #1#2 - \exp_after:wN \@@_basics_pack_low:NNNNNw - \int_value:w \@@_int_eval:w 1 #3#4 + 1 ; - } -% \end{macrocode} -% \end{macro} -% -% \begin{macrocode} -%</package> -% \end{macrocode} -% -% \end{implementation} -% -% \PrintChanges -% -% \PrintIndex |