summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx')
-rw-r--r--Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx499
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