path: root/Master/texmf-dist/source/latex/l3kernel/l3fp-basics.dtx

% \iffalse meta-comment
%
%% File: l3fp-basics.dtx Copyright (C) 2011-2012 The LaTeX3 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
%%
%%    http://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 released version of this bundle is available from CTAN.
%%
%% -----------------------------------------------------------------------
%%
%% The development version of the bundle can be found at
%%
%%    http://www.latex-project.org/svnroot/experimental/trunk/
%%
%% for those people who are interested.
%%
%%%%%%%%%%%
%% NOTE: %%
%%%%%%%%%%%
%%
%%   Snapshots taken from the repository represent work in progress and may
%%   not work or may contain conflicting material!  We therefore ask
%%   people _not_ to put them into distributions, archives, etc. without
%%   prior consultation with the LaTeX Project Team.
%%
%% -----------------------------------------------------------------------
%%
%
%<*driver>
\RequirePackage{l3names}
\GetIdInfo$Id: l3fp-basics.dtx 3986 2012-07-15 19:23:51Z joseph $
  {L3 Floating-point arithmetic}
\documentclass[full]{l3doc}
\begin{document}
  \DocInput{\jobname.dtx}
\end{document}
%</driver>
% \fi
%
% \title{The \textsf{l3fp-basics} package\thanks{This file
%         has version number \ExplFileVersion, last
%         revised \ExplFileDate.}\\
% Floating point arithmetic}
% \author{^^A
%  The \LaTeX3 Project\thanks
%    {^^A
%      E-mail:
%        \href{mailto:latex-team@latex-project.org}
%          {latex-team@latex-project.org}^^A
%    }^^A
% }
% \date{Released \ExplFileDate}
%
% \maketitle
%
% \begin{documentation}
%
% \end{documentation}
%
% \begin{implementation}
%
% \section{\pkg{l3fp-basics} Implementation}
%
%    \begin{macrocode}
%<*initex|package>
%    \end{macrocode}
%
%    \begin{macrocode}
%<@@=fp>
%    \end{macrocode}
%
% All operations implemented in this module yield the outcome of
% rounding the infinitely precise result of the operation to the
% nearest representable number.
%
% ^^A begin[todo]: move
% \section{Internal storage of floating points numbers}
%
% A floating point number \meta{X} is stored as
% \begin{quote}
%   \cs{s_@@} \cs{@@_chk:w} \meta{case} \meta{sign} \meta{body} |;|
% \end{quote}
% Here, \meta{case} is 0 for $\pm 0$, 1 for normal numbers, 2 for $\pm
% \infty$, and 3 for \texttt{nan}, and \meta{sign} is $0$ for positive
% numbers, $1$ for \texttt{nan}s, and $2$ for negative numbers. The
% \meta{body} of normal numbers is \Arg{exponent} \Arg{X_1} \Arg{X_2}
% \Arg{X_3} \Arg{X_4}, with
% \[
% \meta{X} = (-1)^{\meta{sign}} 10^{-\meta{exponent}} \sum_i
% \meta{X_i} 10^{-4i}.
% \]
% Calculations are done in base $10000$, \emph{i.e.} one myriad.  The
% \meta{exponent} lies between $\pm\cs{c_@@_max_exponent_int} = \pm
% \the\csname\detokenize{c__fp_max_exponent_int}\endcsname$ inclusive.
%
% Additionally, positive and negative floating point numbers may only be
% stored with $1000\leq\meta{X_1}<10000$. This requirement is necessary
% in order to preserve accuracy and speed.
%
% ^^A end[todo]
%
% ^^A begin[todo]
%
% Some algorithms used below end up being quite similar to some
% described in \enquote{What Every Computer Scientist Should Know About
%   Floating Point Arithmetic}, by David Goldberg, which can be found at
% \texttt{http://cr.yp.to/2005-590/goldberg.pdf}. I need to compare them
% very carefully.
%
% ^^A end[todo]
%
%^^A todo sanitize, pack.
%
% \subsection{Common to several operations}
%
% \begin{macro}[EXP]
%   {
%     \@@_basics_pack_low:NNNNNw  ,
%     \@@_basics_pack_high:NNNNNw ,
%     \@@_basics_pack_high_carry:w
%   }
%   Addition and multiplication of mantissas are done in two steps:
%   first compute a (more or less) exact result,
%   then round and pack digits in the final (braced) form.
%   These functions take care of the packing, with special attention
%   given to the case where rounding has caused a carry.
%   In \cs{@@_basics_pack_high_carry:w}, |#1| should
%   always be $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} ;
  }
\cs_new:Npn \@@_basics_pack_high:NNNNNw #1 #2#3#4#5 #6;
  {
    \if_meaning:w 2 #1
      \@@_basics_pack_high_carry:w
    \fi:
    ; {#2#3#4#5} {#6}
  }
\cs_new:Npn \@@_basics_pack_high_carry:w \fi: ; #1
  { \fi: + \c_one ; {1000} }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}
%   {
%     \@@_basics_return_nan_nan:NNww ,
%     \@@_basics_return_zero:NNww    ,
%     \@@_basics_return_inf:NNww     ,
%     \@@_basics_return_i:NNNNww     ,
%     \@@_basics_return_ii:NNNNww    ,
%     \@@_basics_return_nan:NNNNww
%   }
%   \begin{syntax}
%     \cs{@@_basics_return_...:NNww}
%     ~~\meta{sign_1} \meta{sign_2} \meta{body_1} |;| \meta{body_2} |;|
%   \end{syntax}
%   Used for binary operations, to return a value for some special
%   cases (common to several operations). All functions expand once
%   after their arguments.
%
%   ^^A todo: redoc, changed.
%   The \texttt{nan_nan} function combines the \texttt{info} fields
%   of the two \texttt{nan}.\footnote{Bruno: check that messages are
%     kept.}
%   The \texttt{zero} and \texttt{inf} functions return $\pm 0$ or
%   $\pm\infty$ with a sign equal to the product of the two signs:
%   three \cs{exp_after:wN} are needed to escape out of the
%   conditional, and expand once after.
%   The \texttt{i} and \texttt{ii} functions return one of their
%   operands and expand after using \cs{@@_exp_after_o:w}.
%   In some cases, this could be optimized, since we know in advance
%   what case of number we have. However, it seems better to keep the
%   number of control sequences low: these functions are called only
%   in special cases anyways, so performance is not an issue.
%    \begin{macrocode}
\cs_new:Npn \@@_basics_return_nan_nan:NNww #1#2 #3; #4;
  { \@@_exp_after_o:w \s_@@ \@@_chk:w 3 1 #3 ; }
\cs_new:Npn \@@_basics_return_zero:NNww #1#2 #3; #4;
  {
    \if_meaning:w #1 #2
      \exp_after:wN \exp_after:wN \exp_after:wN \c_zero_fp
    \else:
      \exp_after:wN \exp_after:wN \exp_after:wN \c_minus_zero_fp
    \fi:
  }
\cs_new:Npn \@@_basics_return_inf:NNww #1#2 #3; #4;
  {
    \if_meaning:w #1 #2
      \exp_after:wN \exp_after:wN \exp_after:wN \c_inf_fp
    \else:
      \exp_after:wN \exp_after:wN \exp_after:wN \c_minus_inf_fp
    \fi:
  }
\cs_new:Npn \@@_basics_return_i:NNNNww #1#2 #3#4 #5; #6;
  { \@@_exp_after_o:w \s_@@ \@@_chk:w #1 #3 #5; }
\cs_new:Npn \@@_basics_return_ii:NNNNww #1#2 #3#4 #5; #6;
  { \@@_exp_after_o:w \s_@@ \@@_chk:w #2 #4 #6; }
\cs_new:Npn \@@_basics_return_nan:NNww #1#2
  {
    \if_meaning:w 1 #1
      \exp_after:wN \@@_basics_return_i:NNNNww
    \else:
      \exp_after:wN \@@_basics_return_ii:NNNNww
    \fi:
    3 3 #1 #2
  }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Addition and subtraction}
%
% \begin{macro}[int, EXP]{\@@_+_o:ww}
%   For addition, everything is easy.  No need to grab the
%   \meta{body_2}.
%    \begin{macrocode}
\cs_new:cpn { @@_+_o:ww }
    \s_@@ \@@_chk:w #1 #2 #3 ; \s_@@ \@@_chk:w #4 #5
  { \@@_add_cases:NN #1 #4 #2 #5 #3 ; }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[int, EXP]{\@@_-_o:ww}
%   Change the sign of the second argument.
%    \begin{macrocode}
\cs_new:cpn { @@_-_o:ww }
    \s_@@ \@@_chk:w #1 #2 #3 ; \s_@@ \@@_chk:w #4 #5
  {
    \exp_after:wN \@@_add_cases:NN
    \exp_after:wN #1
    \exp_after:wN #4
    \exp_after:wN #2
    \int_use:N \__int_eval:w \c_two - #5 \__int_eval_end:
    #3 ;
  }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Sign, exponent, and special numbers}
%
% \begin{macro}[EXP,aux]{\@@_add_cases:NN}
%   \begin{syntax}
%     \cs{@@_add_cases:NN} \meta{case_1} \meta{case_2}
%     ~~\meta{sign_1} \meta{sign_2} \meta{body_1} |;| \meta{body_2} |;|
%   \end{syntax}
%   This performs the addition. it also expands the following tokens
%   on the input stream once.
%
%   Whenever \meta{case_1} is different from \meta{case_2}, the result
%   is simply the floating point number with the highest \meta{case}.
%   For instance, adding a normal number to a zero gives the normal
%   number, and adding a \texttt{nan} to any non-\texttt{nan} gives
%   that \texttt{nan}. Optimizing for addition of normal numbers,
%   we test for equality and then separate the \enquote{greater than}
%   and \enquote{less than} branches.
%    \begin{macrocode}
\cs_new:Npn \@@_add_cases:NN #1 #2
  {
    \if_int_compare:w #1 = #2 \exp_stop_f:
      \exp_after:wN \@@_add_cases_eq:N
    \else:
      \if_int_compare:w #1 < #2 \exp_stop_f:
        \exp_after:wN \exp_after:wN
        \exp_after:wN \@@_basics_return_ii:NNNNww
      \else:
        \exp_after:wN \exp_after:wN
        \exp_after:wN \@@_basics_return_i:NNNNww
      \fi:
      \exp_after:wN #1
    \fi:
    #2
  }
%    \end{macrocode}
%   If the first \meta{case} is larger, then the first number remains
%   untouched, while the second number is ignored. On the other hand,
%   if the second \meta{case} is larger, the opposite happens: we retain
%   the second number. In both cases, there needs to be one step of
%   expansion after.
%    \begin{macrocode}
%    \end{macrocode}
%   We are then ready for the equality case: we split according
%   to the \meta{case}.
%    \begin{macrocode}
\cs_new:Npn \@@_add_cases_eq:N #1
  {
    \if_case:w #1 \exp_stop_f:
         \exp_after:wN \@@_add_zeros:NNww
    \or: \exp_after:wN \@@_add_normal:NNww
    \or: \exp_after:wN \@@_add_inf:NNww
    \or: \exp_after:wN \@@_basics_return_nan_nan:NNww
    \fi:
  }
%    \end{macrocode}
%   Adding two zeros yields \cs{c_zero_fp}, except if both
%   zeros were $-0$.\footnote{Bruno: this should depend on the
%     rounding mode.}
%    \begin{macrocode}
\cs_new:Npn \@@_add_zeros:NNww #1#2 #3;
  {
    \if_int_compare:w #1 #2 = 02 \exp_stop_f:
      \@@_case_return_o:Nw \c_zero_fp
    \else:
      \@@_case_return_same_o:w
    \fi:
    \s_@@ \@@_chk:w 0 #2
  }
%    \end{macrocode}
%   If both infinities have the same sign, just return that infinity,
%   otherwise, it is an invalid operation.
%    \begin{macrocode}
\cs_new:Npn \@@_add_inf:NNww #1#2 #3;
  {
    \if_meaning:w #1 #2
      \@@_case_return_same_o:w
    \else:
      \@@_case_use:nw
        {
          \@@_invalid_operation:Nnww \c_nan_fp { + }
            \s_@@ \@@_chk:w 2 #1 #3 ;
        }
    \fi:
    \s_@@ \@@_chk:w 2 #2
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_add_normal:NNww}
%   \begin{syntax}
%     \cs{@@_add_normal:NNww} \meta{sign_1} \meta{sign_2}
%     ~~\Arg{exp_1}  \meta{body_1} |;| \Arg{exp_2} \meta{body_2} |;|
%   \end{syntax}
%   We now have two normal numbers to add, and we have to check signs
%   and exponents more carefully before performing the addition.
%    \begin{macrocode}
\cs_new:Npn \@@_add_normal:NNww #1#2
  {
    \if_meaning:w #1#2
      \exp_after:wN \@@_add_npos:Nnwnw
    \else:
      \exp_after:wN \@@_sub_npos:Nnwnw
    \fi:
    #1
  }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Absolute addition}
%
% In this subsection, we perform the addition
% of two positive normal numbers.
%
% \begin{macro}[EXP]{\@@_add_npos:Nnwnw}
%   \begin{syntax}
%     \cs{@@_add_npos:Nnwnw} \meta{sign}
%     ~~\Arg{exp_1}  \meta{body_1} |;| \Arg{exp_2} \meta{body_2} |;|
%   \end{syntax}
%   Since we are doing an addition, \meta{sign} will be the final sign.
%   The only special case which may arise is the case of an overflow.
%   This will be checked by \cs{@@_sanitize:Nw} at the end of
%   the calculation. We start an \cs{__int_eval:w}, responsible for
%   computing the exponent, which may receive a contribution of |+1|
%   in case of carry. The exponent should be stopped by |;| followed by
%   the overall \meta{sign} for the sanitizing to work properly.
%
%   Grab and compare the exponents. The smaller number is decimated until
%   its exponent reaches that of the bigger number. We need to bring the
%   final sign down in the midst of the calculation to do the rounding
%   correctly.
%    \begin{macrocode}
\cs_new:Npn \@@_add_npos:Nnwnw #1 #2#3; #4
  {
    \exp_after:wN \@@_sanitize:Nw
    \exp_after:wN #1
    \int_use:N \__int_eval:w
      \if_int_compare:w #2 > #4 \exp_stop_f:
        #2
        \exp_after:wN \@@_add_big_i:wNww \__int_value:w -
      \else:
        #4
        \exp_after:wN \@@_add_big_ii:wNww \__int_value:w
      \fi:
      \__int_eval:w #4 - #2 ; #1 #3;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_add_big_i:wNww,
%     \@@_add_big_ii:wNww}
%   \begin{syntax}
%     \cs{@@_add_big_i:wNww} \meta{shift} |;| \meta{sign}
%     ~~\meta{body_1} |;| \meta{body_2} |;|
%   \end{syntax}
%   Shift the mantissa of the small number, and then add with
%   \cs{@@_add_mantissa:NnnwnnnnN}.
%    \begin{macrocode}
\cs_new:Npn \@@_add_big_i:wNww #1; #2 #3; #4;
  {
    \@@_decimate:nNnnnn {#1}
      \@@_add_mantissa:NnnwnnnnN
      #4
    #3
    #2
  }
\cs_new:Npn \@@_add_big_ii:wNww #1; #2 #3; #4;
  {
    \@@_decimate:nNnnnn {#1}
      \@@_add_mantissa:NnnwnnnnN
      #3
    #4
    #2
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_add_mantissa:NnnwnnnnN}
%   \begin{syntax}
%     \cs{@@_add_mantissa:NnnwnnnnN}
%     ~~\meta{rounding}
%     ~~\Arg{Y'_1} \Arg{Y'_2} \meta{extra-digits} |;|
%     ~~\Arg{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4}
%     ~~\meta{final sign}
%   \end{syntax}
%   To round properly, we must know at which digit the rounding
%   should occur. This requires to know whether the addition
%   produces an overall carry or not. Thus, we do the computation
%   now and check for a carry, then go back and do the rounding.
%   The rounding may cause a carry in very rare cases such as
%   $0.99\cdots 95 \to 1.00\cdots 0$, but this situation always
%   give an exact power of $10$, for which it is easy to correct
%   the result at the end.
%    \begin{macrocode}
\cs_new:Npn \@@_add_mantissa:NnnwnnnnN #1 #2#3 #4; #5#6#7#8
  {
    \exp_after:wN \@@_add_mantissa_test:N
    \int_use:N \__int_eval:w 1#5#6 + #2
      \exp_after:wN \@@_add_mantissa_pack:NNNNNNN
      \int_use:N \__int_eval:w 1#7#8 + #3 ; #1
  }
%    \end{macrocode}
%
%    \begin{macrocode}
\cs_new:Npn \@@_add_mantissa_pack:NNNNNNN #1 #2#3#4#5#6#7
  {
    \if:w 2 #1
      + \c_one
    \fi:
    ; #2 #3 #4 #5 #6 #7 ;
  }
\cs_new:Npn \@@_add_mantissa_test:N #1
  {
    \if:w 2 #1
      \exp_after:wN \@@_add_mantissa_carry:wwNNNN
    \else:
      \exp_after:wN \@@_add_mantissa_no_carry:wwNNNN
    \fi:
  }
%    \end{macrocode}
%
%   \begin{quote}
%     \cs{@@_add_mantissa_no_carry:wwNNNN}
%     ~~\meta{8d} |;| \meta{6d} |;| \meta{2d} |;|
%     ~~\meta{rounding} \meta{sign}
%   \end{quote}
%   If there's no carry, grab all the digits again, and just
%   set the rounding correctly.\footnote{Bruno: an optimization
%     would be to compute whether we need rounding or not,
%     and only grab digits if there is rounding.}
%
%    \begin{macrocode}
\cs_new:Npn \@@_add_mantissa_no_carry:wwNNNN
    #1; #2; #3#4 ; #5#6
  {
    \exp_after:wN \@@_basics_pack_high:NNNNNw
    \int_use:N \__int_eval:w 1 #1
      \exp_after:wN \@@_basics_pack_low:NNNNNw
      \int_use:N \__int_eval:w 1 #2 #3#4
        + \@@_round:NNN #6 #4 #5
        \exp_after:wN ;
  }
%    \end{macrocode}
%
%   The case where there is a carry is very similar: rounding can even
%   raise the first digit from $1$ to $2$ (but we don't need to check that).
%   \begin{quote}
%     \cs{@@_add_mantissa_carry:wwNNNN}
%     ~~\meta{8d} |;| \meta{6d} |;| \meta{2d} |;|
%     ~~\meta{rounding} \meta{sign}
%   \end{quote}
%    \begin{macrocode}
\cs_new:Npn \@@_add_mantissa_carry:wwNNNN
    #1; #2; #3#4; #5#6
  {
    + \c_one
    \exp_after:wN \@@_add_mantissa_carry_pack:NNNNNNNNw
    \int_use:N \__int_eval:w 1 #1
      \exp_after:wN \@@_add_mantissa_carry_pack_ii:NNNNw
      \int_use:N \__int_eval:w 1 #2#3
        + \@@_round:NNNN #6 #3 #4 #5
        \exp_after:wN ;
  }
\cs_new:Npn \@@_add_mantissa_carry_pack_ii:NNNNw #1 #2#3#4 #5;
  {
    \if:w 2 #1
      + \c_one
    \fi:
    \__int_eval_end:
    #2#3#4; {#5} ;
  }
\cs_new:Npn \@@_add_mantissa_carry_pack:NNNNNNNNw
    #1#2#3#4 #5#6#7#8 #9; { ; {#1#2#3#4} {#5#6#7#8} {#9} }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Absolute subtraction}
%
% \begin{macro}[EXP]{\@@_sub_npos:Nnwnw}
%   \begin{syntax}
%     \cs{@@_sub_npos:Nnwnw} \meta{sign}
%     ~~\Arg{exp_1}  \meta{body_1} |;| \Arg{exp_2} \meta{body_2} |;|
%   \end{syntax}
%   Rounding properly in some modes requires to know what the sign
%   of the result will be. For addition, this was easy. Here, besides
%   comparing the exponents to know how to decimate, we need to
%   check carefully which number is bigger when they have the same
%   exponent.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_npos:Nnwnw #1 #2#3; #4 #5;
  {
    \exp_after:wN \@@_sanitize:wN
    \int_use:N \__int_eval:w
      \if_int_compare:w #2 > #4 \exp_stop_f:
        #2
        \exp_after:wN \@@_sub_big_i:wNww \__int_value:w -
      \else:
        #4
        \if_int_compare:w #2 = #4 \exp_stop_f:
          \@@_sub_exponent_eq:nnnnnnnn #3 #5
        \else:
          \exp_after:wN \@@_sub_big_ii:wNww \__int_value:w
        \fi:
      \fi:
      \__int_eval:w #4 - #2 ; #1 #3; #5;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_sub_exponent_eq:nnnnnnnn}
%    \begin{macrocode}
\cs_new:Npn \@@_sub_exponent_eq:nnnnnnnn #1#2#3#4 #5#6#7#8
  {
    \if_int_compare:w #1#2 > #5#6 \exp_stop_f:
      \exp_after:wN \@@_sub_big_i:wNww \__int_value:w
    \else:
      \if_int_compare:w #1#2 < #5#6 \exp_stop_f:
        \exp_after:wN \@@_sub_big_ii:wNww \__int_value:w
      \else:
        \if_int_compare:w #3#4 > #7#8 \exp_stop_f:
          \exp_after:wN \@@_sub_big_i:wNww \__int_value:w
        \else:
          \if_int_compare:w #3#4 < #7#8 \exp_stop_f:
            \exp_after:wN \@@_sub_big_ii:wNww \__int_value:w
          \else:
            \exp_after:wN \@@_sub_eq:wNww \__int_value:w
          \fi:
        \fi:
      \fi:
    \fi:
  }
\cs_new:Npn \@@_sub_eq:wNww #1; #2 #3; #4;
  { \exp_after:wN ; \exp_after:wN 1 \exp_after:wN ; }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_sub_big_i:wNww,\@@_sub_big_ii:wNww}
%   \begin{syntax}
%     \cs{@@_sub_big_i:wNww} \meta{shift} |;| \meta{sign}
%     ~~\meta{body_1} |;| \meta{body_2} |;|
%   \end{syntax}
%   Shift the mantissa of the small number, and then subtract with
%   \cs{@@_sub_back_mantissa:NnnwNnnnn}.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_big_i:wNww #1; #2 #3; #4;
  {
    \@@_decimate:nNnnnn {#1}
      \@@_sub_back_mantissa:NnnwNnnnn
      #4
    #2
    #3
  }
\cs_new:Npn \@@_sub_big_ii:wNww #1; #2 #3; #4;
  {
    \exp_after:wN \@@_sub_big_i:wNww
    \__int_value:w #1 \exp_after:wN ;
    \int_use:N \__int_eval:w 2 - #2 \__int_eval_end:
    #4; #3;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_sub_back_mantissa:NnnwNnnnn}
%   \begin{syntax}
%     \cs{@@_sub_back_mantissa:NnnwNnnnn}
%     ~~\meta{rounding} \Arg{Y'_1} \Arg{Y'_2} \meta{extra-digits}
%     ~~\meta{final sign}
%     ~~\Arg{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4}
%   \end{syntax}
%   At this stage, we know that \meta{Y} is less than \meta{X},
%   and we know the final sign.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_back_mantissa:NnnwNnnnn #1 #2#3 #4; #5 #6#7#8#9
  {
    \exp_after:wN \@@_sub_back_mantissa_i:NNwNNNNwN
    \exp_after:wN #1
    \exp_after:wN #5
    \int_use:N \__int_eval:w 2#6#7 - #2 - \c_two +
      \exp_after:wN \@@_sub_back_mantissa_round:wNN
      \int_use:N \__int_eval:w 2#8#9 - #3 ; #1 #5
  }
%    \end{macrocode}
%   After the computation, we need to check whether the first digit of
%   the result is zero. This can only happen if the two numbers had the
%   same exponent, or exponents differing by $1$. In the latter case,
%   the \meta{rounding} digit is not quite enough to let us retrieve
%   the exact result (consider $\cdots25$ and $\cdots15$, both rounded
%   to $\cdots2$ in the usual mode), so we also move the result of
%   \cs{@@_round_neg:NNN} upstream as the digit $0$ or $1$.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_back_mantissa_round:wNN #1; #2 #3
  {
    \exp_after:wN \@@_sub_back_mantissa_iii:N
    \__int_value:w
        \exp_after:wN \@@_round_neg:NNN
        \exp_after:wN #3
        \use_none:nnnnnnnn #1 #2
      + #1
    \exp_after:wN ;
  }
\cs_new:Npn \@@_sub_back_mantissa_iii:N #1
  {
    \exp_after:wN \@@_sub_back_mantissa_ii:NNNNNNw
    \exp_after:wN #1
    \int_use:N \__int_eval:w
      - #1
  }
\cs_new:Npn \@@_sub_back_mantissa_ii:NNNNNNw #1 #2 #3#4#5#6 #7;
  { #2 ; #1 {#3#4#5#6} {#7} ; }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_sub_back_mantissa_i:NNwNNNNwN}
%   Here, |#3| should always be $2$, but we have to take is
%   as a normal undelimited argument since that would break
%   if |#2| is $2$.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_back_mantissa_i:NNwNNNNwN #1#2 #3 #4#5#6#7 #8; #9
  {
    \if:w 0 #4
      \exp_after:wN \@@_sub_back_carry:NNwNnnnn
      \exp_after:wN #1
      \exp_after:wN #9
    \fi:
    ; #2
    {#4#5#6#7} {#8}
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_sub_back_carry:NNwNnnnn}
%   \begin{syntax}
%     \cs{@@_sub_back_carry:NNwNnnnn}
%     ~~\meta{rounding} \meta{0 or 1} |;| \meta{final sign}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} |;|
%   \end{syntax}
%   This function is called when $\meta{Z_1}\leq 999$. We revert
%   the carry, which is given by \meta{0 or 1}, and subtract the
%   \meta{rounding} digit as appropriate, then feed the result,
%   of the form \meta{$\leq$ 7d} |;| \meta{9d} |;| to
%   \cs{@@_sub_back_carry_i:wwN}. The result is always exact.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_back_carry:NNwNnnnn #1#2 ; #3 #4#5#6#7 ;
  {
    \exp_after:wN \@@_sub_back_carry_i:wwN
    \int_use:N \__int_eval:w #4 #5 - 1 + \exp_after:wN \@@_use_s:n
    \int_use:N \__int_eval:w 1 #6 #7 0 + #2 0 - #1 ; #3
  }
%    \end{macrocode}
%   Unless the first block is zero, check how many digits is has,
%   and shift the exponent down by the corresponding amount. Then
%   pack digits into blocks of $4$ (there are between $10$ and $16$
%   digits in front of \cs{@@_sub_back_carry_large:NNNNNNNNw}).
%    \begin{macrocode}
\cs_new:Npn \@@_sub_back_carry_i:wwN #1 ;
  {
    \if:w 0 #1
      - 8
      \exp_after:wN \@@_sub_back_carry_small:wN \__int_value:w
    \else:
      - \@@_sub_back_carry_ii:NNNNNNNNw #1 1234567;
      \exp_after:wN \@@_sub_back_carry_large:NNNNNNNNw
    \fi:
    #1
  }
%    \end{macrocode}
%   The case where the number is non-zero is slightly easier.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_back_carry_ii:NNNNNNNNw #1#2#3#4#5#6#7#8#9; {#8}
\cs_new:Npn \@@_sub_back_carry_large:NNNNNNNNw #1#2#3#4 #5#6#7#8 #9;
  {
    \@@_sub_back_carry_large_ii:NNNNNNNNw
    #9 000000 ; {#1#2#3#4} {#5#6#7#8}
  }
\cs_new:Npn \@@_sub_back_carry_large_ii:NNNNNNNNw #1#2#3#4 #5#6#7#8 #9;
  { \@@_sub_back_carry_large_iii:nnnnN {#1#2#3#4} {#5#6#7#8} }
\cs_new:Npn \@@_sub_back_carry_large_iii:nnnnN #1#2 #3#4 #5
  { ; #5 {#3}{#4} {#1}{#2} ; }
%    \end{macrocode}
%   In the case of a \enquote{small} result, what comes after
%   \cs{@@_sub_back_carry_small:wN} has between $1$
%   and $9$ digits, and is not zero.
%    \begin{macrocode}
\cs_new:Npn \@@_sub_back_carry_small:wN #1;
  {
    - \exp_after:wN \@@_use_i_until_s:nw
      \use_none:nnnnnnnnn #1 012345678;
    \@@_sub_back_carry_small_ii:NNNNNNNN #1 00000000 ;
  }
\cs_new:Npn \@@_sub_back_carry_small_ii:NNNNNNNN #1#2#3#4 #5#6#7#8
  { \@@_sub_back_carry_small_iii:nnNwN {#1#2#3#4} {#5#6#7#8} }
\cs_new:Npn \@@_sub_back_carry_small_iii:nnNwN #1 #2 #3 #4; #5
  { ; #5 {#1} {#2} {#3000} {0000} ; }
%    \end{macrocode}
% \end{macro}
%
%
% \subsection{Multiplication}
%
% \begin{macro}[int, EXP]{\@@_*_o:ww}
%   For multiplication, everything is easy.  No need to grab the
%   \meta{body_2}.
%    \begin{macrocode}
\cs_new:cpn { @@_*_o:ww }
    \s_@@ \@@_chk:w #1 #2 #3 ; \s_@@ \@@_chk:w #4 #5
  { \@@_mul_cases:NN #1 #4 #2 #5 #3 ; }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Signs, and special numbers}
%
% \begin{macro}[EXP,aux]{\@@_mul_cases:NN}
%   \begin{syntax}
%     \cs{@@_mul_cases:NN} \meta{case_1} \meta{case_2}
%     ~~\meta{sign_1} \meta{sign_2} \meta{body_1} |;| \meta{body_2} |;|
%   \end{syntax}
%   Expands the following tokens on the input stream once.
%   The special cases are coded at the start of this module,
%   and identical to the ones for division.\footnote{Bruno: \texttt{nan}
%     are not treated properly: $\infty\times 0$ should signal.}
%^^A todo: use the faster \if_meaning:w.
%    \begin{macrocode}
\cs_new:Npn \@@_mul_cases:NN #1 #2
  {
    \if_case:w \if_meaning:w 1 #1 #2 \else:
               \if_meaning:w 1 #2 #1 \else:
               \if_meaning:w #1#2 #1 \else:
               \if_int_compare:w \__int_eval:w #1 + #2 > \c_two
                 3 \else: 4 \fi: \fi: \fi: \fi:
               \exp_stop_f:
         \exp_after:wN \@@_basics_return_zero:NNww
    \or: \exp_after:wN \@@_mul_normal:NNww
    \or: \exp_after:wN \@@_basics_return_inf:NNww
    \or: \exp_after:wN \@@_basics_return_nan:NNww
    \or:
      \exp_after:wN \@@_mul_invalid:NNNNww
      \exp_after:wN #1
      \exp_after:wN #2
    \fi:
  }
\cs_new:Npn \@@_mul_invalid:NNNNww #1#2#3#4#5; #6;
  {
    \@@_invalid_operation:Nnww \c_nan_fp { * }
      \s_@@ \@@_chk:w #1 #3 #5 ;
      \s_@@ \@@_chk:w #2 #4 #6 ;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_mul_normal:NNww}
%   \begin{syntax}
%     \cs{@@_mul_normal:NNww} \meta{sign_1} \meta{sign_2}
%     ~~\Arg{exp_1}  \meta{body_1} |;| \Arg{exp_2} \meta{body_2} |;|
%   \end{syntax}
%   We now have two normal numbers to multiply. Combine the signs.
%    \begin{macrocode}
\cs_new:Npn \@@_mul_normal:NNww #1#2
  {
    \if:w #1#2
      \exp_after:wN \@@_mul_npos:Nnwnw
      \exp_after:wN 0
    \else:
      \exp_after:wN \@@_mul_npos:Nnwnw
      \exp_after:wN 2
    \fi:
  }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Absolute multiplication}
%
% In this subsection, we perform the multiplication
% of two positive normal numbers.
%
% \begin{macro}[EXP]{\@@_mul_npos:Nnwnw}
%   \begin{syntax}
%     \cs{@@_mul_npos:Nnwnw} \meta{sign}
%     ~~\Arg{exp_1}  \meta{body_1} |;| \Arg{exp_2} \meta{body_2} |;|
%   \end{syntax}
%   As for addition, \meta{sign} is the final sign. After the computation,
%   \cs{@@_sanitize:Nw} checks for overflow or underflow.
%   As before, \cs{__int_eval:w} computes the exponent, catching any
%   shift coming from the computation in the mantissa. Again, the
%   \meta{sign} is needed for rounding to be done properly, so we move
%   it around with us. We setup the post-expansion here, triggered by
%   \cs{@@_mul_mantissa:nnnnNnnnn}.
%    \begin{macrocode}
\cs_new:Npn \@@_mul_npos:Nnwnw #1 #2#3; #4 #5;
  {
    \exp_after:wN \@@_sanitize:Nw
    \exp_after:wN #1
    \int_use:N \__int_eval:w
      #2 + #4
      \@@_mul_mantissa:nnnnNnnnn #3 #1 #5
    \exp_after:wN ;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_mul_mantissa:nnnnNnnnn}
%   \begin{syntax}
%     \cs{@@_mul_mantissa:nnnnNnnnn}
%     ~~\Arg{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4} \meta{sign}
%     ~~\Arg{Y_1} \Arg{Y_2} \Arg{Y_3} \Arg{Y_4} \meta{;}
%   \end{syntax}
%   After one expansion, the token following \meta{Y_4} must be a semicolon
%   (represented by \meta{;}).
%
%   The product of two $16$ digit integers has $31$ or $32$ digits,
%   but it is impossible to know which one before computing. The place
%   where we round depends on that number of digits, and may depend
%   on all digits until the last in some rare cases. The approach is
%   thus to compute the $5$ first blocks of $4$ digits (the first one
%   is between $100$ and $9999$ inclusive), and a compact version of
%   the remaining $3$ blocks. Afterwards, the number of digits is
%   known, and we can do the rounding within yet another set of
%   \cs{__int_eval:w}.
%    \begin{macrocode}
\cs_new:Npn \@@_mul_mantissa:nnnnNnnnn #1#2#3#4 #5 #6#7#8#9
  {
    \exp_after:wN \@@_mul_mantissa_after:NNN
    \exp_after:wN #5
    \int_use:N \__int_eval:w 99990000 + #1*#6 +
      \exp_after:wN \@@_mul_mantissa_keep:NNNNNw
      \int_use:N \__int_eval:w 99990000 + #1*#7 + #2*#6 +
        \exp_after:wN \@@_mul_mantissa_keep:NNNNNw
        \int_use:N \__int_eval:w 99990000 + #1*#8 + #2*#7 + #3*#6 +
          \exp_after:wN \@@_mul_mantissa_drop:NNNNNw
          \int_use:N \__int_eval:w  99990000 + #1*#9 + #2*#8 + #3*#7 + #4*#6 +
            \exp_after:wN \@@_mul_mantissa_drop:NNNNNw
            \int_use:N \__int_eval:w  99990000 + #2*#9 + #3*#8 + #4*#7 +
              \exp_after:wN \@@_mul_mantissa_drop:NNNNNw
              \int_use:N \__int_eval:w  99990000 + #3*#9 + #4*#8 +
                \exp_after:wN \@@_mul_mantissa_drop:NNNNNw
                \int_use:N \__int_eval:w 100000000 + #4*#9 \exp_after:wN ;
  }
\cs_new:Npn \@@_mul_mantissa_drop:NNNNNw #1#2#3#4#5 #6;
  { #1#2#3#4#5 ; + #6 }
\cs_new:Npn \@@_mul_mantissa_keep:NNNNNw #1#2#3#4#5 #6;
  { #1#2#3#4#5 ; #6 ; }
%    \end{macrocode}
%   Once the first \cs{int_use:N} \cs{__int_eval:w}, and all the
%   \cs{@@_mul_mantissa_...:NNNNNw} have been expanded,
%   we get
%   \begin{quote}
%     \cs{@@_mul_mantissa_after:NNN} \meta{sign} |1|
%     ~~\meta{digits 1--8} |;| \meta{digits 9--12} |;| \meta{digits 13--16} |;|
%     ~~|+| \meta{digits 17--20} |+| \meta{digits 21--24}
%     ~~|+| \meta{digits 25--28} |+| \meta{digits 29--32} |;|
%   \end{quote}
%   If the \meta{digit 1} is non-zero, then for rounding we only care
%   about the digits $16$ and $17$, and whether all other digits are zero
%   or not (check for exact ties). On the other hand, if \meta{digit 1}
%   is zero, we care about digits $17$ and $18$, and whether all others are
%   zero.
%    \begin{macrocode}
\cs_new:Npn \@@_mul_mantissa_after:NNN #1 #2 #3
  {
    \if:w 0 #3
      \exp_after:wN \@@_mul_mantissa_small:NNwwwN
    \else:
      \exp_after:wN \@@_mul_mantissa_large:NwwNNNN
    \fi:
    #1 #3
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_mul_mantissa_large:NwwNNNN}
%   In this branch, \meta{digit 1} is non-zero. The result is thus
%   \meta{digits 1--16}, plus some rounding which depends on the digits
%   $16$, $17$, and whether all subsequent digits are zero or not.
%   Here, \cs{@@_round_s:NNNw} takes the \meta{sign}, followed by
%   digits $16$, $17$, and an integer expression which is zero if and
%   only if all further digits are zero.
%    \begin{macrocode}
\cs_new:Npn \@@_mul_mantissa_large:NwwNNNN #1 #2; #3; #4#5#6#7; +
  {
    \exp_after:wN \@@_basics_pack_high:NNNNNw
    \int_use:N \__int_eval:w 1#2
      \exp_after:wN \@@_basics_pack_low:NNNNNw
      \int_use:N \__int_eval:w 1#3#4#5#6#7 + \@@_round_s:NNNw #1 #7
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_mul_mantissa_small:NNwwwN}
%   In this branch, \meta{digit 1} is zero. Our result will thus be
%   \meta{digits 2--17}, plus some rounding which depends on the digits
%   $17$, $18$, and whether all subsequent digits are zero or not.
%   The $8$ digits |1#3| are followed, after expansion of the
%   \texttt{small_pack} auxiliary, by the next digit, to form a $9$
%   digit number. Also, rounding may have caused a carry, which is
%   then converted to \cs{c_ten} rather than the usual \cs{c_one},
%   because of the shift.
%    \begin{macrocode}
\cs_new:Npn \@@_mul_mantissa_small:NNwwwN #1 #2#3; #4; #5; + #6
  {
    - \c_one
    \exp_after:wN \@@_basics_pack_high:NNNNNw
    \int_use:N \__int_eval:w 1#3
      \exp_after:wN \@@_mul_mantissa_small_pack:NNNNNNw
      \int_use:N \__int_eval:w 1#4#5#6 + \@@_round_s:NNNw #1 #6
  }
\cs_new:Npn \@@_mul_mantissa_small_pack:NNNNNNw #1#2 #3#4#5#6 #7;
  {
    #2
    \if:w 2 #1
      + \c_ten
    \fi:
    ; {#3#4#5#6} {#7} ;
  }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Division}
%
% Time is now ripe to tackle the hardest of the four elementary
% operations: division.
%
% \begin{macro}[EXP]{\@@_/_o:ww}
%   For division we swap the two floating point numbers.
%    \begin{macrocode}
\cs_new:cpn { @@_/_o:ww }
    \s_@@ \@@_chk:w #1 #2 #3 ; \s_@@ \@@_chk:w #4 #5 #6 ;
  { \@@_div_back_cases:NN #4 #1 #5 #2 #6 ; #3 ; }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Signs, and special numbers}
%
% In the case of division, the order of the operands matters,
% and it turns out to be slightly simpler if we internally
% compute the \enquote{backwards} division.
%
% \begin{macro}[EXP,aux]{\@@_div_back_cases:NN}
%   \begin{syntax}
%     \cs{@@_div_back_cases:NN} \meta{case_2} \meta{case_1}
%     ~~\meta{sign_2} \meta{sign_1} \meta{body_2} |;| \meta{body_1} |;|
%   \end{syntax}
%   Expands the following tokens on the input stream once.
%    \begin{macrocode}
\cs_new:Npn \@@_div_back_cases:NN #1 #2
  {
    \if_case:w \if_int_compare:w #1 = #2 \exp_stop_f:
                 #1 \exp_stop_f:
               \else:
                 \if_int_compare:w #1 < #2 \exp_stop_f:
                   \if:w 3 #2 \c_four \else: \c_five \fi:
                 \else:
                   \if:w 3 #1 \c_six \else: \c_seven \fi:
                 \fi:
               \fi:
         \exp_after:wN \@@_div_back_invalid:NNNww \exp_after:wN 0
    \or: \exp_after:wN \@@_div_back_normal:NNww
    \or: \exp_after:wN \@@_div_back_invalid:NNNww \exp_after:wN 2
    \or: \exp_after:wN \@@_basics_return_nan_nan:NNww
    \or:
      \exp_after:wN \@@_basics_return_ii:NNNNww
      \exp_after:wN #1
      \exp_after:wN #2
    \or: \exp_after:wN \@@_basics_return_inf:NNww
    \or:
      \exp_after:wN \@@_basics_return_i:NNNNww
      \exp_after:wN #1
      \exp_after:wN #2
    \or: \exp_after:wN \@@_basics_return_zero:NNww
    \fi:
  }
%    \end{macrocode}
%   Most of the special cases are common with some
%   previous operations. We only need to write the cases of
%   $0/0$ and $\infty/\infty$.
%    \begin{macrocode}
\cs_new:Npn \@@_div_back_invalid:NNNww #1#2#3 #4; #5;
  {
    \@@_invalid_operation:Nnww \c_nan_fp { / }
      \s_@@ \@@_chk:w #1 #3 #5 ;
      \s_@@ \@@_chk:w #1 #2 #4 ;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_div_back_normal:NNww}
%   \begin{syntax}
%     \cs{@@_div_back_normal:NNww} \meta{sign_1} \meta{sign_2}
%     ~~\Arg{exp_1}  \meta{body_1} |;| \Arg{exp_2} \meta{body_2} |;|
%   \end{syntax}
%   We now have two normal numbers to divide. Combine the signs.
%    \begin{macrocode}
\cs_new:Npn \@@_div_back_normal:NNww #1#2
  {
    \if:w #1#2
      \exp_after:wN \@@_div_back_npos:Nnwnw
      \exp_after:wN 0
    \else:
      \exp_after:wN \@@_div_back_npos:Nnwnw
      \exp_after:wN 2
    \fi:
  }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Absolute (backwards) division}
%
% In this subsection, we perform the division
% of two positive normal numbers.
%
% \begin{macro}[EXP]{\@@_div_back_npos:Nnwnw}
%   \begin{syntax}
%     \cs{@@_div_back_npos:Nnwnw} \meta{sign}
%     ~~\Arg{exp Z} \Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} |;|
%     ~~\Arg{exp A} \Arg{A_1} \Arg{A_2} \Arg{A_3} \Arg{A_4} |;|
%   \end{syntax}
%   We want to compute $A/Z$. As for addition and multiplication,
%   \meta{sign} is the final sign. Checking for underflow and
%   overflow is done using the same auxiliary as for multiplication.
%   As explained just below, we first compute $y$, which is
%   the $5$ first digits of $Z$, plus $1$, and then compute pieces
%   of the quotient roughly $4$ digits at a time. Here, |#1| is
%   a single digit, |#2| and |#7| are the exponents (integers),
%   |#8| is three brace groups, and all others are each $4$ digits.
%    \begin{macrocode}
\cs_new:Npn \@@_div_back_npos:Nnwnw #1 #2 #3#4#5#6; #7 #8;
  {
    \exp_after:wN \@@_sanitize:Nw
    \exp_after:wN #1
    \int_use:N \__int_eval:w
      #7 - #2
      \@@_div_mantissa_i:wNwnn #3; #4;
      #8 {#3}{#4}{#5}{#6} #1
  }
%    \end{macrocode}
% \end{macro}
%
%   We are given two numbers, $A=0.A_{1}A_{2}A_{3}A_{4}$
%   and $Z=0.Z_{1}Z_{2}Z_{3}Z_{4}$, in blocks of $4$ digits,
%   and we know that the first digits of $A_{1}$ and of $Z_{1}$
%   are non-zero. To compute $A/Z$, we proceed as follows.
%   \begin{itemize}
%   \item Find an integer $Q_{1} \simeq 10^4 A / Z$.
%   \item Replace $A$ by $B = 10^4 A - Q_{1} Z$.
%   \item Find an integer $Q_{2} \simeq 10^4 B / Z$.
%   \item Replace $B$ by $C = 10^4 B - Q_{2} Z$.
%   \item Find an integer $Q_{3} \simeq 10^4 C / Z$.
%   \item Replace $C$ by $D = 10^4 C - Q_{3} Z$.
%   \item Find an integer $Q_{4} \simeq 10^4 D / Z$.
%   \item Consider $E = 10^4 D - Q_{4} Z$, and ensure
%     correct rounding.
%   \end{itemize}
%   The calculations of $B$, $C$, $D$, and $E$ can be done
%   exactly with only $16$ (or $17$) digits.
%
%   Unfortunately, things are not as easy as they seem.
%   Firstly, we make sure that all intermediate steps are positive,
%   since negative results would require extra calculations at the end.
%   This requires that $Q_{1} \leq 10^4 A / Z$ etc. A reasonable
%   attempt would be to define $Q_{1}$ as
%   \[
%   \cs{int_eval:n} \left\{
%     \frac{ A_{1} A_{2} }{ Z_{1} + 1 } - 1 \right\}.
%   \]
%   Subtracting $1$ at the end takes care of the fact that e\TeX{}'s
%   \cs{__int_eval:w} rounds instead of truncating. We add $1$ to $Z_{1}$
%   because $ Z_{1} \leq 10^4 Z < Z_{1}+1$ and we need $Q_{1}$
%   to be an underestimate. However, we are now underestimating
%   $Q_{1}$ too much: it can be wrong by up to $100$, for instance
%   when $Z = 0.1$ and $A \simeq 1$. Then $B$ could take values up to
%   $10$ (maybe more), and a few steps down the line, we would run into
%   arithmetic overflow, since \TeX{} can only handle integers less than
%   roughly $2\cdot 10^9$.
%
%   A better formula is to take
%   \[
%   Q_{1} = \cs{int_eval:n} \left\{
%     \frac{ 10 \cdot A_{1} A_{2} }
%       { \left\lfloor 10^{-3} \cdot Z_{1} Z_{2} \right\rfloor + 1 }
%     - 1 \right\}.
%   \]
%   This is always less than $10^9 A / (10^5 Z)$, as we wanted.
%   In words, we take the $5$ first digits of $Z$ into account,
%   and the $8$ first digits of $A$, using $0$ as a $9$-th digit
%   rather than the true digit for efficiency reasons. We shall
%   prove that using this formula to define all the $Q_{i}$
%   leads to no overflow. For convenience, let us denote
%   \[
%   y = \left\lfloor 10^{-3} \cdot Z_{1} Z_{2} \right\rfloor + 1,
%   \]
%   so that, taking into account the fact that e\TeX{} rounds ties
%   away from zero,
%   \[
%   Q_{1} = \left\lfloor A_{1}A_{2}0/y - 1/2 \right\rfloor.
%   \]
%   Note that $10^4<y\leq 10^5$, and $999 \leq Q_{1} \leq 99989$.
%   Also note that this formula does not cause an overflow as long as
%   $A<2.147\cdots$, since the numerator involves an integer slightly
%   smaller than $10^9A$.
%
%   Let us bound $B$:
%   \begin{align*}
%   10^5 B
%   &=
%   A_{1}A_{2}0 + 10 \cdot 0.A_{3}A_{4}
%   - 10 \cdot Z_{1}.Z_{2}Z_{3}Z_{4}
%   \cdot \left\lfloor A_{1}A_{2}0/y - 1/2 \right\rfloor
%   \\
%   &<
%   A_{1}A_{2}0
%   \cdot \left( 1 - 10 \frac{Z_{1}.Z_{2}Z_{3}Z_{4}}{y} \right)
%   + \frac{3}{2} \cdot 10 \cdot Z_{1}.Z_{2}Z_{3}Z_{4} + 10
%   \\
%   &\leq
%   \frac{A_{1}A_{2}0 \cdot (y - 10 Z_{1}.Z_{2}Z_{3}Z_{4})}{y}
%   + \frac{3}{2} y + 10
%   \\
%   &\leq
%   \frac{A_{1}A_{2}0}{y} + \frac{3}{2} y + 10
%   \leq
%   \frac{10^9 A}{y} + 1.6 y
%   \end{align*}
%   At the last step, we hide $10$ into the second term
%   for later convenience. The same reasoning yields\footnote{Bruno:
%     I need to find much better notations. These are not great.}
%   \begin{align*}
%     10^5 B &< 10^9 A/y + 1.6 y, \\
%     10^5 C &< 10^9 B/y + 1.6 y, \\
%     10^5 D &< 10^9 C/y + 1.6 y, \\
%     10^5 E &< 10^9 D/y + 1.6 y. \\
%   \end{align*}
%   The goal is now to prove that none of $B$, $C$, $D$, and $E$
%   can go beyond $2.147\cdots$. Simply bounding each term on the
%   right-hand side separately will not be tight enough: for instance,
%   we would get $10^5 B < 10^5 + 1.6\cdot 10^5 = 2.6 \cdot 10^5$,
%   which is too large.
%
%   Combining the various inequalities together with $A<1$, we get
%   \begin{align*}
%     10^5 B &< 10^9/y + 1.6 y, \\
%     10^5 C &< 10^{13}/y^2 + 1.6 (y + 10^4), \\
%     10^5 D &< 10^{17}/y^3 + 1.6 (y + 10^4 + 10^8/y), \\
%     10^5 E &< 10^{21}/y^4 + 1.6 (y + 10^4 + 10^8/y + 10^{12}/y^2). \\
%   \end{align*}
%   All of those bounds are convex functions of $y$ (since every power
%   of $y$ involved is convex, and hte coefficients are positive), and
%   thus maximal at one of the end-points of the allowed range
%   $10^4<y\leq 10^5$. Thus,
%   \begin{align*}
%     10^5 B &< \mathrm{max} ( 1.16\cdot 10^5, 1.7 \cdot 10^5), \\
%     10^5 C &< \mathrm{max} ( 1.32\cdot 10^5, 1.77 \cdot 10^5), \\
%     10^5 D &< \mathrm{max} ( 1.48\cdot 10^5, 1.777 \cdot 10^5), \\
%     10^5 E &< \mathrm{max} ( 1.64\cdot 10^5, 1.7777 \cdot 10^5). \\
%   \end{align*}
%   All of those bounds are less than $2.147\cdot 10^5$, and
%   we are thus within \TeX{}'s bounds in all cases!\footnote{Bruno:
%     but I need to check this very carefully again.}
%
%   We will later need to have a bound on the $Q_{i}$. Their
%   definitions imply that $Q_{1} < 10^9 A/y - 1/2 < 10^5 A$ and
%   similarly for the other $Q_{i}$. Thus each of them is at most
%   $177770$.
%
%   The last step is to ensure correct rounding. We have
%   \[
%   A/Z = \sum_{i=1}^4 \left(10^{-4i} Q_{i}\right) + 10^{-16} E/Z
%   \]
%   exactly. Furthermore, we know that the result will be between
%   $0.1$ (inclusive) and $10$, so we only need to know the integer
%   part of $E/Z$, and a \enquote{rounding} digit encoding the rest
%   (see maybe addition for an explanation of why). Equivalently,
%   we need to find the integer part of $2E/Z$, and determine whether
%   it was an exact integer or not (this serves to detect ties). Since
%   \[
%   \frac{2E}{Z} = 2\frac{10^5 E}{10^5 Z}
%   \leq 2\frac{10^5 E}{10^4} < 36,
%   \]
%   this integer part is between $0$ and $35$ inclusive. We let
%   e\TeX{} round
%   \[
%   P = \cs{int_eval:n} \left\{
%     2 \frac{E_{1}E_{2}}{Z_{1}Z_{2}} \right\},
%   \]
%   which differs from $2E/Z$ by at most
%   \[
%   \frac{1}{2}
%   + 2 \left\lvert \frac{E}{Z} - \frac{E}{10^{-8} Z_{1}Z_{2}}\right\rvert
%   + 2 \left\lvert \frac{10^8 E - E_{1}E_{2}}{Z_{1}Z_{2}}\right\rvert
%   < 1,
%   \]
%   ($1/2$ comes from e\TeX{}'s rounding) because each absolute value
%   is less than $10^{-7}$. Thus $P$ is either the correct integer part,
%   or an overestimate by $1$ (impossible if $2E/Z$ is an integer). It
%   then suffices to compare $PZ$ with $2E$ to get the integer part of
%   $2E/Z$ and the information of whether it was an exact quotient or not.
%
% \begin{macro}[EXP]{\@@_div_mantissa_i:wNwnn}
%   \begin{syntax}
%     \cs{@@_div_mantissa_i:wNwnn}
%     ~~\meta{Z_1} |;| \meta{Z_2} |;|
%     ~~\Arg{A_1} \Arg{A_2} \Arg{A_3} \Arg{A_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{syntax}
%   First compute $y$ from the first $5$ digits of $Z$, and
%   unbrace \meta{A_1} and \meta{A_2}.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_i:wNwnn #1; #2 #3; #4 #5
  {
    \exp_after:wN \@@_div_mantissa_ii:ww
    \int_use:N \__int_eval:w #1#2 + \c_one ;
    #4 #5 ;
  }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_ii:ww}
%     ~~\meta{y} |;| \meta{A_1} \meta{A_2} |;| \Arg{A_3} \Arg{A_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%   Compute $Q_{1}$ by evaluating $\meta{A_1}\meta{A_2}0/y - 1$.
%   The result will be output to the left, in an \cs{__int_eval:w}
%   which we start now.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_ii:ww #1; #2;
  {
    \exp_after:wN \@@_div_mantissa_iii:www
    \__int_value:w #1 \exp_after:wN ;
    \__int_value:w
      \exp_after:wN \@@_div_mantissa_calc:Nwwnnnnnn
      \int_use:N \__int_eval:w #20/#1 + 999999 ; #2 ;
  }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_calc:Nwwnnnnnn} \meta{$10^6+{}$Q_1} |;|
%     ~~\meta{A_1} \meta{A_2} |;| \Arg{A_3} \Arg{A_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%   The goal here is to expand to
%   \begin{quote}
%     \meta{$10^6+{}$Q_1} |;| \meta{B_1} \meta{B_2} |;| \Arg{B_3} \Arg{B_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%   where $B = 10^4 A - Q_{1} \cdot Z$. More generally, this function
%   is used with $A\to B$, $B\to C$ and $Q_{1} \to Q_{2}$, etc.
%
%   Computing the product $Q_{1} \cdot Z$ is almost simple, since
%   $Q_{1}$ is rather small, but not quite:
%   the product of $Q_{1}$ with each block of four digits $Z_{i}$
%   is within \TeX{}'s bounds, but we wouldn't be able to use the
%   usual trick of adding a large power of $10$ to ensure that the
%   number of digits is fixed (see other operations for many examples
%   of this). Instead, we split off the digit of $10^5$ in $Q_{1}$
%   (and more generally $Q_{i}$), and do something similar to the
%   case of the full multiplication.
%
%   We know that $0<Q_{i}<1.8\cdot 10^5$, so $10^6+Q_{i}$ starts
%   with the digit $1$, followed by $\#1 = 1$ or $0$, then |#2|, which
%   is $5$ more digits. It would be somewhat simpler if we got |#1|
%   to be two digits, and |#2| four, but we are constrained by the $9$
%   arguments limit.
%
%   The result we want is then (the overall power of $10$ is arbitrary):
%   \begin{align*}
%   &10^{-4} ( \#3 - \#2 \cdot \#6 - 10 \cdot \#1 \cdot \#6\#7 )
%   + 10^{-8} ( \#4 - \#2 \cdot \#7 - 10 \cdot \#1 \cdot \#8 ) \\
%   &+ 10^{-12}( \#5 - \#2 \cdot \#8 - 10 \cdot \#1 \cdot \#9 )
%   + 10^{-16}(     - \#2 \cdot \#9 ).
%   \end{align*}
%   The factors of $10$ come from the fact that
%   $Q_{i} = 10\cdot 10^4 \cdot \#1 + \#2$. As usual, to combine
%   all the terms, we need to choose some shifts which must ensure
%   that the number of digits of the second, third, and fourth terms
%   are each fixed. Here, a good choice is $2\cdot 10^9$.
%   We are flirting with \TeX{}'s limits once more.
%
%   If $\#1=0$, then each term in parentheses (omitting the first)
%   is in the open interval $(-10^9, 10^4)$. Thus, adding
%   $2\cdot 10^9$ to it gives a $10$ digits number.\footnote{Bruno:
%     check that the carry from below does not screw that up. This
%   requires slightly tighter bounds.}
%
%   If $\#1=1$, then $\#2 < 7.8 \cdot 10^4$, and each term
%   in parentheses (omitting the first) is in the interval
%   $(-8\cdot 10^8, 10^4)$, and we are even safer.
%
%   We add the terms containing $\#1$ in a slightly tricky way
%   for efficiency reasons: if $\#1=0$, no need to do any computation,
%   while if $\#1=1$ we want $10$ times some number, simply obtained
%   by appending a $0$ digit.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_calc:Nwwnnnnnn 1#1#2; #3;#4#5 #6#7#8#9
  {
    1 #1 #2 \exp_after:wN ;
    \int_use:N \__int_eval:w
      - 200000 + #3 - #2 * #6
      \if_meaning:w 1 #1
        - #6#70
      \fi:
      +
      \exp_after:wN \@@_div_mantissa_calc_last:NNNNNN
      \int_use:N \__int_eval:w
        1999800000 + #4 - #2*#7
        \if_meaning:w 1 #1
          - #80
        \fi:
        +
        \exp_after:wN \@@_div_mantissa_calc_pack:NNNNNNw
        \int_use:N \__int_eval:w
          1999800000 + #5 - #2*#8
          \if_meaning:w 1 #1
            - #90
          \fi:
          +
          \exp_after:wN \@@_div_mantissa_calc_pack:NNNNNNw
          \int_use:N \__int_eval:w 2000000000 - #2*#9 ;
    {#6}{#7}{#8}{#9}
  }
\cs_new:Npn \@@_div_mantissa_calc_pack:NNNNNNw #1#2#3#4#5#6 #7;
  { #1#2#3#4#5#6 ; {#7} }
\cs_new:Npn \@@_div_mantissa_calc_last:NNNNNN #1#2#3#4#5#6
  { #1#2#3#4#5#6 \__int_eval_end: }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_iii:www} \meta{y} |;| \meta{$10^6+{}$Q_1} |;|
%     ~~\meta{B_1} \meta{B_2} |;| \Arg{B_3} \Arg{B_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%   \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_iii:www #1; #2; #3;
  {
    \exp_after:wN \@@_div_mantissa_iii_after:w
    \int_use:N \__int_eval:w #2
      \exp_after:wN \@@_div_mantissa_iv:www
      \__int_value:w #1 \exp_after:wN ;
      \__int_value:w
        \exp_after:wN \@@_div_mantissa_calc:Nwwnnnnnn
        \int_use:N \__int_eval:w #30/#1 + 999999 ;
        #3 ;
  }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_iv:www} \meta{y} |;| \meta{$10^6+{}$Q_2} |;|
%     ~~\meta{C_1} \meta{C_2} |;| \Arg{C_3} \Arg{C_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_iv:www #1; #2; #3;
  {
    \exp_after:wN \@@_div_mantissa_pack:NNN
    \int_use:N \__int_eval:w #2
      \exp_after:wN \@@_div_mantissa_v:www
      \__int_value:w #1 \exp_after:wN ;
      \__int_value:w
        \exp_after:wN \@@_div_mantissa_calc:Nwwnnnnnn
        \int_use:N \__int_eval:w #30/#1 + 999999 ;
        #3 ;
  }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_v:www} \meta{y} |;| \meta{$10^6+{}$Q_3} |;|
%     ~~\meta{D_1} \meta{D_2} |;| \Arg{D_3} \Arg{D_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_v:www #1; #2; #3;
  {
    \exp_after:wN \@@_div_mantissa_pack:NNN
    \int_use:N \__int_eval:w #2
      \exp_after:wN \@@_div_mantissa_vi:wwnnnn
      \__int_value:w
        \exp_after:wN \@@_div_mantissa_calc:Nwwnnnnnn
        \int_use:N \__int_eval:w #30/#1 + 999999 ;
        #3 ;
  }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_vi:wwnnnn} \meta{$10^6+{}$Q_4} |;|
%     ~~\meta{E_1} \meta{E_2} |;| \Arg{E_3} \Arg{E_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%   We compute $P$ by rounding $2 E_{1} E_{2}/Z_{1}Z_{2}$.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_vi:wwnnnn #1; #2;#3#4 #5#6
  {
    \exp_after:wN \@@_div_mantissa_pack:NNN
    \int_use:N \__int_eval:w #10
      \exp_after:wN \@@_div_mantissa_vii:wwnnnnnn
      \int_use:N \__int_eval:w (\c_two*#2)/#5#6 ; % <- P
      #2;{#3}{#4}
      {#5}{#6}
  }
%    \end{macrocode}
%   Note that we used |#10| instead of |#2| which we had previously.
%   Two reasons: firstly, since we dropped $y$, the argument which
%   holds $Q_{i}$ has changed, and secondly, we will want the
%   fourth piece of the result to have $5$ digits, including the
%   \meta{rounding} digit, which we shall compute now from $P$.
%   \begin{quote}
%     \cs{@@_div_mantissa_vii:wwnnnnnn} \meta{P} |;|
%     ~~\meta{E_1} \meta{E_2} |;| \Arg{E_3} \Arg{E_4}
%     ~~\Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \meta{sign}
%   \end{quote}
%   Then compute $2E-PZ$. Once more, we need to be careful and show
%   that the calculation $\#1\cdot\#5\#6$ below does not cause an
%   overflow: naively, $P$ can be up to $35$, and $\#5\#6$ up to
%   $10^8$, but both cannot happen simultaneously. To show that things
%   are fine, we split in two (non-disjoint) cases.
%   \begin{itemize}
%   \item For small $P$, say, $P< 10$, the product obeys
%     $P\cdot\#5\#6 < 10^8 \cdot P < 10^9 $.
%   \item For large $P$, say, $P\geq 3$, the rounding error on $P$,
%     which is at most $1$, is less than a factor of $2$, hence
%     $P\leq 4E/Z$, and $P\cdot \#5\#6 \leq 4E\cdot 10^8 < 10^9$.
%   \end{itemize}
%   Both inequalities could be made tighter if needed.
%
%   Note however that $P\cdot \#7\#8$ may overflow,
%   since the two factors are now independent, and the result may reach
%   $3.5\cdot 10^9$.
%
%   Also, we add $10\cdot P/2$ to the \enquote{fourth piece} of the result
%   as a first estimate of $10$ times $E/Z$. The goal is that the last digit
%   (for now $0$ or $5$) should be the \meta{rounding} digit. More precisely,
%   it will be corrected later by adding or subtracting $1$ depending on
%   whether $F$ was the correct integer part, or an overestimate (and nothing
%   is added when the quotient was exact). This does not give the
%   \enquote{correct} \meta{rounding} digit, but it always gives a digit
%   in the right \enquote{class} ($0$, $[1,4]$, $5$, or $[6-9]$), enough
%   for rounding purposes.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_vii:wwnnnnnn #1; #2;#3#4 #5#6#7#8
  {
    + \c_five * #1
    \exp_after:wN \@@_div_mantissa_ix:Nww
    \int_use:N \__int_eval:w -20 + 2*#2 - #1*#5#6 +
      \exp_after:wN \@@_div_mantissa_viii:NNw
      \int_use:N \__int_eval:w 199980 + 2*#3 - #1*#7 +
        \exp_after:wN \@@_div_mantissa_viii:NNw
        \int_use:N \__int_eval:w 200000 + 2*#4 - #1*#8 ; ;
  }
\cs_new:Npn \@@_div_mantissa_viii:NNw #1#2#3; { #1#2 ; + #3 }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_ix:Nww}
%     ~~\meta{F_1} \meta{F_2} |;| |+| \meta{F_3} |+| \meta{F_4} |;| \meta{sign}
%   \end{quote}
%   where $F=2E-PZ$. We only need to know whether it is positive,
%   negative, or exactly zero.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_ix:Nww #1#2;#3;
  {
    \if_meaning:w 0 #1
      \exp_after:wN \@@_div_mantissa_x:w
      \int_use:N \__int_eval:w #3
    \else:
      \if_meaning:w - #1
        -
      \else:
        +
      \fi:
      \c_one
    \fi:
    ;
  }
\cs_new:Npn \@@_div_mantissa_x:w #1;
  {
    \if_int_compare:w #1 > \c_zero
      + \c_one
    \fi:
    ;
  }
%    \end{macrocode}
%   We now obtain the following code, where \TeX{} is in the process
%   of expanding each of the integer expressions, and thus expands
%   the function at the bottom before the ones above it.
%   \begin{quote}
%     \cs{@@_div_mantissa_iii_after:w}  $10^6 + Q_{1}$
%     \cs{@@_div_mantissa_pack:NNN} $10^6 + Q_{2}$
%     \cs{@@_div_mantissa_pack:NNN} $10^6 + Q_{3}$
%     \cs{@@_div_mantissa_pack:NNN}
%     $10^7 + 10\cdot Q_{4} + 5 \cdot P + \varepsilon$ |;| \meta{sign}
%   \end{quote}
%   Here, $\varepsilon$ is $0$ in case $2E=PZ$ (\emph{i.e.}, $F=0$),
%   $1$ in case $2E>PZ$, which means that $P$ was the correct value,
%   but not with an exact quotient, and $-1$ if $2E<PZ$, \emph{i.e.},
%   $P$ was an overestimate.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_pack:NNN 1 #1 #2 { + #1 #2 ; }
%    \end{macrocode}
%   Once those have been expanded, we get
%   \begin{quote}
%     \cs{@@_div_mantissa_iii_after:w} |1| |0| \meta{5d} |;|
%     ~~\meta{4d} |;| \meta{4d} |;| \meta{5d} |;| \meta{sign}
%   \end{quote}
%   The reason we know that the first two digits are |1| and |0|
%   is that the final result is known to be between $0.1$ (inclusive)
%   and $10$, hence $\widetilde{Q_{1}}$ (the tilde denoting the
%   contribution from the other $Q_{i}$) is at most $99999$,
%   and $10^6+\widetilde{Q_{1}} = 10\cdots$.
%
%   It is now time to round. This depends on how many digits
%   the final result will have.
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_iii_after:w 10 #1
  {
    \if_meaning:w 0 #1
      \exp_after:wN \@@_div_mantissa_small:wwwNNNNwN
    \else:
      \exp_after:wN \@@_div_mantissa_large:wwwNNNNwN
    \fi:
    #1
  }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_small:wwwNNNNwN} |0| \meta{4d} |;|
%     ~~\meta{4d} |;| \meta{4d} |;| \meta{5d} |;| \meta{sign}
%   \end{quote}
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_small:wwwNNNNwN
    0 #1; #2; #3; #4#5#6#7#8; #9
  {
    \exp_after:wN \@@_basics_pack_high:NNNNNw
    \int_use:N \__int_eval:w 1 #1#2
      \exp_after:wN \@@_basics_pack_low:NNNNNw
      \int_use:N \__int_eval:w 1 #3#4#5#6#7
        + \@@_round:NNN #9 #7 #8
        \exp_after:wN ;
  }
%    \end{macrocode}
%   \begin{quote}
%     \cs{@@_div_mantissa_large:wwwNNNNwN} \meta{5d} |;|
%     ~~\meta{4d} |;| \meta{4d} |;| \meta{5d} |;| \meta{sign}
%   \end{quote}
%   \footnote{Bruno: rename the \enquote{add mantissa carry pack} function.}
%    \begin{macrocode}
\cs_new:Npn \@@_div_mantissa_large:wwwNNNNwN
    #1; #2; #3; #4#5#6#7#8; #9
  {
    + \c_one
    \exp_after:wN \@@_div_mantissa_large_pack:NNNNNNNNw
    \int_use:N \__int_eval:w 1 #1 #2 %<- 1+9d
      \exp_after:wN \@@_add_mantissa_carry_pack_ii:NNNNw
      \int_use:N \__int_eval:w 1 #3 #4 #5 #6
        + \@@_round:NNNN #9 #6 #7 #8
        \exp_after:wN ;
  }
\cs_new:Npn \@@_div_mantissa_large_pack:NNNNNNNNw
    1 #1#2#3#4 #5#6#7#8 #9; { ; {#1#2#3#4} {#5#6#7#8} {#9} }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Unary operations}
%
% \begin{macro}{\@@_neg:w}
%   This function flips the sign of the \meta{floating point}
%   and expands after it in the input stream, just like
%   \cs{@@_+_o:ww} etc.
%    \begin{macrocode}
\cs_new:Npn \@@_neg:w \s_@@ \@@_chk:w #1 #2
  {
    \exp_after:wN \@@_exp_after_o:w
    \exp_after:wN \s_@@
    \exp_after:wN \@@_chk:w
    \exp_after:wN #1
    \int_use:N \__int_eval:w \c_two - #2 \__int_eval_end:
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_abs:w}
%   This function sets the sign of the \meta{floating point} to be
%   positive, and expands after itself in the input stream, just like
%   \cs{@@_neg:w}.
%    \begin{macrocode}
\cs_new:Npn \@@_abs:w \s_@@ \@@_chk:w #1 #2
  { \@@_exp_after_o:w \s_@@ \@@_chk:w #1 0 }
%    \end{macrocode}
% \end{macro}
%
%    \begin{macrocode}
%</initex|package>
%    \end{macrocode}
%
% \end{implementation}
%
% \PrintChanges
%
% \PrintIndex