% \iffalse meta-comment
%
%% File: l3fp-extended.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{l3bootstrap}
\GetIdInfo$Id: l3fp-extended.dtx 4339 2012-11-24 19:16:43Z joseph $
  {L3 Floating-point extended precision fixed-points}
\documentclass[full]{l3doc}
\begin{document}
  \DocInput{\jobname.dtx}
\end{document}
%</driver>
% \fi
%
% \title{The \textsf{l3fp-extended} package\thanks{This file
%         has version number \ExplFileVersion, last
%         revised \ExplFileDate.}\\
% Fixed points with extended precision for internal use}
% \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-extended} implementation}
%
%    \begin{macrocode}
%<*initex|package>
%    \end{macrocode}
%
%    \begin{macrocode}
%<@@=fp>
%    \end{macrocode}
%
% \subsection{Description of extended fixed points}
%
% In this module, we work on (almost) fixed-point numbers with
% extended ($24$ digits) precision.  This is 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 trailing 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:Nw}.  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 extended fixed points}
%
% \begin{variable}[int]{\c_@@_one_fixed_tl}
%   The extended 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}[int, EXP]{\@@_fixed_continue:wn}
%   This function does nothing.  Of course, there is no bound on
%   $a_1$ (except \TeX{}'s own $2^{31}-1$).
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_continue:wn #1; #2 { #2 #1; }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[int, 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 \leq 2^{31} - 10001$.
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_add_one:wN #1#2; #3
  {
    \exp_after:wN #3 \exp_after:wN
      { \int_use:N \__int_eval:w \c_ten_thousand + #1 } #2 ;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, EXP]{\@@_fixed_mul_after:wn}
%   The fixed point operations which involve multiplication end by
%   calling this auxiliary.  It braces the last block of digits, and
%   places the \meta{continuation} |#2| in front.  The
%   \meta{continuation} was brought up through the expansions by
%   the packing functions.
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_mul_after:wn #1; #2 { #2 {#1} }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Dividing a fixed point number by a small integer}
%
% \begin{macro}[int, EXP]{\@@_fixed_div_int:wwN}
% \begin{macro}[aux, 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 will have $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_use:N \__int_eval:w \c_minus_one
      \@@_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_use:N \__int_eval:w #1 / #2 - \c_one ;
    {#2}
    {#1}
  }
\cs_new:Npn \@@_fixed_div_int_auxi:wnn #1; #2 #3
  {
    + #1
    \exp_after:wN \@@_fixed_div_int_pack:Nw
    \int_use:N \__int_eval:w 9999
      \exp_after:wN \@@_fixed_div_int:wnN
      \int_use:N \__int_eval:w #3 - #1*#2 \__int_eval_end:
  }
\cs_new:Npn \@@_fixed_div_int_auxii:wnn #1; #2 #3 { + #1 + \c_two ; }
\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}[int, EXP]{\@@_fixed_add:wwn, \@@_fixed_sub:wwn}
% \begin{macro}[aux, 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}<50000$, and requires the result to be positive (this
%   happens automatically for addition).  The two functions only differ
%   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 two signs, $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_nopar:Npn \@@_fixed_add:wwn { \@@_fixed_add:Nnnnnwnn + }
\cs_new_nopar: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_use:N \__int_eval:w 9 9999 9998 + #2#3 #1 #7#8
      \exp_after:wN \@@_fixed_add_pack:NNNNNwn
      \int_use:N \__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_use:N \__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}
%
% \begin{macro}[int, EXP]{\@@_fixed_mul:wwn}
% \begin{macro}[aux, EXP]{\@@_fixed_mul:nnnnnnnwn}
%   \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} \\
%     & + \left(a_{2} \cdot b_{4} + a_{3} \cdot b_{3} + a_{4} \cdot b_{2}
%           + \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}}
%           + a_{1} \cdot b_{5} + a_{5} \cdot b_{1}\right) \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 packing auxiliaries bring the
%   \meta{continuation} up through the expansion chain, as |#7|, and it
%   is finally placed in front of the $6$ brace groups by
%   \cs{@@_fixed_mul_after:wn}.
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_mul:wwn #1#2#3#4 #5; #6#7#8#9
  {
    \exp_after:wN \@@_fixed_mul_after:wn
    \int_use:N \__int_eval:w \c_@@_leading_shift_int
      \exp_after:wN \@@_pack:NNNNNwn
      \int_use:N \__int_eval:w \c_@@_middle_shift_int
        + #1*#6
        \exp_after:wN \@@_pack:NNNNNwn
        \int_use:N \__int_eval:w \c_@@_middle_shift_int
          + #1*#7 + #2*#6
          \exp_after:wN \@@_pack:NNNNNwn
          \int_use:N \__int_eval:w \c_@@_middle_shift_int
            + #1*#8 + #2*#7 + #3*#6
            \exp_after:wN \@@_pack:NNNNNwn
            \int_use:N \__int_eval:w \c_@@_middle_shift_int
              + #1*#9 + #2*#8 + #3*#7 + #4*#6
              \exp_after:wN \@@_pack:NNNNNwn
              \int_use:N \__int_eval:w \c_@@_trailing_shift_int
                + #2*#9 + #3*#8 + #4*#7
                + ( #3*#9 + #4*#8
                  + \@@_fixed_mul:nnnnnnnwn #5 {#6}{#7}  {#1}{#2}
  }
\cs_new:Npn \@@_fixed_mul:nnnnnnnwn #1#2 #3#4 #5#6 #7#8 ; #9
  {
    #1*#4 + #2*#3 + #5*#8 + #6*#7 ) / \c_ten_thousand
    + #1*#3 + #5*#7 ;
    {#9} ;
  }
%    \end{macrocode}
% \end{macro}
% \end{macro}
%
% \subsection{Combining product and sum of fixed points}
%
% \begin{macro}[int, 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}
%   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 will 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}
%           + \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}}
%           + 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}, will be taken in the
%   integer expression for the $10^{-24}$ level.  The
%   \meta{continuation} is placed correctly to be taken upstream by
%   packing auxiliaries.
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_mul_add:wwwn #1; #2; #3#4#5#6#7#8; #9
  {
    \exp_after:wN \@@_fixed_mul_after:wn
    \int_use:N \__int_eval:w \c_@@_big_leading_shift_int
      \exp_after:wN \@@_pack_big:NNNNNNwn
      \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #3 #4
        \@@_fixed_mul_add:Nwnnnwnnn +
          + #5 #6 ; #2 ; #1 ; #2 ; +
          + #7 #8 ; {#9} ;
  }
\cs_new:Npn \@@_fixed_mul_sub_back:wwwn #1; #2; #3#4#5#6#7#8; #9
  {
    \exp_after:wN \@@_fixed_mul_after:wn
    \int_use:N \__int_eval:w \c_@@_big_leading_shift_int
      \exp_after:wN \@@_pack_big:NNNNNNwn
      \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #3 #4
        \@@_fixed_mul_add:Nwnnnwnnn -
          + #5 #6 ; #2 ; #1 ; #2 ; -
          + #7 #8 ; {#9} ;
  }
\cs_new:Npn \@@_fixed_one_minus_mul:wwn #1; #2; #3
  {
    \exp_after:wN \@@_fixed_mul_after:wn
    \int_use:N \__int_eval:w \c_@@_big_leading_shift_int
      \exp_after:wN \@@_pack_big:NNNNNNwn
      \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + 1 0000 0000
        \@@_fixed_mul_add:Nwnnnwnnn -
          ; #2 ; #1 ; #2 ; -
          ; {#3} ;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, 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 huse 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:NNNNNNwn
    \int_use:N \__int_eval:w \c_@@_big_middle_shift_int
      #1 #7*#4 #1 #8*#3
      \exp_after:wN \@@_pack_big:NNNNNNwn
      \int_use:N \__int_eval:w \c_@@_big_middle_shift_int
        #1 #7*#5 #1 #8*#4 #1 #9*#3 #2
        \exp_after:wN \@@_pack_big:NNNNNNwn
        \int_use:N \__int_eval:w \c_@@_big_middle_shift_int
          #1 \@@_fixed_mul_add:nnnnwnnnn {#7}{#8}{#9}
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, 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 will
%   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:NNNNNNwn
    \int_use:N \__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}[aux, 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_ten_thousand
  }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Converting from fixed point to floating point}
%
% \begin{macro}[int, rEXP]
%   {\@@_fixed_to_float:wN, \@@_fixed_to_float:Nw}
%   \begin{syntax}
%     \ldots{} \cs{__int_eval:w} \meta{exponent} \cs{@@_fixed_to_float: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 to know the sign.
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_to_float:Nw #1#2; { \@@_fixed_to_float:wN #2; #1 }
\cs_new:Npn \@@_fixed_to_float:wN #1#2#3#4#5#6; #7
  {
    + \c_four % for the 8-digit-at-the-start thing.
    \exp_after:wN \exp_after:wN
    \exp_after:wN \@@_fixed_to_loop:N
    \exp_after:wN \use_none:n
    \int_use:N \__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
      - \c_one
      \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 ;
  {
    - \c_two * \c_@@_max_exponent_int ;
    {0000} {0000} {0000} {0000} ;
  }
\cs_new:Npn \@@_fixed_to_float_pack:ww #1 ; #2#3 ; ;
  {
    \if_int_compare:w #2 > \c_four
      \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_use:N \__int_eval:w 1 #1#2
      \exp_after:wN \@@_basics_pack_low:NNNNNw
      \int_use:N \__int_eval:w 1 #3#4 + \c_one ;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[rEXP, int]{\@@_fixed_inv_to_float:wN, \@@_fixed_div_to_float:ww}
%   Starting from \texttt{fixed_dtf} $A$ |;| $B$ |;| we want to compute
%   $A/B$, and express it as a floating point number.  Normalize both
%   numbers by removing leading brace groups of zeros and leaving the
%   appropriate exponent shift in the input stream.
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_inv_to_float:wN #1#2; #3
  {
    + \__int_eval:w % ^^A todo: remove the +?
        \if_int_compare:w #1 < \c_one_thousand
          \@@_fixed_dtf_zeros:wNnnnnnn
        \fi:
        \@@_fixed_dtf_no_zero:Nwn + {#1} #2 \s_@@
        \@@_fixed_dtf_approx:n
        {10000} {0000} {0000} {0000} {0000} {0000} ;
  }
\cs_new:Npn \@@_fixed_div_to_float:ww #1#2; #3#4;
  {
    \if_int_compare:w #1 < \c_one_thousand
      \@@_fixed_dtf_zeros:wNnnnnnn
    \fi:
    \@@_fixed_dtf_no_zero:Nwn - {#1} #2 \s_@@
    {
      \if_int_compare:w #3 < \c_one_thousand
        \@@_fixed_dtf_zeros:wNnnnnnn
      \fi:
      \@@_fixed_dtf_no_zero:Nwn + {#3} #4 \s_@@
      \@@_fixed_dtf_approx:n
    }
  }
\cs_new:Npn \@@_fixed_dtf_no_zero:Nwn #1#2 \s_@@ #3 { #3 #2; }
\cs_new:Npn \@@_fixed_dtf_zeros:wNnnnnnn
    \fi: \@@_fixed_dtf_no_zero:Nwn #1#2#3#4#5#6#7
  {
    \fi:
    #1 \c_minus_one
    \exp_after:wN \use_i_ii:nnn
    \exp_after:wN \@@_fixed_dtf_zeros:NN
    \exp_after:wN #1
    \int_use:N \__int_eval:w 10 0000 + #2 \__int_eval_end: #3#4#5#6#7
    ; 1 ;
  }
\cs_new:Npn \@@_fixed_dtf_zeros:NN #1#2
  {
    \if_meaning:w 0 #2
      #1 \c_one
    \else:
      \@@_fixed_dtf_zeros_end:wNww #2
    \fi:
    \@@_fixed_dtf_zeros:NN #1
  }
\cs_new:Npn \@@_fixed_dtf_zeros_end:wNww
    #1 \fi: \@@_fixed_dtf_zeros:NN #2 #3; #4 \s_@@
  {
    \fi:
    \if_meaning:w ; #1
      #2 \c_two * \c_@@_max_exponent_int
      \use_i_ii:nnn
    \fi:
    \@@_fixed_dtf_zeros_auxi:ww
    #1#3 0000 0000 0000 0000 0000 0000 ;
  }
\cs_new:Npn \@@_fixed_dtf_zeros_auxi:ww
  {
    \@@_pack_twice_four:wNNNNNNNN
    \@@_pack_twice_four:wNNNNNNNN
    \@@_pack_twice_four:wNNNNNNNN
    \@@_fixed_dtf_zeros_auxii:ww
    ;
  }
\cs_new:Npn \@@_fixed_dtf_zeros_auxii:ww #1; #2; #3 { #3 #1; }
%    \end{macrocode}
%   \newcommand{\eTeXfrac}[2]{\left[\frac{#1}{#2}\right]}
%   We get
%   \begin{quote}
%     \cs{@@_fixed_dtf_approx:n} \meta{B'} |;| \meta{A'} |;|
%   \end{quote}
%   where \meta{B'} and \meta{A'} are each $6$ brace groups,
%   representing fixed point numbers in the range $[0.1,1)$.  Denote by
%   $x\in[1000,9999]$ and $y\in[0,9999]$ the first two groups of
%   \meta{B'}.  We first find an estimate $a$ for the inverse of $B'$ by
%   computing
%   \begin{align*}
%     \alpha &= \eTeXfrac{10^{9}}{x+1} \\
%     \beta  &= \eTeXfrac{10^{9}}{x} \\
%     a &= 10^{3} \alpha + (\beta-\alpha) \cdot
%       \left(10^{3}-\eTeXfrac{y}{10}\right) - 1750,
%   \end{align*}
%   where $\eTeXfrac{\bullet}{\bullet}$ denotes \eTeX{}'s rounding
%   division.  The idea is to interpolate between $\alpha$ and $\beta$
%   with a parameter $y/10^{4}$.  The shift by $1750$ helps to ensure
%   that $a$ is an underestimate of the correct value.  We will prove
%   that
%   \[
%   1 - 2.255\cdot 10^{-5} < \frac{B'a}{10^{8}} < 1 .
%   \]
%   We can then compute the inverse $B'a/10^{8}$ using $1/(1-\epsilon)
%   \simeq (1+\epsilon)(1+\epsilon^{2})$, which is correct up to a
%   relative error of $\epsilon^4 < 2.6\cdot 10^{-19}$.  Since we target
%   a $16$-digit value, this is small enough.
%
%   Let us prove the upper bound first.
%   \begin{align}\label{l3fp-fixed-eTeXfrac}
%     10^{7} B'a
%       & < \left(10^{3} x + \eTeXfrac{y}{10} + \frac{3}{2}\right)
%       \left(\left(10^{3}-\eTeXfrac{y}{10}\right) \beta
%         + \eTeXfrac{y}{10} \alpha - 1750\right)
%     \\& < \left(10^{3} x + \eTeXfrac{y}{10} + \frac{3}{2}\right)
%       \left(\left(10^{3}-\eTeXfrac{y}{10}\right)
%         \left(\frac{10^{9}}{x} + \frac{1}{2} \right)
%         + \eTeXfrac{y}{10} \left(\frac{10^{9}}{x+1} + \frac{1}{2} \right)
%         - 1750\right)
%     \\& < \left(10^{3} x + \eTeXfrac{y}{10} + \frac{3}{2}\right)
%       \left(\frac{10^{12}}{x}
%         - \eTeXfrac{y}{10} \frac{10^{9}}{x(x+1)}
%         - 1250\right)
%   \end{align}
%   We recognize a quadratic polynomial in $[y/10]$ with a negative
%   leading coefficient, $([y/10]+a)(b-c[y/10]) \leq (b+ca)^2/(4c)$.
%   Hence,
%   \[
%   10^{7} B'a
%   < \frac{10^{15}}{x(x+1)} \left(
%     x + \frac{1}{2} + \frac{3}{4} 10^{-3}
%     - 6.25\cdot 10^{-10} x(x+1) \right)^2
%   \]
%   We want to prove that the squared expression is less than $x(x+1)$,
%   which we do by simplifying the difference, and checking its sign,
%   \[
%   x(x+1) - \left(x + \frac{1}{2} + \frac{3}{4} 10^{-3}
%     - 6.25\cdot 10^{-10} x(x+1) \right)^2
%   > - \frac{1}{4} (1+1.5\cdot 10^{-3})^2 - 10^{-3} x
%     + 1.25\cdot 10^{-9} x(x+1)(x+0.5)
%   > 0.
%   \]
%
%   Now, the lower bound.  The same computation as
%   \eqref{l3fp-fixed-eTeXfrac} imply
%   \[
%     10^{7} B'a
%       > \left(10^{3} x + \eTeXfrac{y}{10} - \frac{1}{2}\right)
%       \left(\frac{10^{12}}{x} - \eTeXfrac{y}{10} \frac{10^{9}}{x(x+1)}
%         - 2250\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=0$ or $y=9999$, and we
%   easily check the bound for those values.
%
%   We have proven that the algorithm will give us a precise enough
%   answer.  Incidentally, the upper bound that we derived tells us that
%   $a < 10^{8}/B \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.
%    \begin{macrocode}
\cs_new:Npn \@@_fixed_dtf_approx:n #1
  {
    \exp_after:wN \@@_fixed_dtf_approx:wnn
    \int_use:N \__int_eval:w 10 0000 0000 / ( #1 + \c_one ) ;
    {#1}
  }
\cs_new:Npn \@@_fixed_dtf_approx:wnn #1; #2#3
  {
%<assert>    \assert:n { \tl_count:n {#1} = 6 }
    \exp_after:wN \@@_fixed_dtf_approx:NNNNNw
    \int_use:N \__int_eval:w 10 0000 0000 - 1750
      + #1000 + (10 0000 0000/#2-#1) * (1000-#3/10) ;
    {#2}{#3}
  }
\cs_new:Npn \@@_fixed_dtf_approx:NNNNNw 1#1#2#3#4#5#6; #7; #8;
  {
    + \c_four % because of the line below "dtf_epsilon" here.
    \@@_fixed_mul:wwn {000#1}{#2#3#4#5}{#6}{0000}{0000}{0000} ; #7;
    \@@_fixed_dtf_epsilon:wN
    \@@_fixed_mul:wwn {000#1}{#2#3#4#5}{#6}{0000}{0000}{0000} ;
    \@@_fixed_mul:wwn #8;
    \@@_fixed_to_float:wN ?
  }
\cs_new:Npn \@@_fixed_dtf_epsilon:wN #1#2#3#4#5#6;
  {
%<assert>    \assert:n { #1 = 0000 }
%<assert>    \assert:n { #2 = 9999 }
    \exp_after:wN \@@_fixed_dtf_epsilon:NNNNNww
    \int_use:N \__int_eval:w 1 9999 9998 - #3#4 +
      \exp_after:wN \@@_fixed_dtf_epsilon_pack:NNNNNw
      \int_use:N \__int_eval:w 2 0000 0000 - #5#6 ; {0000} ;
  }
\cs_new:Npn \@@_fixed_dtf_epsilon_pack:NNNNNw #1#2#3#4#5#6;
  { #1 ; {#2#3#4#5} {#6} }
\cs_new:Npn \@@_fixed_dtf_epsilon:NNNNNww #1#2#3#4#5#6; #7;
  {
    \@@_fixed_mul:wwn %^^A todo: optimize to use \@@_mul_significand.
      {0000} {#2#3#4#5} {#6} #7 ;
      {0000} {#2#3#4#5} {#6} #7 ;
    \@@_fixed_add_one:wN
    \@@_fixed_mul:wwn {10000} {#2#3#4#5} {#6} #7 ;
  }
%    \end{macrocode}
% \end{macro}
%
%    \begin{macrocode}
%</initex|package>
%    \end{macrocode}
%
% \end{implementation}
%
% \PrintChanges
%
% \PrintIndex
