% \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