diff options
Diffstat (limited to 'Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx')
| -rw-r--r-- | Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx | 499 |
1 files changed, 491 insertions, 8 deletions
diff --git a/Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx b/Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx index 104a118ca3b..d5ddae31ab0 100644 --- a/Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx +++ b/Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx @@ -36,7 +36,7 @@ % %<*driver> \RequirePackage{l3bootstrap} -\GetIdInfo$Id: l3fp-basics.dtx 4601 2013-11-18 23:13:28Z bruno $ +\GetIdInfo$Id: l3fp-basics.dtx 4619 2013-12-14 12:53:54Z bruno $ {L3 Floating-point arithmetic} \documentclass[full]{l3doc} \begin{document} @@ -107,12 +107,7 @@ % |{0000}|. % \begin{macrocode} \cs_new:Npn \@@_basics_pack_low:NNNNNw #1 #2#3#4#5 #6; - { - \if_meaning:w 2 #1 - + \c_one - \fi: - ; {#2#3#4#5} {#6} ; - } + { + #1 - \c_one ; {#2#3#4#5} {#6} ; } \cs_new:Npn \@@_basics_pack_high:NNNNNw #1 #2#3#4#5 #6; { \if_meaning:w 2 #1 @@ -1599,7 +1594,495 @@ % \end{macrocode} % \end{macro} % -% \subsection{Unary operations} +% \subsection{Square root} +% +% \begin{macro}[int, EXP]{\@@_sqrt_o:w} +% Zeros are unchanged: $\sqrt{-0} = -0$ and $\sqrt{+0} = +0$. +% Negative numbers (other than $-0$) have no real square root. +% Positive infinity, and \texttt{nan}, are unchanged. Finally, for +% normal positive numbers, there is some work to do. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_o:w #1 \s_@@ \@@_chk:w #2#3#4; @ + { + \if_meaning:w 0 #2 \@@_case_return_same_o:w \fi: + \if_meaning:w 2 #3 + \@@_case_use:nw { \@@_invalid_operation_o:nw { sqrt } } + \fi: + \if_meaning:w 1 #2 \else: \@@_case_return_same_o:w \fi: + \@@_sqrt_npos_o:w + \s_@@ \@@_chk:w #2 #3 #4; + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, EXP]{\@@_sqrt_npos_o:w} +% \begin{macro}[aux, rEXP] +% {\@@_sqrt_npos_auxi_o:wwnnN, \@@_sqrt_npos_auxii_o:wNNNNNNNN} +% Prepare \cs{@@_sanitize:Nw} to receive the final sign~|0| (the +% result is always positive) and the exponent, equal to half of the +% exponent~|#1| of the argument. If the exponent~|#1| is even, find a +% first approximation of the square root of the significand $10^{8} +% a_1 + a_2 = 10^{8} |#2#3| + |#4#5|$ through Newton's method, +% starting at $x = 57234133 \simeq 10^{7.75}$. Otherwise, first shift +% the significand of of the argument by one digit, getting +% $a_1'\in[10^{6}, 10^{7})$ instead of $[10^{7}, 10^{8})$, then use +% Newton's method starting at $17782794 \simeq 10^{7.25}$. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_npos_o:w \s_@@ \@@_chk:w 1 0 #1#2#3#4#5; + { + \exp_after:wN \@@_sanitize:Nw + \exp_after:wN 0 + \int_use:N \__int_eval:w + \if_int_odd:w #1 \exp_stop_f: + \exp_after:wN \@@_sqrt_npos_auxi_o:wwnnN + \fi: + #1 / \c_two + \@@_sqrt_Newton_o:wwn 56234133; 0; {#2#3} {#4#5} 0 + } +\cs_new:Npn \@@_sqrt_npos_auxi_o:wwnnN #1 / \c_two #2; 0; #3#4#5 + { + ( #1 + \c_one ) / \c_two + \@@_pack_eight:wNNNNNNNN + \@@_sqrt_npos_auxii_o:wNNNNNNNN + ; + 0 #3 #4 + } +\cs_new:Npn \@@_sqrt_npos_auxii_o:wNNNNNNNN #1; #2#3#4#5#6#7#8#9 + { \@@_sqrt_Newton_o:wwn 17782794; 0; {#1} {#2#3#4#5#6#7#8#9} } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP]{\@@_sqrt_Newton_o:wwn} +% Newton's method maps $x\mapsto\bigl[(x + [10^{8} a_1 / x])/2\bigr]$ +% in each iteration, where $[b/c]$ denotes \eTeX{}'s division. This +% division rounds the real number $b/c$ to the closest integer, +% rounding ties away from zero, hence when $c$~is even, +% $b/c - 1/2 + 1/c \leq [b/c] \leq b/c + 1/2$ +% and when $c$~is odd, +% $b/c - 1/2 + 1/(2c) \leq [b/c] \leq b/c + 1/2 - 1/(2c)$. +% For all~$c$, $b/c - 1/2 + 1/(2c) \leq [b/c] \leq b/c + 1/2$. +% +% Let us prove that the method converges when implemented with \eTeX{} +% integer division, for any $10^{6} \leq a_1 < 10^{8}$ and starting +% value $10^{6} \leq x < 10^{8}$. Using the inequalities above and +% the arithmetic--geometric inequality $(x+t)/2 \geq \sqrt{xt}$ for $t +% = 10^{8} a_1 / x$, we find +% \[ +% x' +% = \left[\frac{x + [10^{8} a_1 / x]}{2}\right] +% \geq \frac{x + 10^{8} a_1 / x - 1/2 + 1/(2x)}{2} +% \geq \sqrt{10^{8} a_1} - \frac{1}{4} + \frac{1}{4x} \,. +% \] +% After any step of iteration, we thus have $\delta = x - \sqrt{10^{8} +% a_1} \geq -0.25 + 0.25 \cdot 10^{-8}$. The new difference +% $\delta' = x' - \sqrt{10^{8} a_1}$ after one step is bounded above +% as +% \[ +% x' - \sqrt{10^{8} a_1} +% \leq \frac{x + 10^{8} a_1 / x + 1/2}{2} + \frac{1}{2} +% - \sqrt{10^{8} a_1} +% \leq \frac{\delta}{2} \frac{\delta}{\sqrt{10^{8} a_1} + \delta} +% + \frac{3}{4} \,. +% \] +% For $\delta > 3/2$, this last expression is +% $\leq\delta/2+3/4<\delta$, hence $\delta$~decreases at each step: +% since all~$x$ are integers, $\delta$~must reach a value +% $-1/4<\delta\leq 3/2$. In this range of values, we get $\delta' +% \leq \frac{3}{4} \frac{3}{2\sqrt{10^{8} a_1}} + \frac{3}{4} \leq +% 0.75 + 1.125 \cdot 10^{-7}$. We deduce that the difference $\delta +% = x - \sqrt{10^{8} a_1}$ eventually reaches a value in the interval +% $[-0.25 + 0.25\cdot 10^{-8}, 0.75 + 11.25 \cdot 10^{-8}]$, whose +% width is $1 + 11 \cdot 10^{-8}$. The corresponding interval for~$x$ +% may contain two integers, hence $x$~might oscillate between those +% two values. +% +% However, the fact that $x\mapsto x-1$ and $x-1 \mapsto x$ puts +% stronger constraints, which are not compatible: the first implies +% \[ +% x + [10^{8} a_1 / x] \leq 2x - 2 +% \] +% hence $10^{8} a_1 / x \leq x - 3/2$, while the second implies +% \[ +% x - 1 + [10^{8} a_1 / (x - 1)] \geq 2x - 1 +% \] +% hence $10^{8} a_1 / (x - 1) \geq x - 1/2$. Combining the two +% inequalities yields $x^2 - 3x/2 \geq 10^{8} a_1 \geq x - 3x/2 + +% 1/2$, which cannot hold. Therefore, the iteration always converges +% to a single integer~$x$. To stop the iteration when two consecutive +% results are equal, the function \cs{@@_sqrt_Newton_o:wwn} receives +% the newly computed result as~|#1|, the previous result as~|#2|, and +% $a_1$ as~|#3|. Note that \eTeX{} combines the computation of a +% multiplication and a following division, thus avoiding overflow in +% |#3 * 100000000 / #1|. In any case, the result is within $[10^{7}, +% 10^{8}]$. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_Newton_o:wwn #1; #2; #3 + { + \if_int_compare:w #1 = #2 \exp_stop_f: + \exp_after:wN \@@_sqrt_auxi_o:NNNNwnnN + \int_use:N \__int_eval:w 9999 9999 + + \exp_after:wN \@@_use_none_until_s:w + \fi: + \exp_after:wN \@@_sqrt_Newton_o:wwn + \int_use:N \__int_eval:w (#1 + #3 * 1 0000 0000 / #1) / \c_two ; + #1; {#3} + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP]{\@@_sqrt_auxi_o:NNNNwnnN} +% This function is followed by $10^{8}+x-1$, which has~$9$ digits +% starting with~$1$, then |;| \Arg{a_1} \Arg{a_2} \meta{a'}. Here, $x +% \simeq \sqrt{10^{8} a_1}$ and we want to estimate the square root of +% $a = 10^{-8} a_1 + 10^{-16} a_2 + 10^{-17} a'$. We set up an +% initial underestimate +% \[ +% y = (x - 1) 10^{-8} + 0.2499998875 \cdot 10^{-8} \lesssim \sqrt{a}\,. +% \] +% From the inequalities shown earlier, we know that $y \leq +% \sqrt{10^{-8} a_1} \leq \sqrt{a}$ and that $\sqrt{10^{-8} a_1} \leq +% y + 10^{-8} + 11\cdot 10^{-16}$ hence (using $0.1\leq y\leq +% \sqrt{a}\leq 1$) +% \[ +% a - y^2 \leq 10^{-8} a_1 + 10^{-8} - y^2 +% \leq (y + 10^{-8} + 11\cdot 10^{-16})^2 - y^2 + 10^{-8} +% < 3.2 \cdot 10^{-8} \,, +% \] +% and $\sqrt{a} - y = (a - y^2)/(\sqrt{a} + y) \leq 16 \cdot 10^{-8}$. +% Next, \cs{@@_sqrt_auxii_o:NnnnnnnnN} will be called several times to +% get closer and closer underestimates of~$\sqrt{a}$. By +% construction, the underestimates~$y$ are always increasing, $a - y^2 +% < 3.2 \cdot 10^{-8}$ for all. Also, $y<1$. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_auxi_o:NNNNwnnN 1 #1#2#3#4#5; + { + \@@_sqrt_auxii_o:NnnnnnnnN + \@@_sqrt_auxiii_o:wnnnnnnnn + {#1#2#3#4} {#5} {2499} {9988} {7500} + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP]{\@@_sqrt_auxii_o:NnnnnnnnN} +% This receives a continuation function~|#1|, then five blocks of~$4$ +% digits for~$y$, then two $8$-digit blocks and a single digit +% for~$a$. A common estimate of $\sqrt{a} - y = (a - y^2) / (\sqrt{a} +% + y)$ is $(a - y^2)/(2y)$, which leads to alternating overestimates +% and underestimates. We tweak this, to only work with underestimates +% (no need then to worry about signs in the computation). Each step +% finds the largest integer $j\leq 6$ such that $10^{4j}(a-y^2) < +% 2\cdot 10^{8}$, then computes the integer (with \eTeX{}'s rounding +% division) +% \[ +% 10^{4j} z = +% \Bigl[\bigl(\lfloor 10^{4j}(a-y^2)\rfloor - 257\bigr) +% \cdot (0.5\cdot 10^{8}) +% \Bigm/ \lfloor 10^{8} y + 1\rfloor\Bigr] \,. +% \] +% The choice of~$j$ ensures that $10^{4j} z < 2\cdot 10^{8} \cdot +% 0.5\cdot 10^{8} / 10^{7} = 10^{9}$, thus $10^{9} + 10^{4j} z$ has +% exactly $10$~digits, does not overflow \TeX{}'s integer range, and +% starts with~$1$. Incidentally, since all $a - y^2 \leq 3.2\cdot +% 10^{-8}$, we know that $j\geq 3$. +% +% Let us show that $z$ is an underestimate of $\sqrt{a} - y$. On the +% one hand, $\sqrt{a} - y \leq 16\cdot 10^{-8}$ because this holds for +% the initial~$y$ and values of~$y$ can only increase. On the other +% hand, the choice of~$j$ implies that $\sqrt{a} - y \leq +% 5(\sqrt{a}+y)(\sqrt{a}-y) = 5(a - y^2) < 10^{9-4j}$. For $j=3$, the +% first bound is better, while for larger~$j$, the second bound is +% better. For all $j\in[3,6]$, we find $\sqrt{a}-y < 16\cdot +% 10^{-2j}$. From this, we deduce that +% \[ +% 10^{4j} (\sqrt{a}-y) +% = \frac{10^{4j}\bigl(a-y^2-(\sqrt{a}-y)^2\bigr)}{2y} +% \geq \frac{\bigl\lfloor 10^{4j}(a-y^2)\bigr\rfloor-257} +% {2\cdot 10^{-8} \lfloor 10^{8}y+1\rfloor} +% + \frac{1}{2} +% \] +% where we have replaced the bound $10^{4j}(16\cdot 10^{-2j}) = 256$ +% by~$257$ and extracted the corresponding term $1/\bigl(2\cdot +% 10^{-8} \lfloor 10^{8}y+1\rfloor\bigr) \geq 1/2$. Given that +% \eTeX{}'s integer division obeys $[b/c] \leq b/c + 1/2$, we deduce +% that $10^{4j} z \leq 10^{4j} (\sqrt{a}-y)$, hence $y+z\leq\sqrt{a}$ +% is an underestimate of~$\sqrt{a}$, as claimed. One implementation +% detail: because the computation involves |-#4*#4| |-| |2*#3*#5| |-| +% |2*#2*#6| which may be as low as $-5\cdot 10^{8}$, we need to use +% the \texttt{pack_big} functions, and the \texttt{big} shifts. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_auxii_o:NnnnnnnnN #1 #2#3#4#5#6 #7#8#9 + { + \exp_after:wN #1 + \int_use:N \__int_eval:w \c_@@_big_leading_shift_int + + #7 - #2 * #2 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + - 2 * #2 * #3 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + + #8 - #3 * #3 - 2 * #2 * #4 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + - 2 * #3 * #4 - 2 * #2 * #5 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + + #9 000 0000 - #4 * #4 - 2 * #3 * #5 - 2 * #2 * #6 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + - 2 * #4 * #5 - 2 * #3 * #6 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + - #5 * #5 - 2 * #4 * #6 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w + \c_@@_big_middle_shift_int + - 2 * #5 * #6 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w + \c_@@_big_trailing_shift_int + - #6 * #6 ; + % ( + - 257 ) * 5000 0000 / (#2#3 + 1) + 10 0000 0000 ; + {#2}{#3}{#4}{#5}{#6} {#7}{#8}#9 + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP] +% { +% \@@_sqrt_auxiii_o:wnnnnnnnn, +% \@@_sqrt_auxiv_o:NNNNNw, +% \@@_sqrt_auxv_o:NNNNNw, +% \@@_sqrt_auxvi_o:NNNNNw, +% \@@_sqrt_auxvii_o:NNNNNw +% } +% We receive here the difference $a-y^2=d=\sum_i d_i \cdot 10^{-4i}$, +% as \meta{d_2} |;| \Arg{d_3} \ldots{} \Arg{d_{10}}, where each block +% has~$4$ digits, except \meta{d_2}. This function finds the largest +% $j\leq 6$ such that $10^{4j}(a-y^2) < 2\cdot 10^{8}$, then leaves an +% open parenthesis and the integer +% $\bigl\lfloor 10^{4j}(a-y^2)\bigr\rfloor$ in an integer +% expression. The closing parenthesis is provided by the caller +% \cs{@@_sqrt_auxii_o:NnnnnnnnN}, which completes the expression +% \[ +% 10^{4j} z = +% \Bigl[\bigl(\lfloor 10^{4j}(a-y^2)\rfloor - 257\bigr) +% \cdot (0.5\cdot 10^{8}) +% \Bigm/ \lfloor 10^{8} y + 1\rfloor\Bigr] +% \] +% for an estimate of $10^{4j} (\sqrt{a} - y)$. If $d_2\geq 2$, $j=3$ +% and the \texttt{auxiv} auxiliary receives $10^{12} z$. If $d_2\leq +% 1$ but $10^{4} d_2 + d_3 \geq 2$, $j=4$ and the \texttt{auxv} +% auxiliary is called, and receives $10^{16} z$, and so on. In all +% those cases, the \texttt{auxviii} auxiliary is set up to add~$z$ +% to~$y$, then go back to the \texttt{auxii} step with continuation +% \texttt{auxiii} (the function we are currently describing). The +% maximum value of $j$ is~$6$, regardless of whether $10^{12} d_2 + +% 10^{8} d_3 + 10^{4} d_4 + d_5 \geq 1$. In this last case, we detect +% when $10^{24} z < 10^{7}$, which essentially means $\sqrt{a} - y +% \lesssim 10^{-17}$: once this threshold is reached, there is enough +% information to find the correctly rounded~$\sqrt{a}$ with only one +% more call to \cs{@@_sqrt_auxii_o:NnnnnnnnN}. Note that the +% iteration cannot be stuck before reaching $j=6$, because for $j<6$, +% one has $2\cdot 10^{8}\leq 10^{4(j+1)}(a-y^2)$, hence +% \[ +% 10^{4j} z +% \geq \frac{(20000-257)(0.5\cdot 10^{8})}{\lfloor 10^{8} y + 1\rfloor} +% \geq (20000-257)\cdot 0.5 > 0 \,. +% \] +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_auxiii_o:wnnnnnnnn + #1; #2#3#4#5#6#7#8#9 + { + \if_int_compare:w #1 > \c_one + \exp_after:wN \@@_sqrt_auxiv_o:NNNNNw + \int_use:N \__int_eval:w (#1#2 %) + \else: + \if_int_compare:w #1#2 > \c_one + \exp_after:wN \@@_sqrt_auxv_o:NNNNNw + \int_use:N \__int_eval:w (#1#2#3 %) + \else: + \if_int_compare:w #1#2#3 > \c_one + \exp_after:wN \@@_sqrt_auxvi_o:NNNNNw + \int_use:N \__int_eval:w (#1#2#3#4 %) + \else: + \exp_after:wN \@@_sqrt_auxvii_o:NNNNNw + \int_use:N \__int_eval:w (#1#2#3#4#5 %) + \fi: + \fi: + \fi: + } +\cs_new:Npn \@@_sqrt_auxiv_o:NNNNNw 1#1#2#3#4#5#6; + { \@@_sqrt_auxviii_o:nnnnnnn {#1#2#3#4#5#6} {00000000} } +\cs_new:Npn \@@_sqrt_auxv_o:NNNNNw 1#1#2#3#4#5#6; + { \@@_sqrt_auxviii_o:nnnnnnn {000#1#2#3#4#5} {#60000} } +\cs_new:Npn \@@_sqrt_auxvi_o:NNNNNw 1#1#2#3#4#5#6; + { \@@_sqrt_auxviii_o:nnnnnnn {0000000#1} {#2#3#4#5#6} } +\cs_new:Npn \@@_sqrt_auxvii_o:NNNNNw 1#1#2#3#4#5#6; + { + \if_int_compare:w #1#2 = \c_zero + \exp_after:wN \@@_sqrt_auxx_o:Nnnnnnnn + \fi: + \@@_sqrt_auxviii_o:nnnnnnn {00000000} {000#1#2#3#4#5} + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP] +% {\@@_sqrt_auxviii_o:nnnnnnn, \@@_sqrt_auxix_o:wnwnw} +% Simply add the two $8$-digit blocks of~$z$, aligned to the last four +% of the five $4$-digit blocks of~$y$, then call the \texttt{auxii} +% auxiliary to evaluate $y'^{2} = (y+z)^{2}$. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_auxviii_o:nnnnnnn #1#2 #3#4#5#6#7 + { + \exp_after:wN \@@_sqrt_auxix_o:wnwnw + \int_use:N \__int_eval:w #3 + \exp_after:wN \@@_basics_pack_low:NNNNNw + \int_use:N \__int_eval:w #1 + 1#4#5 + \exp_after:wN \@@_basics_pack_low:NNNNNw + \int_use:N \__int_eval:w #2 + 1#6#7 ; + } +\cs_new:Npn \@@_sqrt_auxix_o:wnwnw #1; #2#3; #4#5; + { + \@@_sqrt_auxii_o:NnnnnnnnN + \@@_sqrt_auxiii_o:wnnnnnnnn {#1}{#2}{#3}{#4}{#5} + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP] +% {\@@_sqrt_auxx_o:Nnnnnnnn, \@@_sqrt_auxxi_o:wwnnN} +% At this stage, $j=6$ and $10^{24} z < 10^{7}$, hence +% \[ +% 10^{7} + 1/2 > 10^{24} z + 1/2 \geq +% \bigl(10^{24}(a-y^2) - 258\bigr) \cdot (0.5\cdot 10^{8}) +% \Bigm/ (10^{8} y + 1) \,, +% \] +% then $10^{24}(a-y^2) - 258 < 2 (10^{7} + 1/2) (y + 10^{-8})$, and +% \[ +% 10^{24}(a-y^2) +% < (10^{7} + 1290.5) (1 + 10^{-8}/y) (2y) +% < (10^{7} + 1290.5) (1 + 10^{-7}) (y + \sqrt{a}) \,, +% \] +% which finally implies $0\leq\sqrt{a}-y < 0.2\cdot 10^{-16}$. In +% particular, $y$~is an underestimate of~$\sqrt{a}$ and $y+0.5\cdot +% 10^{-16}$ is a (strict) overestimate. There is at exactly one +% multiple $m$~of $0.5\cdot 10^{-16}$ in the interval $[y, y+0.5\cdot +% 10^{-16})$. If $m^2>a$, then the square root is inexact and is +% obtained by rounding $m-\epsilon$ to a multiple of $10^{-16}$ (the +% precise shift $0<\epsilon<0.5\cdot 10^{-16}$ is irrelevant for +% rounding). If $m^2=a$ then the square root is exactly~$m$, and +% there is no rounding. If $m^2<a$ then we round $m+\epsilon$. For +% now, discard a few irrelevant arguments |#1|, |#2|, |#3|, and find +% the multiple of $0.5\cdot 10^{-16}$ within $[y, y+0.5\cdot +% 10^{-16})$; rather, only the last $4$~digits |#8| of~$y$ are +% considered, and we do not perform any carry yet. The \texttt{auxxi} +% auxiliary sets up \texttt{auxii} with a continuation function +% \texttt{auxxii} instead of \texttt{auxiii} as before. To prevent +% \texttt{auxii} from giving a negative results $a-m^2$, we compute +% $a+10^{-16}-m^2$ instead, always positive since $m<\sqrt{a}+0.5\cdot +% 10^{-16}$ and $a\leq 1-10^{-16}$. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_auxx_o:Nnnnnnnn #1#2#3 #4#5#6#7#8 + { + \exp_after:wN \@@_sqrt_auxxi_o:wwnnN + \int_use:N \__int_eval:w + (#8 + 2499) / 5000 * 5000 ; + {#4} {#5} {#6} {#7} ; + } +\cs_new:Npn \@@_sqrt_auxxi_o:wwnnN #1; #2; #3#4#5 + { + \@@_sqrt_auxii_o:NnnnnnnnN + \@@_sqrt_auxxii_o:nnnnnnnnw + #2 {#1} + {#3} { #4 + \c_one } #5 + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP] +% {\@@_sqrt_auxxii_o:nnnnnnnnw, \@@_sqrt_auxxiii_o:w} +% The difference $0\leq a+10^{-16}-m^2\leq +% 10^{-16}+(\sqrt{a}-m)(\sqrt{a}+m)\leq 2\cdot 10^{-16}$ was just +% computed: its first $8$~digits vanish, as do the next four,~|#1|, +% and most of the following four,~|#2|. The guess~$m$ is an +% overestimate if $a+10^{-16}-m^2 < 10^{-16}$, that is, |#1#2| +% vanishes. Otherwise it is an underestimate, unless +% $a+10^{-16}-m^2=10^{-16}$ exactly. For an underestimate, call the +% \texttt{auxxiv} function with argument~$9998$. For an exact result +% call it with~$9999$, and for an overestimate call it with~$10000$. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_auxxii_o:nnnnnnnnw 0; #1#2#3#4#5#6#7#8 #9; + { + \if_int_compare:w #1#2 > \c_zero + \if_int_compare:w #1#2 = \c_one + \if_int_compare:w #3#4 = \c_zero + \if_int_compare:w #5#6 = \c_zero + \if_int_compare:w #7#8 = \c_zero + \@@_sqrt_auxxiii_o:w + \fi: + \fi: + \fi: + \fi: + \exp_after:wN \@@_sqrt_auxxiv_o:wnnnnnnnN + \__int_value:w 9998 + \else: + \exp_after:wN \@@_sqrt_auxxiv_o:wnnnnnnnN + \__int_value:w 10000 + \fi: + ; + } +\cs_new:Npn \@@_sqrt_auxxiii_o:w \fi: \fi: \fi: \fi: #1 \fi: ; + { + \fi: \fi: \fi: \fi: \fi: + \@@_sqrt_auxxiv_o:wnnnnnnnN 9999 ; + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[aux, rEXP]{\@@_sqrt_auxxiv_o:wnnnnnnnN} +% This receives $9998$, $9999$ or $10000$ as~|#1| when $m$~is an +% underestimate, exact, or an overestimate, respectively. Then +% comes~$m$ as five blocks of~$4$ digits, but where the last +% block~|#6| may be $0$, $5000$, or~$10000$. In the latter case, we +% need to add a carry, unless $m$~is an overestimate (|#1|~is then +% $10000$). Then comes~$a$ as three arguments. Rounding is done by +% \cs{@@_round:NNN}, whose first argument is the final sign~$0$ +% (square roots are positive). We fake its second argument. It +% should be the last digit kept, but this is only used when ties are +% \enquote{rounded to even}, and only when the result is exactly +% half-way between two representable numbers rational square roots of +% numbers with $16$~significant digits have: this situation never +% arises for the square root, as any exact square root of a $16$~digit +% number has at most $8$~significant digits. Finally, the last +% argument is the next digit, possibly shifted by~$1$ when there are +% further nonzero digits. This is achieved by \cs{@@_round_digit:Nw}, +% which receives (after removal of the $10000$'s digit) one of $0000$, +% $0001$, $4999$, $5000$, $5001$, or~$9999$, which it converts to $0$, +% $1$, $4$, $5$, $6$, and~$9$, respectively. +% \begin{macrocode} +\cs_new:Npn \@@_sqrt_auxxiv_o:wnnnnnnnN #1; #2#3#4#5#6 #7#8#9 + { + \exp_after:wN \@@_basics_pack_high:NNNNNw + \int_use:N \__int_eval:w 1 0000 0000 + #2#3 + \exp_after:wN \@@_basics_pack_low:NNNNNw + \int_use:N \__int_eval:w 1 0000 0000 + + #4#5 + \if_int_compare:w #6 > #1 \exp_stop_f: + \c_one \fi: + + \exp_after:wN \@@_round:NNN + \exp_after:wN 0 + \exp_after:wN 0 + \__int_value:w + \exp_after:wN \use_i:nn + \exp_after:wN \@@_round_digit:Nw + \int_use:N \__int_eval:w #6 + 19999 - #1 ; + \exp_after:wN ; + } +% \end{macrocode} +% \end{macro} +% +% \subsection{Setting the sign} % % \begin{macro}[int, EXP]{\@@_set_sign_o:w} % This function is used for the unary minus and for \texttt{abs}. It |