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

% \iffalse meta-comment
%
%% File: l3fp-expo.dtx Copyright (C) 2011-2014,2016,2017 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 development version of the bundle can be found at
%
%    https://github.com/latex3/latex3
%
% for those people who are interested.
%
%<*driver>
\documentclass[full]{l3doc}
\def\ExplFileDate{2017/03/18}
\def\ExplFileVersion{7019}
\begin{document}
  \DocInput{\jobname.dtx}
\end{document}
%</driver>
% \fi
%
% \title{The \textsf{l3fp-expo} package\thanks{This file
%         has version number \ExplFileVersion, last
%         revised \ExplFileDate.}\\
% Floating point exponential-related functions}
% \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-expo} implementation}
%
%    \begin{macrocode}
%<*initex|package>
%    \end{macrocode}
%
%    \begin{macrocode}
%<@@=fp>
%    \end{macrocode}
%
% \begin{macro}[aux, EXP]
%   {
%     \@@_parse_word_exp:N   ,
%     \@@_parse_word_ln:N    ,
%   }
%   Unary functions.
%    \begin{macrocode}
\cs_new:Npn \@@_parse_word_exp:N
  { \@@_parse_unary_function:NNN \@@_exp_o:w ? }
\cs_new:Npn \@@_parse_word_ln:N
  { \@@_parse_unary_function:NNN \@@_ln_o:w ? }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Logarithm}
%
% \subsubsection{Work plan}
%
% As for many other functions, we filter out special cases in
% \cs{@@_ln_o:w}.  Then \cs{@@_ln_npos_o:w} receives a positive normal
% number, which we write in the form $a\cdot 10^{b}$ with $a\in[0.1,1)$.
%
% \emph{The rest of this section is actually not in sync with the code.
%   Or is the code not in sync with the section?  In the current code,
%   $c\in [1,10]$ will be such that $0.7\leq ac < 1.4$.}
%
% We are given a positive normal number, of the form $a\cdot 10^{b}$
% with $a\in[0.1,1)$.  To compute its logarithm, we find a small integer
% $5\leq c < 50$ such that $0.91 \leq a c / 5 < 1.1$, and use the
% relation
% \begin{equation*}
%   \ln (a \cdot 10^b) = b \cdot \ln (10) - \ln (c/5) + \ln (ac/5).
% \end{equation*}
% The logarithms $\ln(10)$ and $\ln(c/5)$ are looked up in a table.  The
% last term is computed using the following Taylor series of $\ln$ near
% $1$:
% \begin{equation*}
%   \ln\left(\frac{ac}{5}\right)
%   = \ln\left(\frac{1+t}{1-t}\right)
%   = 2t\left(1 + t^2 \left(\frac{1}{3} + t^2 \left(\frac{1}{5}
%         + t^2 \left(\frac{1}{7} + t^2 \left( \frac{1}{9} + \cdots
%           \right)\right)\right)\right)\right)
% \end{equation*}
% where $t=1-10/(ac+5)$.  We can now see one reason for the choice of
% $ac\sim 5$: then $ac+5=10(1-\epsilon)$ with $-0.05<\epsilon\leq
% 0.045$, hence
% \begin{equation*}
%   t = \frac{\epsilon}{1-\epsilon}
%   = \epsilon (1+\epsilon)(1+\epsilon^2)(1+\epsilon^4)\ldots,
% \end{equation*}
% is not too difficult to compute.
%
% \subsubsection{Some constants}
%
% \begin{variable}[aux]
%   {
%     \c_@@_ln_i_fixed_tl ,
%     \c_@@_ln_ii_fixed_tl ,
%     \c_@@_ln_iii_fixed_tl ,
%     \c_@@_ln_iv_fixed_tl ,
%     \c_@@_ln_vi_fixed_tl ,
%     \c_@@_ln_vii_fixed_tl ,
%     \c_@@_ln_viii_fixed_tl ,
%     \c_@@_ln_ix_fixed_tl ,
%     \c_@@_ln_x_fixed_tl,
%   }
%   A few values of the logarithm as extended fixed point numbers.
%   Those are needed in the implementation.  It turns out that we don't
%   need the value of $\ln(5)$.
%    \begin{macrocode}
\tl_const:Nn \c_@@_ln_i_fixed_tl   { {0000}{0000}{0000}{0000}{0000}{0000};}
\tl_const:Nn \c_@@_ln_ii_fixed_tl  { {6931}{4718}{0559}{9453}{0941}{7232};}
\tl_const:Nn \c_@@_ln_iii_fixed_tl {{10986}{1228}{8668}{1096}{9139}{5245};}
\tl_const:Nn \c_@@_ln_iv_fixed_tl  {{13862}{9436}{1119}{8906}{1883}{4464};}
\tl_const:Nn \c_@@_ln_vi_fixed_tl  {{17917}{5946}{9228}{0550}{0081}{2477};}
\tl_const:Nn \c_@@_ln_vii_fixed_tl {{19459}{1014}{9055}{3133}{0510}{5353};}
\tl_const:Nn \c_@@_ln_viii_fixed_tl{{20794}{4154}{1679}{8359}{2825}{1696};}
\tl_const:Nn \c_@@_ln_ix_fixed_tl  {{21972}{2457}{7336}{2193}{8279}{0490};}
\tl_const:Nn \c_@@_ln_x_fixed_tl   {{23025}{8509}{2994}{0456}{8401}{7991};}
%    \end{macrocode}
% \end{variable}
%
% \subsubsection{Sign, exponent, and special numbers}
%
% \begin{macro}[EXP, int]{\@@_ln_o:w}
%   The logarithm of negative numbers (including $-\infty$ and $-0$)
%   raises the \enquote{invalid} exception.  The logarithm of $+0$ is
%   $-\infty$, raising a division by zero exception.  The logarithm of
%   $+\infty$ or a \texttt{nan} is itself.  Positive normal numbers call
%   \cs{@@_ln_npos_o:w}.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
  {
    \if_meaning:w 2 #3
      \@@_case_use:nw { \@@_invalid_operation_o:nw { ln } }
    \fi:
    \if_case:w #2 \exp_stop_f:
      \@@_case_use:nw
        { \@@_division_by_zero_o:Nnw \c_minus_inf_fp { ln } }
    \or:
    \else:
      \@@_case_return_same_o:w
    \fi:
    \@@_ln_npos_o:w \s_@@ \@@_chk:w #2#3#4;
  }
%    \end{macrocode}
% \end{macro}
%
% \subsubsection{Absolute ln}
%
% \begin{macro}[aux, EXP]{\@@_ln_npos_o:w}
%   We catch the case of a significand very close to $0.1$ or to $1$.
%   In all other cases, the final result is at least $10^{-4}$, and
%   then an error of $0.5\cdot 10^{-20}$ is acceptable.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_npos_o:w \s_@@ \@@_chk:w 10#1#2#3;
  { %^^A todo: ln(1) should be "exact zero", not "underflow"
    \exp_after:wN \@@_sanitize:Nw
    \__int_value:w % for the overall sign
      \if_int_compare:w #1 < 1 \exp_stop_f:
        2
      \else:
        0
      \fi:
      \exp_after:wN \exp_stop_f:
      \__int_value:w \__int_eval:w % for the exponent
        \@@_ln_significand:NNNNnnnN #2#3
        \@@_ln_exponent:wn {#1}
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, EXP]{\@@_ln_significand:NNNNnnnN}
%   \begin{quote}
%     \cs{@@_ln_significand:NNNNnnnN} \meta{X_1} \Arg{X_2} \Arg{X_3}
%     \Arg{X_4} \meta{continuation}
%   \end{quote}
%   This function expands to
%   \begin{quote}
%     \meta{continuation} \Arg{Y_1} \Arg{Y_2} \Arg{Y_3} \Arg{Y_4}
%     \Arg{Y_5} \Arg{Y_6} |;|
%   \end{quote}
%   where $Y = - \ln(X)$ as an extended fixed point.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_significand:NNNNnnnN #1#2#3#4
  {
    \exp_after:wN \@@_ln_x_ii:wnnnn
    \__int_value:w
      \if_case:w #1 \exp_stop_f:
      \or:
        \if_int_compare:w #2 < 4 \exp_stop_f:
          \__int_eval:w 10 - #2
        \else:
          6
        \fi:
      \or: 4
      \or: 3
      \or: 2
      \or: 2
      \or: 2
      \else: 1
      \fi:
    ; { #1 #2 #3 #4 }
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, EXP]{\@@_ln_x_ii:wnnnn}
%   We have thus found $c \in [1,10]$ such that $0.7\leq ac < 1.4$
%   in all cases. Compute $ 1 + x = 1 + ac \in [1.7,2.4)$.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_x_ii:wnnnn #1; #2#3#4#5
  {
    \exp_after:wN \@@_ln_div_after:Nw
    \cs:w c_@@_ln_ \__int_to_roman:w #1 _fixed_tl \exp_after:wN \cs_end:
    \__int_value:w
      \exp_after:wN \@@_ln_x_iv:wnnnnnnnn
      \__int_value:w \__int_eval:w
        \exp_after:wN \@@_ln_x_iii_var:NNNNNw
        \__int_value:w \__int_eval:w 9999 9990 + #1*#2#3 +
          \exp_after:wN \@@_ln_x_iii:NNNNNNw
          \__int_value:w \__int_eval:w 10 0000 0000 + #1*#4#5 ;
    {20000} {0000} {0000} {0000}
  } %^^A todo: reoptimize (a generalization attempt failed).
\cs_new:Npn \@@_ln_x_iii:NNNNNNw #1#2 #3#4#5#6 #7;
  { #1#2; {#3#4#5#6} {#7} }
\cs_new:Npn \@@_ln_x_iii_var:NNNNNw #1 #2#3#4#5 #6;
  {
    #1#2#3#4#5 + 1 ;
    {#1#2#3#4#5} {#6}
  }
%    \end{macrocode}
%   The Taylor series will be expressed in terms of
%   $t = (x-1)/(x+1) = 1 - 2/(x+1)$. We now compute the
%   quotient with extended precision, reusing some code
%   from \cs{@@_/_o:ww}. Note that $1+x$ is known exactly.
%
%   To reuse notations from \pkg{l3fp-basics}, we want to
%   compute $ A / Z $ with $ A = 2 $ and $ Z = x + 1 $.
%   In \pkg{l3fp-basics}, we considered the case where
%   both $A$ and $Z$ are arbitrary, in the range $[0.1,1)$,
%   and we had to monitor the growth of the sequence of
%   remainders $A$, $B$, $C$, etc. to ensure that no overflow
%   occurred during the computation of the next quotient.
%   The main source of risk was our choice to define the
%   quotient as roughly $10^9 \cdot A / 10^5 \cdot Z$: then
%   $A$ was bound to be below $2.147\cdots$, and this limit
%   was never far.
%
%   In our case, we can simply work with $10^8 \cdot A$ and
%   $10^4 \cdot Z$, because our reason to work with higher
%   powers has gone: we needed the integer $y \simeq 10^5 \cdot Z$
%   to be at least $10^4$, and now, the definition
%   $y \simeq 10^4 \cdot Z$ suffices.
%
%   Let us thus define $y = \left\lfloor 10^4 \cdot Z \right\rfloor + 1
%   \in ( 1.7 \cdot 10^4, 2.4 \cdot 10^4 ] $, and
%   \[
%   Q_{1}
%   =
%   \left\lfloor
%     \frac {\left\lfloor 10^8 \cdot A\right\rfloor} {y} - \frac{1}{2}
%   \right\rfloor.
%   \]
%   (The $1/2$ comes from how e\TeX{} rounds.) As for division, it is
%   easy to see that $Q_{1} \leq 10^4 A / Z$, \emph{i.e.}, $Q_{1}$
%   is an underestimate.
%
%   Exactly as we did for division, we set $B = 10^4 A - Q_{1}Z$. Then
%   \begin{align*}
%     10^4 B
%     & \leq
%     A_{1}A_{2}.A_{3}A_{4}
%     - \left( \frac{A_{1}A_{2}}{y} - \frac{3}{2} \right) 10^4 Z
%     \\
%     & \leq
%     A_{1}A_{2} \left( 1 - \frac{10^4 Z}{y} \right) + 1 + \frac{3}{2} y
%     \\
%     & \leq
%     10^8 \frac{A}{y} + 1 + \frac{3}{2} y
%   \end{align*}
%   In the same way, and using $1.7\cdot 10^4\leq y\leq 2.4\cdot 10^4$,
%   and convexity, we get
%   \begin{align*}
%     10^4 A &= 2\cdot 10^4 \\
%     10^4 B &\leq 10^8 \frac{A}{y} + 1.6 y \leq 4.7\cdot 10^4\\
%     10^4 C &\leq 10^8 \frac{B}{y} + 1.6 y \leq 5.8\cdot 10^4\\
%     10^4 D &\leq 10^8 \frac{C}{y} + 1.6 y \leq 6.3\cdot 10^4\\
%     10^4 E &\leq 10^8 \frac{D}{y} + 1.6 y \leq 6.5\cdot 10^4\\
%     10^4 F &\leq 10^8 \frac{E}{y} + 1.6 y \leq 6.6\cdot 10^4\\
%   \end{align*}
%   Note that we compute more steps than for division: since $t$ is
%   not the end result, we need to know it with more accuracy
%   (on the other hand, the ending is much simpler, as we don't
%   need an exact rounding for transcendental functions, but just
%   a faithful rounding).
%   ^^A todo: doc
%
%   \begin{quote}
%     \cs{@@_ln_x_iv:wnnnnnnnn}
%     \meta{1 or 2} \meta{8d} |;| \Arg{4d} \Arg{4d} \meta{fixed-tl}
%   \end{quote}
%   The number is $x$. Compute $y$ by adding 1 to the five first digits.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_x_iv:wnnnnnnnn #1; #2#3#4#5 #6#7#8#9
  {
    \exp_after:wN \@@_div_significand_pack:NNN
    \__int_value:w \__int_eval:w
    \@@_ln_div_i:w #1 ;
      #6 #7 ; {#8} {#9}
      {#2} {#3} {#4} {#5}
      { \exp_after:wN \@@_ln_div_ii:wwn \__int_value:w #1 }
      { \exp_after:wN \@@_ln_div_ii:wwn \__int_value:w #1 }
      { \exp_after:wN \@@_ln_div_ii:wwn \__int_value:w #1 }
      { \exp_after:wN \@@_ln_div_ii:wwn \__int_value:w #1 }
      { \exp_after:wN \@@_ln_div_vi:wwn \__int_value:w #1 }
  }
\cs_new:Npn \@@_ln_div_i:w #1;
  {
    \exp_after:wN \@@_div_significand_calc:wwnnnnnnn
    \__int_value:w \__int_eval:w 999999 + 2 0000 0000 / #1 ; % Q1
  }
\cs_new:Npn \@@_ln_div_ii:wwn #1; #2;#3 % y; B1;B2 <- for k=1
  {
    \exp_after:wN \@@_div_significand_pack:NNN
    \__int_value:w \__int_eval:w
      \exp_after:wN \@@_div_significand_calc:wwnnnnnnn
      \__int_value:w \__int_eval:w 999999 + #2 #3 / #1 ; % Q2
      #2 #3 ;
  }
\cs_new:Npn \@@_ln_div_vi:wwn #1; #2;#3#4#5 #6#7#8#9 %y;F1;F2F3F4x1x2x3x4
  {
    \exp_after:wN \@@_div_significand_pack:NNN
    \__int_value:w \__int_eval:w 1000000 + #2 #3 / #1 ; % Q6
  }
%    \end{macrocode}
%   We now have essentially
%   ^^A todo: determine error on $Q_{6}$ (probably $6.7$),
%   ^^A todo: conclude the final result is off by $<10^{-23}$
%   \begin{quote}
%     \cs{@@_ln_div_after:Nw} \meta{fixed tl}
%     \cs{@@_div_significand_pack:NNN} $10^6 + Q_{1}$
%     \cs{@@_div_significand_pack:NNN} $10^6 + Q_{2}$
%     \cs{@@_div_significand_pack:NNN} $10^6 + Q_{3}$
%     \cs{@@_div_significand_pack:NNN} $10^6 + Q_{4}$
%     \cs{@@_div_significand_pack:NNN} $10^6 + Q_{5}$
%     \cs{@@_div_significand_pack:NNN} $10^6 + Q_{6}$ |;|
%     \meta{exponent} |;| \meta{continuation}
%   \end{quote}
%   where \meta{fixed tl} holds the logarithm of a number
%   in $[1,10]$, and \meta{exponent} is
%   the exponent. Also, the expansion is done backwards. Then
%   \cs{@@_div_significand_pack:NNN} puts things in the
%   correct order to add the $Q_{i}$ together and put semicolons
%   between each piece. Once those have been expanded, we get
%   \begin{quote}
%     \cs{@@_ln_div_after:Nw} \meta{fixed-tl} \meta{1d} |;| \meta{4d} |;| \meta{4d} |;| \meta{4d} |;| \meta{4d} |;| \meta{4d} |;| \meta{4d} |;| \meta{exponent} |;|
%   \end{quote}
%   ^^A todo: redoc.
%   Just as with division, we know that the first two digits
%   are |1| and |0| because of bounds on the final result of
%   the division $2/(x+1)$, which is between roughly $0.8$ and $1.2$.
%   We then compute $1-2/(x+1)$, after testing whether $2/(x+1)$ is
%   greater than or smaller than $1$.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_div_after:Nw #1#2;
  {
    \if_meaning:w 0 #2
      \exp_after:wN \@@_ln_t_small:Nw
    \else:
      \exp_after:wN \@@_ln_t_large:NNw
      \exp_after:wN -
    \fi:
    #1
  }
\cs_new:Npn \@@_ln_t_small:Nw #1 #2; #3; #4; #5; #6; #7;
  {
    \exp_after:wN \@@_ln_t_large:NNw
    \exp_after:wN + % <sign>
    \exp_after:wN #1
    \__int_value:w \__int_eval:w 9999 - #2 \exp_after:wN ;
    \__int_value:w \__int_eval:w 9999 - #3 \exp_after:wN ;
    \__int_value:w \__int_eval:w 9999 - #4 \exp_after:wN ;
    \__int_value:w \__int_eval:w 9999 - #5 \exp_after:wN ;
    \__int_value:w \__int_eval:w 9999 - #6 \exp_after:wN ;
    \__int_value:w \__int_eval:w 1 0000 - #7 ;
  }
%    \end{macrocode}
%
%   \begin{quote}
%     \cs{@@_ln_t_large:NNw} \meta{sign}\meta{fixed tl}  \meta{t_1}|;| \meta{t_2} |;| \meta{t_3}|;| \meta{t_4}|;| \meta{t_5} |;| \meta{t_6}|;| \meta{exponent} |;| \meta{continuation}
%   \end{quote}
%   Compute the square $|t|^2$, and keep $|t|$ at the end with its
%   sign. We know that $|t|<0.1765$, so every piece has at most $4$
%   digits. However, since we were not careful in \cs{@@_ln_t_small:w},
%   they can have less than $4$ digits.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_t_large:NNw #1 #2 #3; #4; #5; #6; #7; #8;
  {
    \exp_after:wN \@@_ln_square_t_after:w
    \__int_value:w \__int_eval:w 9999 0000 + #3*#3
      \exp_after:wN \@@_ln_square_t_pack:NNNNNw
      \__int_value:w \__int_eval:w 9999 0000 + 2*#3*#4
        \exp_after:wN \@@_ln_square_t_pack:NNNNNw
        \__int_value:w \__int_eval:w 9999 0000 + 2*#3*#5 + #4*#4
          \exp_after:wN \@@_ln_square_t_pack:NNNNNw
          \__int_value:w \__int_eval:w 9999 0000 + 2*#3*#6 + 2*#4*#5
            \exp_after:wN \@@_ln_square_t_pack:NNNNNw
            \__int_value:w \__int_eval:w 1 0000 0000 + 2*#3*#7 + 2*#4*#6 + #5*#5
              + (2*#3*#8 + 2*#4*#7 + 2*#5*#6) / 1 0000
              % ; ; ;
    \exp_after:wN \@@_ln_twice_t_after:w
    \__int_value:w \__int_eval:w -1 + 2*#3
      \exp_after:wN \@@_ln_twice_t_pack:Nw
      \__int_value:w \__int_eval:w 9999 + 2*#4
        \exp_after:wN \@@_ln_twice_t_pack:Nw
        \__int_value:w \__int_eval:w 9999 + 2*#5
          \exp_after:wN \@@_ln_twice_t_pack:Nw
          \__int_value:w \__int_eval:w 9999 + 2*#6
            \exp_after:wN \@@_ln_twice_t_pack:Nw
            \__int_value:w \__int_eval:w 9999 + 2*#7
              \exp_after:wN \@@_ln_twice_t_pack:Nw
              \__int_value:w \__int_eval:w 10000 + 2*#8 ; ;
    { \@@_ln_c:NwNw #1 }
    #2
  }
\cs_new:Npn \@@_ln_twice_t_pack:Nw #1 #2; { + #1 ; {#2} }
\cs_new:Npn \@@_ln_twice_t_after:w #1; { ;;; {#1} }
\cs_new:Npn \@@_ln_square_t_pack:NNNNNw #1 #2#3#4#5 #6;
  { + #1#2#3#4#5 ; {#6} }
\cs_new:Npn \@@_ln_square_t_after:w 1 0 #1#2#3 #4;
  { \@@_ln_Taylor:wwNw {0#1#2#3} {#4} }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_ln_Taylor:wwNw}
%   Denoting $T=t^2$, we get
%   \begin{quote}
%     \cs{@@_ln_Taylor:wwNw}
%     \Arg{T_1} \Arg{T_2} \Arg{T_3} \Arg{T_4} \Arg{T_5} \Arg{T_6} |;| |;|
%     \Arg{(2t)_1} \Arg{(2t)_2} \Arg{(2t)_3} \Arg{(2t)_4} \Arg{(2t)_5} \Arg{(2t)_6} |;|
%     |{| \cs{@@_ln_c:NwNw} \meta{sign} |}|
%     \meta{fixed tl} \meta{exponent} |;| \meta{continuation}
%   \end{quote}
%   And we want to compute
%   \[
%   \ln\left(\frac{1+t}{1-t}\right)
%   = 2t\left(1 + T \left(\frac{1}{3} + T \left(\frac{1}{5}
%         + T \left(\frac{1}{7} + T \left( \frac{1}{9} + \cdots
%           \right)\right)\right)\right)\right)
%   \]
%   The process looks as follows
%   \begin{verbatim}
%     \loop 5; A;
%     \div_int 5; 1.0; \add A; \mul T; {\loop \eval 5-2;}
%     \add 0.2; A; \mul T; {\loop \eval 5-2;}
%     \mul B; T; {\loop 3;}
%     \loop 3; C;
%   \end{verbatim}
%   ^^A todo: doc
%
%   This uses the routine for dividing a number by a small integer
%   (${}<10^4$).
%    \begin{macrocode}
\cs_new:Npn \@@_ln_Taylor:wwNw
  { \@@_ln_Taylor_loop:www 21 ; {0000}{0000}{0000}{0000}{0000}{0000} ; }
\cs_new:Npn \@@_ln_Taylor_loop:www #1; #2; #3;
  {
    \if_int_compare:w #1 = 1 \exp_stop_f:
      \@@_ln_Taylor_break:w
    \fi:
    \exp_after:wN \@@_fixed_div_int:wwN \c_@@_one_fixed_tl #1;
    \@@_fixed_add:wwn #2;
    \@@_fixed_mul:wwn #3;
    {
      \exp_after:wN \@@_ln_Taylor_loop:www
      \__int_value:w \__int_eval:w #1 - 2 ;
    }
    #3;
  }
\cs_new:Npn \@@_ln_Taylor_break:w \fi: #1 \@@_fixed_add:wwn #2#3; #4 ;;
  {
    \fi:
    \exp_after:wN \@@_fixed_mul:wwn
    \exp_after:wN { \__int_value:w \__int_eval:w 10000 + #2 } #3;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_ln_c:NwNw}
%   \begin{quote}
%     \cs{@@_ln_c:NwNw} \meta{sign}
%     \Arg{r_1} \Arg{r_2} \Arg{r_3} \Arg{r_4} \Arg{r_5} \Arg{r_6} |;|
%     \meta{fixed tl} \meta{exponent} |;| \meta{continuation}
%   \end{quote}
%   We are now reduced to finding $\ln(c)$ and $\meta{exponent}\ln(10)$
%   in a table, and adding it to the mixture. The first step is to
%   get $\ln(c) - \ln(x) = - \ln(a)$, then we get $|b|\ln(10)$ and add
%   or subtract.
%
%   For now, $\ln(x)$ is given as $\cdot 10^0$. Unless both the exponent
%   is $1$ and $c=1$, we shift to working in units of $\cdot 10^4$,
%   since the final result will be at least $\ln(10/7) \simeq
%   0.35$.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_c:NwNw #1 #2; #3
  {
    \if_meaning:w + #1
      \exp_after:wN \exp_after:wN \exp_after:wN \@@_fixed_sub:wwn
    \else:
      \exp_after:wN \exp_after:wN \exp_after:wN \@@_fixed_add:wwn
    \fi:
    #3 #2 ;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}{\@@_ln_exponent:wn}
%   \begin{quote}\raggedright
%     \cs{@@_ln_exponent:wn}
%     \Arg{s_1} \Arg{s_2} \Arg{s_3} \Arg{s_4} \Arg{s_5} \Arg{s_6} |;|
%     \Arg{exponent}
%   \end{quote}
%   Compute \meta{exponent} times $\ln(10)$. Apart from the cases where
%   \meta{exponent} is $0$ or $1$, the result will necessarily be at
%   least $\ln(10) \simeq 2.3$ in magnitude. We can thus drop the least
%   significant $4$ digits. In the case of a very large (positive or
%   negative) exponent, we can (and we need to) drop $4$ additional
%   digits, since the result is of order $10^4$. Naively, one would
%   think that in both cases we can drop $4$ more digits than we do,
%   but that would be slightly too tight for rounding to happen correctly.
%   Besides, we already have addition and subtraction for $24$ digits
%   fixed point numbers.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_exponent:wn #1; #2
  {
    \if_case:w #2 \exp_stop_f:
      0 \@@_case_return:nw { \@@_fixed_to_float_o:Nw 2 }
    \or:
      \exp_after:wN \@@_ln_exponent_one:ww \__int_value:w
    \else:
      \if_int_compare:w #2 > 0 \exp_stop_f:
        \exp_after:wN \@@_ln_exponent_small:NNww
        \exp_after:wN 0
        \exp_after:wN \@@_fixed_sub:wwn \__int_value:w
      \else:
        \exp_after:wN \@@_ln_exponent_small:NNww
        \exp_after:wN 2
        \exp_after:wN \@@_fixed_add:wwn \__int_value:w -
      \fi:
    \fi:
    #2; #1;
  }
%    \end{macrocode}
%   Now we painfully write all the cases.\footnote{Bruno: do rounding.}
%   No overflow nor underflow can happen, except when computing \texttt{ln(1)}.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_exponent_one:ww 1; #1;
  {
    0
    \exp_after:wN \@@_fixed_sub:wwn \c_@@_ln_x_fixed_tl #1;
    \@@_fixed_to_float_o:wN 0
  }
%    \end{macrocode}
%   For small exponents, we just drop one block of digits, and set the
%   exponent of the log to $4$ (minus any shift coming from leading zeros
%   in the conversion from fixed point to floating point). Note that here
%   the exponent has been made positive.
%    \begin{macrocode}
\cs_new:Npn \@@_ln_exponent_small:NNww #1#2#3; #4#5#6#7#8#9;
  {
    4
    \exp_after:wN \@@_fixed_mul:wwn
      \c_@@_ln_x_fixed_tl
      {#3}{0000}{0000}{0000}{0000}{0000} ;
    #2
      {0000}{#4}{#5}{#6}{#7}{#8};
    \@@_fixed_to_float_o:wN #1
  }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Exponential}
%
% \subsubsection{Sign, exponent, and special numbers}
%
% \begin{macro}[int, EXP]{\@@_exp_o:w}
%    \begin{macrocode}
\cs_new:Npn \@@_exp_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
  {
    \if_case:w #2 \exp_stop_f:
      \@@_case_return_o:Nw \c_one_fp
    \or:
      \exp_after:wN \@@_exp_normal_o:w
    \or:
      \if_meaning:w 0 #3
        \exp_after:wN \@@_case_return_o:Nw
        \exp_after:wN \c_inf_fp
      \else:
        \exp_after:wN \@@_case_return_o:Nw
        \exp_after:wN \c_zero_fp
      \fi:
    \or:
      \@@_case_return_same_o:w
    \fi:
    \s_@@ \@@_chk:w #2#3#4;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, EXP]{\@@_exp_normal_o:w, \@@_exp_pos_o:Nnwnw, \@@_exp_overflow:NN}
%   \begin{macrocode}
\cs_new:Npn \@@_exp_normal_o:w \s_@@ \@@_chk:w 1#1
  {
    \if_meaning:w 0 #1
      \@@_exp_pos_o:NNwnw + \@@_fixed_to_float_o:wN
    \else:
      \@@_exp_pos_o:NNwnw - \@@_fixed_inv_to_float_o:wN
    \fi:
  }
\cs_new:Npn \@@_exp_pos_o:NNwnw #1#2#3 \fi: #4#5;
  {
    \fi:
    \if_int_compare:w #4 > \c_@@_max_exp_exponent_int
      \token_if_eq_charcode:NNTF + #1
        { \@@_exp_overflow:NN \@@_overflow:w \c_inf_fp }
        { \@@_exp_overflow:NN \@@_underflow:w \c_zero_fp }
      \exp:w
    \else:
      \exp_after:wN \@@_sanitize:Nw
      \exp_after:wN 0
      \__int_value:w #1 \__int_eval:w
        \if_int_compare:w #4 < 0 \exp_stop_f:
          \exp_after:wN \use_i:nn
        \else:
          \exp_after:wN \use_ii:nn
        \fi:
        {
          0
          \@@_decimate:nNnnnn { - #4 }
            \@@_exp_Taylor:Nnnwn
        }
        {
          \@@_decimate:nNnnnn { \c_@@_prec_int - #4 }
            \@@_exp_pos_large:NnnNwn
        }
        #5
        {#4}
        #1 #2 0
        \exp:w
    \fi:
    \exp_after:wN \exp_end:
  }
\cs_new:Npn \@@_exp_overflow:NN #1#2
  {
    \exp_after:wN \exp_after:wN
    \exp_after:wN #1
    \exp_after:wN #2
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[int, EXP]{\@@_exp_Taylor:Nnnwn}
% \begin{macro}[aux, EXP]{\@@_exp_Taylor_loop:www, \@@_exp_Taylor_break:Nww}
%   This function is called for numbers in the range $[10^{-9},
%   10^{-1})$.  We compute $10$ terms of the Taylor series.  The
%   first argument is irrelevant (rounding digit used by some other
%   functions).  The next three arguments, at least $16$ digits,
%   delimited by a semicolon, form a fixed point number, so we pack it
%   in blocks of $4$ digits.
%    \begin{macrocode}
\cs_new:Npn \@@_exp_Taylor:Nnnwn #1#2#3 #4; #5 #6
  {
    #6
    \@@_pack_twice_four:wNNNNNNNN
    \@@_pack_twice_four:wNNNNNNNN
    \@@_pack_twice_four:wNNNNNNNN
    \@@_exp_Taylor_ii:ww
    ; #2#3#4 0000 0000 ;
  }
\cs_new:Npn \@@_exp_Taylor_ii:ww #1; #2;
  { \@@_exp_Taylor_loop:www 10 ; #1 ; #1 ; \s__stop }
\cs_new:Npn \@@_exp_Taylor_loop:www #1; #2; #3;
  {
    \if_int_compare:w #1 = 1 \exp_stop_f:
      \exp_after:wN \@@_exp_Taylor_break:Nww
    \fi:
    \@@_fixed_div_int:wwN #3 ; #1 ;
    \@@_fixed_add_one:wN
    \@@_fixed_mul:wwn #2 ;
    {
      \exp_after:wN \@@_exp_Taylor_loop:www
      \__int_value:w \__int_eval:w #1 - 1 ;
      #2 ;
    }
  }
\cs_new:Npn \@@_exp_Taylor_break:Nww #1 #2; #3 \s__stop
  { \@@_fixed_add_one:wN #2 ; }
%    \end{macrocode}
% \end{macro}
% \end{macro}
%
% \begin{macro}[aux, rEXP]
%   {
%     \@@_exp_pos_large:NnnNwn ,
%     \@@_exp_large_after:wwn ,
%     \@@_exp_large:w ,
%     \@@_exp_large_v:wN,
%     \@@_exp_large_iv:wN,
%     \@@_exp_large_iii:wN,
%     \@@_exp_large_ii:wN,
%     \@@_exp_large_i:wN,
%     \@@_exp_large_:wN,
%   }
%   The first two arguments are irrelevant (a rounding digit, and a
%   brace group with $8$ zeros).  The third argument is the integer part
%   of our number, then we have the decimal part delimited by a
%   semicolon, and finally the exponent, in the range $[0,5]$.  Remove
%   leading zeros from the integer part: putting |#4| in there too
%   ensures that an integer part of $0$ is also removed.  Then read
%   digits one by one, looking up $\exp(\meta{digit}\cdot
%   10^{\meta{exponent}})$ in a table, and multiplying that to the
%   current total.  The loop is done by having the auxiliary for one
%   exponent call the auxiliary for the next exponent.  The current
%   total is expressed by leaving the exponent behind in the input
%   stream (we are currently within an \cs{__int_eval:w}), and keeping
%   track of a fixed point number, |#1| for the numbered auxiliaries.
%   Our usage of \cs{if_case:w} is somewhat dirty for optimization:
%   \TeX{} jumps to the appropriate case, but we then close the
%   \cs{if_case:w} \enquote{by hand}, using \cs{or:} and \cs{fi:} as
%   delimiters.
%    \begin{macrocode}
\cs_new:Npn \@@_exp_pos_large:NnnNwn #1#2#3 #4#5; #6
  {
    \exp_after:wN \exp_after:wN
    \cs:w @@_exp_large_ \__int_to_roman:w #6 :wN \exp_after:wN \cs_end:
    \exp_after:wN \c_@@_one_fixed_tl
    \__int_value:w #3 #4 \exp_stop_f:
    #5 00000 ;
  }
\cs_new:Npn \@@_exp_large:w #1 \or: #2 \fi:
  { \fi: \@@_fixed_mul:wwn #1; }
\cs_new:Npn \@@_exp_large_v:wN #1; #2
  {
    \if_case:w #2 ~           \exp_after:wN \@@_fixed_continue:wn  \or:
      +  4343 \@@_exp_large:w {8806}{8182}{2566}{2921}{5872}{6150} \or:
      +  8686 \@@_exp_large:w {7756}{0047}{2598}{6861}{0458}{3204} \or:
      + 13029 \@@_exp_large:w {6830}{5723}{7791}{4884}{1932}{7351} \or:
      + 17372 \@@_exp_large:w {6015}{5609}{3095}{3052}{3494}{7574} \or:
      + 21715 \@@_exp_large:w {5297}{7951}{6443}{0315}{3251}{3576} \or:
      + 26058 \@@_exp_large:w {4665}{6719}{0099}{3379}{5527}{2929} \or:
      + 30401 \@@_exp_large:w {4108}{9724}{3326}{3186}{5271}{5665} \or:
      + 34744 \@@_exp_large:w {3618}{6973}{3140}{0875}{3856}{4102} \or:
      + 39087 \@@_exp_large:w {3186}{9209}{6113}{3900}{6705}{9685} \or:
    \fi:
    #1;
    \@@_exp_large_iv:wN
  }
\cs_new:Npn \@@_exp_large_iv:wN #1; #2
  {
    \if_case:w #2 ~          \exp_after:wN \@@_fixed_continue:wn  \or:
      +  435 \@@_exp_large:w {1970}{0711}{1401}{7046}{9938}{8888} \or:
      +  869 \@@_exp_large:w {3881}{1801}{9428}{4368}{5764}{8232} \or:
      + 1303 \@@_exp_large:w {7646}{2009}{8905}{4704}{8893}{1073} \or:
      + 1738 \@@_exp_large:w {1506}{3559}{7005}{0524}{9009}{7592} \or:
      + 2172 \@@_exp_large:w {2967}{6283}{8402}{3667}{0689}{6630} \or:
      + 2606 \@@_exp_large:w {5846}{4389}{5650}{2114}{7278}{5046} \or:
      + 3041 \@@_exp_large:w {1151}{7900}{5080}{6878}{2914}{4154} \or:
      + 3475 \@@_exp_large:w {2269}{1083}{0850}{6857}{8724}{4002} \or:
      + 3909 \@@_exp_large:w {4470}{3047}{3316}{5442}{6408}{6591} \or:
    \fi:
    #1;
    \@@_exp_large_iii:wN
  }
\cs_new:Npn \@@_exp_large_iii:wN #1; #2
  {
    \if_case:w #2 ~         \exp_after:wN \@@_fixed_continue:wn  \or:
      +  44 \@@_exp_large:w {2688}{1171}{4181}{6135}{4484}{1263} \or:
      +  87 \@@_exp_large:w {7225}{9737}{6812}{5749}{2581}{7748} \or:
      + 131 \@@_exp_large:w {1942}{4263}{9524}{1255}{9365}{8421} \or:
      + 174 \@@_exp_large:w {5221}{4696}{8976}{4143}{9505}{8876} \or:
      + 218 \@@_exp_large:w {1403}{5922}{1785}{2837}{4107}{3977} \or:
      + 261 \@@_exp_large:w {3773}{0203}{0092}{9939}{8234}{0143} \or:
      + 305 \@@_exp_large:w {1014}{2320}{5473}{5004}{5094}{5533} \or:
      + 348 \@@_exp_large:w {2726}{3745}{7211}{2566}{5673}{6478} \or:
      + 391 \@@_exp_large:w {7328}{8142}{2230}{7421}{7051}{8866} \or:
    \fi:
    #1;
    \@@_exp_large_ii:wN
  }
\cs_new:Npn \@@_exp_large_ii:wN #1; #2
  {
    \if_case:w #2 ~        \exp_after:wN \@@_fixed_continue:wn  \or:
      +  5 \@@_exp_large:w {2202}{6465}{7948}{0671}{6516}{9579} \or:
      +  9 \@@_exp_large:w {4851}{6519}{5409}{7902}{7796}{9107} \or:
      + 14 \@@_exp_large:w {1068}{6474}{5815}{2446}{2146}{9905} \or:
      + 18 \@@_exp_large:w {2353}{8526}{6837}{0199}{8540}{7900} \or:
      + 22 \@@_exp_large:w {5184}{7055}{2858}{7072}{4640}{8745} \or:
      + 27 \@@_exp_large:w {1142}{0073}{8981}{5684}{2836}{6296} \or:
      + 31 \@@_exp_large:w {2515}{4386}{7091}{9167}{0062}{6578} \or:
      + 35 \@@_exp_large:w {5540}{6223}{8439}{3510}{0525}{7117} \or:
      + 40 \@@_exp_large:w {1220}{4032}{9431}{7840}{8020}{0271} \or:
    \fi:
    #1;
    \@@_exp_large_i:wN
  }
\cs_new:Npn \@@_exp_large_i:wN #1; #2
  {
    \if_case:w #2 ~       \exp_after:wN \@@_fixed_continue:wn  \or:
      + 1 \@@_exp_large:w {2718}{2818}{2845}{9045}{2353}{6029} \or:
      + 1 \@@_exp_large:w {7389}{0560}{9893}{0650}{2272}{3043} \or:
      + 2 \@@_exp_large:w {2008}{5536}{9231}{8766}{7740}{9285} \or:
      + 2 \@@_exp_large:w {5459}{8150}{0331}{4423}{9078}{1103} \or:
      + 3 \@@_exp_large:w {1484}{1315}{9102}{5766}{0342}{1116} \or:
      + 3 \@@_exp_large:w {4034}{2879}{3492}{7351}{2260}{8387} \or:
      + 4 \@@_exp_large:w {1096}{6331}{5842}{8458}{5992}{6372} \or:
      + 4 \@@_exp_large:w {2980}{9579}{8704}{1728}{2747}{4359} \or:
      + 4 \@@_exp_large:w {8103}{0839}{2757}{5384}{0077}{1000} \or:
    \fi:
    #1;
    \@@_exp_large_:wN
  }
\cs_new:Npn \@@_exp_large_:wN #1; #2
  {
    \if_case:w #2 ~       \exp_after:wN \@@_fixed_continue:wn  \or:
      + 1 \@@_exp_large:w {1105}{1709}{1807}{5647}{6248}{1171} \or:
      + 1 \@@_exp_large:w {1221}{4027}{5816}{0169}{8339}{2107} \or:
      + 1 \@@_exp_large:w {1349}{8588}{0757}{6003}{1039}{8374} \or:
      + 1 \@@_exp_large:w {1491}{8246}{9764}{1270}{3178}{2485} \or:
      + 1 \@@_exp_large:w {1648}{7212}{7070}{0128}{1468}{4865} \or:
      + 1 \@@_exp_large:w {1822}{1188}{0039}{0508}{9748}{7537} \or:
      + 1 \@@_exp_large:w {2013}{7527}{0747}{0476}{5216}{2455} \or:
      + 1 \@@_exp_large:w {2225}{5409}{2849}{2467}{6045}{7954} \or:
      + 1 \@@_exp_large:w {2459}{6031}{1115}{6949}{6638}{0013} \or:
    \fi:
    #1;
    \@@_exp_large_after:wwn
  }
\cs_new:Npn \@@_exp_large_after:wwn #1; #2; #3
  {
    \@@_exp_Taylor:Nnnwn ? { } { } 0 #2; {} #3
    \@@_fixed_mul:wwn #1;
  }
%    \end{macrocode}
% \end{macro}
%
% \subsection{Power}
%
% Raising a number $a$ to a power $b$ leads to many distinct situations.
% \begin{center}\def\abs#1{\lvert #1\rvert}
%   \begin{tabular}{>{$}c<{$}|*8{>{$}l<{$}}}
%     a^b          &-\infty &(-\infty,-0)  &-p/5^k              &\pm 0 &+p/5^k            &(0,\infty)      &+\infty &\nan \\ \hline
%     +\infty      &+0      &\multicolumn{2}{c}{$+0$}           &+1    &\multicolumn{2}{c}{$+\infty$}      &+\infty &\nan \\
%     (1,\infty)   &+0      &\multicolumn{2}{c}{$+\abs{a}^{b}$} &+1    &\multicolumn{2}{c}{$+\abs{a}^{b}$} &+\infty &\nan \\
%     +1           &+1      &\multicolumn{2}{c}{$+1$}           &+1    &\multicolumn{2}{c}{$+1$}           &+1      &+1   \\
%     (0,1)        &+\infty &\multicolumn{2}{c}{$+\abs{a}^{b}$} &+1    &\multicolumn{2}{c}{$+\abs{a}^{b}$} &+0      &\nan \\
%     +0           &+\infty &\multicolumn{2}{c}{$+\infty$}      &+1    &\multicolumn{2}{c}{$+0$}           &+0      &\nan \\
%     -0           &+\infty &\nan          &(-1)^p\infty        &+1    &(-1)^p 0          &+0              &+0      &\nan \\
%     (-1,0)       &+\infty &\nan          &(-1)^p\abs{a}^{b}   &+1    &(-1)^p\abs{a}^{b} &\nan            &+0      &\nan \\
%     -1           &+1      &\nan          &(-1)^p              &+1    &(-1)^p            &\nan            &+1      &\nan \\
%     (-\infty,-1) &+0      &\nan          &(-1)^p\abs{a}^{b}   &+1    &(-1)^p\abs{a}^{b} &\nan            &+\infty &\nan \\
%     -\infty      &+0      &+0            &(-1)^p 0            &+1    &(-1)^p\infty      &\nan            &+\infty &\nan \\
%     \nan         &\nan    &\nan          &\nan                &+1    &\nan              &\nan            &\nan    &\nan \\
%   \end{tabular}
% \end{center}
% We distinguished in this table the cases of finite (positive or
% negative) exponents of the form $b=p/q$ with $q$~odd (hence
% necessarily a power of~$5$), as $(-1)^{p/q}=(-1)^p$ is defined in that
% case.
% One peculiarity of this operation is that $\nan^0 = 1^\nan = 1$,
% because this relation is obeyed for any number, even $\pm\infty$.
%
% \begin{macro}[int, EXP]+\@@_^_o:ww+
%   We cram most of the tests into a single function to save csnames.
%   First treat the case $b=0$: $a^0=1$ for any $a$, even \texttt{nan}.
%   Then test the sign of $a$.
%   \begin{itemize}
%   \item If it is positive, and $a$ is a normal number, call
%     \cs{@@_pow_normal_o:ww} followed by the two \texttt{fp} $a$ and $b$.
%     For $a=+0$ or $+\inf$, call \cs{@@_pow_zero_or_inf:ww} instead, to
%     return either $+0$ or $+\infty$ as appropriate.
%   \item If $a$ is a \texttt{nan}, then skip to the next semicolon
%     (which happens to be conveniently the end of $b$) and return
%     \texttt{nan}.
%   \item Finally, if $a$ is negative, compute $|a|^b$
%     (\cs{@@_pow_normal_o:ww} which ignores the sign of its first
%     operand), and keep an extra copy of $a$ and $b$ (the second brace
%     group, containing \{~$b$~$a$~\}, is inserted between $a$ and $b$).
%     Then do some tests to find the final sign of the result if it
%     exists.
%   \end{itemize}
%    \begin{macrocode}
\cs_new:cpn { @@_ \iow_char:N \^ _o:ww }
    \s_@@ \@@_chk:w #1#2#3; \s_@@ \@@_chk:w #4#5#6;
  {
    \if_meaning:w 0 #4
      \@@_case_return_o:Nw \c_one_fp
    \fi:
    \if_case:w #2 \exp_stop_f:
      \exp_after:wN \use_i:nn
    \or:
      \@@_case_return_o:Nw \c_nan_fp
    \else:
      \exp_after:wN \@@_pow_neg:www
      \exp:w \exp_end_continue_f:w \exp_after:wN \use:nn
    \fi:
    {
      \if_meaning:w 1 #1
        \exp_after:wN \@@_pow_normal_o:ww
      \else:
        \exp_after:wN \@@_pow_zero_or_inf:ww
      \fi:
      \s_@@ \@@_chk:w #1#2#3;
    }
    { \s_@@ \@@_chk:w #4#5#6; \s_@@ \@@_chk:w #1#2#3; }
    \s_@@ \@@_chk:w #4#5#6;
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, EXP]{\@@_pow_zero_or_inf:ww}
%   Raising $-0$ or $-\infty$ to \texttt{nan} yields \texttt{nan}.  For
%   other powers, the result is $+0$ if $0$ is raised to a positive
%   power or $\infty$ to a negative power, and $+\infty$ otherwise.
%   Thus, if the type of $a$ and the sign of $b$ coincide, the result
%   is~$0$, since those conveniently take the same possible values, $0$
%   and~$2$.  Otherwise, either $a=\pm\infty$ and $b>0$ and the result
%   is $+\infty$, or $a=\pm 0$ with $b<0$ and we have a division by zero
%   unless $b=-\infty$.
%    \begin{macrocode}
\cs_new:Npn \@@_pow_zero_or_inf:ww
    \s_@@ \@@_chk:w #1#2; \s_@@ \@@_chk:w #3#4
  {
    \if_meaning:w 1 #4
      \@@_case_return_same_o:w
    \fi:
    \if_meaning:w #1 #4
      \@@_case_return_o:Nw \c_zero_fp
    \fi:
    \if_meaning:w 2 #1
      \@@_case_return_o:Nw \c_inf_fp
    \fi:
    \if_meaning:w 2 #3
      \@@_case_return_o:Nw \c_inf_fp
    \else:
      \@@_case_use:nw
        {
          \@@_division_by_zero_o:NNww \c_inf_fp ^
            \s_@@ \@@_chk:w #1 #2 ;
        }
    \fi:
    \s_@@ \@@_chk:w #3#4
  }
%    \end{macrocode}
% \end{macro}
%
% \begin{macro}[aux, EXP]{\@@_pow_normal_o:ww}
%   We have in front of us $a$, and $b\neq 0$, we know that $a$ is a
%   normal number, and we wish to compute $\lvert a\rvert^{b}$.  If
%   $\lvert a\rvert=1$, we return $1$, unless $a=-1$ and $b$ is
%   \texttt{nan}.  Indeed, returning $1$ at this point would wrongly
%   raise \enquote{invalid} when the sign is considered.  If $\lvert
%   a\rvert\neq 1$, test the type of $b$:
%   \begin{itemize}
%   \item[0] Impossible, we already filtered $b=\pm 0$.
%   \item[1] Call \cs{@@_pow_npos_o:Nww}.
%   \item[2] Return $+\infty$ or $+0$ depending on the sign of $b$ and
%     whether the exponent of $a$ is positive or not.
%   \item[3] Return $b$.
%   \end{itemize}
%    \begin{macrocode}
\cs_new:Npn \@@_pow_normal_o:ww
    \s_@@ \@@_chk:w 1 #1#2#3; \s_@@ \@@_chk:w #4#5
  {
    \if_int_compare:w \__str_if_eq_x:nn { #2 #3 }
              { 1 {1000} {0000} {0000} {0000} } = 0 \exp_stop_f:
      \if_int_compare:w #4 #1 = 32 \exp_stop_f:
        \exp_after:wN \@@_case_return_ii_o:ww
      \fi:
      \@@_case_return_o:Nww \c_one_fp
    \fi:
    \if_case:w #4 \exp_stop_f:
    \or:
      \exp_after:wN \@@_pow_npos_o:Nww
      \exp_after:wN #5
    \or:
      \if_meaning:w 2 #5 \exp_after:wN \reverse_if:N \fi:
      \if_int_compare:w #2 > 0 \exp_stop_f:
        \exp_after:wN \@@_case_return_o:Nww
        \exp_after:wN \c_inf_fp
      \else:
        \exp_after:wN \@@_case_return_o:Nww
        \exp_after:wN \c_zero_fp
      \fi:
    \or:
      \@@_case_return_ii_o:ww
    \fi:
    \s_@@ \@@_chk:w 1 #1 {#2} #3 ;
    \s_@@ \@@_chk:w #4 #5
  }
%    \end{macrocode}
% \end{macro}
%
% ^^A todo: check that we compute ln to 21 digits!
% \begin{macro}[aux, EXP]{\@@_pow_npos_o:Nww}
%   We now know that $a\neq\pm 1$ is a normal number, and $b$ is a
%   normal number too.  We want to compute $\lvert a\rvert^{b} = (\lvert
%   x\rvert\cdot 10^{n})^{y\cdot 10^{p}} = \exp((\ln\lvert x\rvert + n
%   \ln(10))\cdot y \cdot 10^{p}) = \exp(z)$.  To compute the
%   exponential accurately, we need to know the digits of $z$ up to the
%   $16$-th position.  Since the exponential of $10^{5}$ is infinite, we
%   only need at most $21$ digits, hence the fixed point result of
%   \cs{@@_ln_o:w} is precise enough for our needs.  Start an integer
%   expression for the decimal exponent of $e^{\lvert z\rvert}$.  If $z$
%   is negative, negate that decimal exponent, and prepare to take the
%   inverse when converting from the fixed point to the floating point result.
%    \begin{macrocode}
\cs_new:Npn \@@_pow_npos_o:Nww #1 \s_@@ \@@_chk:w 1#2#3
  {
    \exp_after:wN \@@_sanitize:Nw
    \exp_after:wN 0
    \__int_value:w
      \if:w #1 \if_int_compare:w #3 > 0 \exp_stop_f: 0 \else: 2 \fi:
        \exp_after:wN \@@_pow_npos_aux:NNnww
        \exp_after:wN +
        \exp_after:wN \@@_fixed_to_float_o:wN
      \else:
        \exp_after:wN \@@_pow_npos_aux:NNnww
        \exp_after:wN -
        \exp_after:wN \@@_fixed_inv_to_float_o:wN
      \fi:
      {#3}
  }
%    \end{macrocode}
% \end{macro}
%
%^^A begin[todo]
% \begin{macro}[aux, EXP]{\@@_pow_npos_aux:NNnww}
%   The first argument is the conversion function from fixed point to
%   float.  Then comes an exponent and the $4$ brace groups of $x$,
%   followed by $b$.  Compute $-\ln(x)$.
%    \begin{macrocode}
\cs_new:Npn \@@_pow_npos_aux:NNnww #1#2#3#4#5; \s_@@ \@@_chk:w 1#6#7#8;
  {
    #1
    \__int_eval:w
      \@@_ln_significand:NNNNnnnN #4#5
      \@@_pow_exponent:wnN {#3}
      \@@_fixed_mul:wwn #8 {0000}{0000} ;
      \@@_pow_B:wwN #7;
      #1 #2 0 % fixed_to_float_o:wN
  }
\cs_new:Npn \@@_pow_exponent:wnN #1; #2
  {
    \if_int_compare:w #2 > 0 \exp_stop_f:
      \exp_after:wN \@@_pow_exponent:Nwnnnnnw % n\ln(10) - (-\ln(x))
      \exp_after:wN +
    \else:
      \exp_after:wN \@@_pow_exponent:Nwnnnnnw % -(|n|\ln(10) + (-\ln(x)))
      \exp_after:wN -
    \fi:
    #2; #1;
  }
\cs_new:Npn \@@_pow_exponent:Nwnnnnnw #1#2; #3#4#5#6#7#8;
  { %^^A todo: use that in ln.
    \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#2*23025 - #1 #3
        \exp_after:wN \@@_pack:NNNNNw
        \__int_value:w \__int_eval:w \c_@@_middle_shift_int
          #1 #2*8509 - #1 #4
          \exp_after:wN \@@_pack:NNNNNw
          \__int_value:w \__int_eval:w \c_@@_middle_shift_int
            #1 #2*2994 - #1 #5
            \exp_after:wN \@@_pack:NNNNNw
            \__int_value:w \__int_eval:w \c_@@_middle_shift_int
              #1 #2*0456 - #1 #6
              \exp_after:wN \@@_pack:NNNNNw
              \__int_value:w \__int_eval:w \c_@@_trailing_shift_int
                #1 #2*8401 - #1 #7
                #1 ( #2*7991 - #8 ) / 1 0000 ; ;
  }
\cs_new:Npn \@@_pow_B:wwN #1#2#3#4#5#6; #7;
  {
    \if_int_compare:w #7 < 0 \exp_stop_f:
      \exp_after:wN \@@_pow_C_neg:w \__int_value:w -
    \else:
      \if_int_compare:w #7 < 22 \exp_stop_f:
        \exp_after:wN \@@_pow_C_pos:w \__int_value:w
      \else:
        \exp_after:wN \@@_pow_C_overflow:w \__int_value:w
      \fi:
    \fi:
    #7 \exp_after:wN ;
    \__int_value:w \__int_eval:w 10 0000 + #1 \__int_eval_end:
    #2#3#4#5#6 0000 0000 0000 0000 0000 0000 ; %^^A todo: how many 0?
  }
\cs_new:Npn \@@_pow_C_overflow:w #1; #2; #3
  {
    + 2 * \c_@@_max_exponent_int
    \exp_after:wN \@@_fixed_continue:wn \c_@@_one_fixed_tl
  }
\cs_new:Npn \@@_pow_C_neg:w #1 ; 1
  {
    \exp_after:wN \exp_after:wN \exp_after:wN \@@_pow_C_pack:w
    \prg_replicate:nn {#1} {0}
  }
\cs_new:Npn \@@_pow_C_pos:w #1; 1
  { \@@_pow_C_pos_loop:wN #1; }
\cs_new:Npn \@@_pow_C_pos_loop:wN #1; #2
  {
    \if_meaning:w 0 #1
      \exp_after:wN \@@_pow_C_pack:w
      \exp_after:wN #2
    \else:
      \if_meaning:w 0 #2
        \exp_after:wN \@@_pow_C_pos_loop:wN \__int_value:w
      \else:
        \exp_after:wN \@@_pow_C_overflow:w \__int_value:w
      \fi:
      \__int_eval:w #1 - 1 \exp_after:wN ;
    \fi:
  }
\cs_new:Npn \@@_pow_C_pack:w
  { \exp_after:wN \@@_exp_large_v:wN \c_@@_one_fixed_tl }
%    \end{macrocode}
% \end{macro}
%^^A end[todo]
%
% \begin{macro}[aux, EXP]{\@@_pow_neg:www, \@@_pow_neg_aux:wNN}
%   This function is followed by three floating point numbers: $|a|^b$,
%   $a\in[-\infty,-0]$, and $b$.  If $b$ is an even integer (case $-1$),
%   $a^b=|a|^b$.  If $b$ is an odd integer (case $0$), $a^b=-|a|^b$,
%   obtained by a call to \cs{@@_pow_neg_aux:wNN}.  Otherwise, the sign is
%   undefined.  This is invalid, unless $|a|^b$ turns out to be $+0$ or
%   \texttt{nan}, in which case we return that as $a^b$.  In particular,
%   since the underflow detection occurs before \cs{@@_pow_neg:www} is
%   called, |(-0.1)**(12345.67)| will give $+0$ rather than complaining
%   that the sign is not defined.
%    \begin{macrocode}
\cs_new:Npn \@@_pow_neg:www \s_@@ \@@_chk:w #1#2; #3; #4;
  {
    \if_case:w \@@_pow_neg_case:w #4 ;
      \exp_after:wN \@@_pow_neg_aux:wNN
    \or:
      \if_int_compare:w \__int_eval:w #1 / 2 = 1 \exp_stop_f:
        \@@_invalid_operation_o:Nww ^ #3; #4;
        \exp:w \exp_end_continue_f:w
        \exp_after:wN \exp_after:wN
        \exp_after:wN \@@_use_none_until_s:w
      \fi:
    \fi:
    \@@_exp_after_o:w
    \s_@@ \@@_chk:w #1#2;
  }
\cs_new:Npn \@@_pow_neg_aux:wNN #1 \s_@@ \@@_chk:w #2#3
  {
    \exp_after:wN \@@_exp_after_o:w
    \exp_after:wN \s_@@
    \exp_after:wN \@@_chk:w
    \exp_after:wN #2
    \__int_value:w \__int_eval:w 2 - #3 \__int_eval_end:
  }
%    \end{macrocode}
% ^^A todo: is this \@@_exp_after_o:w necessary?  Appropriate?
% \end{macro}
%
% \begin{macro}[aux, rEXP]
%   {
%     \@@_pow_neg_case:w, \@@_pow_neg_case_aux:nnnnn,
%     \@@_pow_neg_case_aux:w
%   }
%   This function expects a floating point number, and determines its
%   \enquote{parity}.  It should be used after \cs{if_case:w} or in an
%   integer expression.  It gives $-1$ if the number is an even integer
%   divided by some power of~$5$, $0$~if the number is an odd integer
%   divided by some power of~$5$, and $1$~otherwise.  Zeros and
%   $\pm\infty$ are even (because very large finite floating points are
%   even), while \texttt{nan} is a non-integer.  The sign of normal
%   numbers is irrelevant to parity.  The idea is to repeatedly multiply
%   the number by~$5$ (by halving the mantissa and shifting the
%   exponent) until the mantissa is odd (this can only happen at most
%   $53$ times since $2^{54}>10^{16}$): if the resulting exponent is
%   larger than $16$ the parity is even, if it is exactly $16$ the
%   parity is odd, and otherwise we should return~$1$.  Of course there
%   is a shortcut: we stop as soon as the exponent exceeds~$16$.
%    \begin{macrocode}
\cs_new:Npn \@@_pow_neg_case:w \s_@@ \@@_chk:w #1#2#3;
  {
    \if_case:w #1 \exp_stop_f:
           -1
    \or:   \@@_pow_neg_case_aux:nnnnn #3
    \or:   -1
    \else: 1
    \fi:
    \exp_stop_f:
  }
\cs_new:Npn \@@_pow_neg_case_aux:nnnnn #1#2#3#4#5
  { \@@_pow_neg_case_aux:w #1 ; #2 #3 ; #4 #5 ; }
\cs_new:Npn \@@_pow_neg_case_aux:w #1 ; #2 ; #3 ;
  {
    \if_int_compare:w #1 > \c_@@_prec_int
      -1
    \else:
      \if_int_odd:w #3 \exp_stop_f:
        \if_int_compare:w #1 = \c_@@_prec_int
          0
        \else:
          1
        \fi:
      \else:
        \exp_after:wN \@@_pow_neg_case_aux:w
        \__int_value:w \__int_eval:w #1 + 1 \exp_after:wN ;
        \__int_value:w \__int_eval:w (#2 + 1) / 2 - 1 \exp_after:wN ;
        \__int_value:w \__int_eval:w
          \if_int_odd:w #2 \exp_stop_f: 5000 0000 + \fi: #3 / 2 ;
      \fi:
    \fi:
  }
%    \end{macrocode}
% \end{macro}
%
%    \begin{macrocode}
%</initex|package>
%    \end{macrocode}
%
% \end{implementation}
%
% \PrintChanges
%
% \PrintIndex