summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx')
-rw-r--r--Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx1252
1 files changed, 1252 insertions, 0 deletions
diff --git a/Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx b/Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx
new file mode 100644
index 00000000000..edf9ac2ad35
--- /dev/null
+++ b/Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx
@@ -0,0 +1,1252 @@
+% \iffalse meta-comment
+%
+%% File: l3fp-expo.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-expo.dtx 3986 2012-07-15 19:23:51Z joseph $
+ {L3 Floating-point exponential-related functions}
+\documentclass[full]{l3doc}
+\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}
+%
+% \subsection{General comments}
+%
+%^^A todo: redoc
+% The algorithm for computing the logarithm of the significand could be
+% made to use a $5$ terms Taylor series instead of $10$ terms by taking
+% $c = 2000/(\lfloor 200x\rfloor +1) \in [10,95]$ instead of $c\in
+% [1,10]$. Also, it would then be possible to simplify the computation
+% of $t$, using methods similar to \cs{@@_fixed_div_to_float:ww}.
+% However, we would then have to hard-code the logarithms of $44$ small
+% integers instead of $9$.
+%
+% \subsection{Some constants}
+%
+% ^^A todo: very odd: can I delete ln_ten?
+% \begin{variable}
+% {
+% \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,
+% \c_@@_ln_ten_fixed_tl,
+% }
+% A few values of the logarithm which are needed in the
+% implementation. It turns out that we don't need the value of
+% $\log(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_v_fixed_tl {{16094}{3791}{2434}{1003}{7460}{0759} }
+\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} }
+ %\tl_const:Nn \c_@@_ln_ten_fixed_tl {{23025}{8509}{2994}{0456}{8401}{7991} }
+% \end{macrocode}
+% \end{variable}
+%
+% \subsection{Logarithm}
+%
+% \subsubsection{Sign, exponent, and special numbers}
+%
+% \begin{macro}[EXP]{\@@_ln:w}
+% The logarithm of $\pm 0$ is $-\infty$, raising a division by zero
+% exception. The logarithm of negative numbers (including $-\infty$,
+% but not $-0$) raises the \enquote{invalid} exception. The logarithm
+% of $+\infty$ or a \texttt{nan} is itself. Positive normal numbers
+% call \cs{@@_ln_npos:w}.
+% \begin{macrocode}
+\cs_new:Npn \@@_ln:w \s_@@ \@@_chk:w #1 #2
+ {
+ \if_meaning:w 0 #1
+ \@@_case_use:nw
+ { \@@_division_by_zero:Nnw \c_minus_inf_fp { ln } }
+ \fi:
+ \if_meaning:w 2 #2
+ \@@_case_use:nw
+ { \@@_invalid_operation:Nnw \c_nan_fp { ln } }
+ \fi:
+ \if_meaning:w 1 #1 \else:
+ \@@_case_return_same_o:w
+ \fi:
+ \@@_ln_npos:w \s_@@ \@@_chk:w #1#2
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsubsection{Absolute ln}
+%
+% 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
+% \[
+% \ln (a \cdot 10^b) = b \cdot \ln (10) - \ln (c/5) + \ln (ac/5).
+% \]
+% The logarithms $\ln(10)$ and $\ln(c/5)$ are looked up in a table. The
+% last term is computed using the following Talor series of $\ln$ near
+% $1$:
+% \[
+% \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)
+% \]
+% 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
+% \[
+% t = \frac{\epsilon}{1-\epsilon}
+% = \epsilon (1+\epsilon)(1+\epsilon^2)(1+\epsilon^4)\ldots,
+% \]
+% is not too difficult to compute.
+%
+% \begin{macro}{\@@_ln_npos: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: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 < \c_one
+ 2
+ \else:
+ 0
+ \fi:
+ \exp_after:wN \exp_stop_f:
+ \int_use:N \__int_eval:w % for the exponent
+ \@@_ln_significand:NNNNnnnN #2#3
+ \@@_ln_exponent:wn {#1}
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[int, EXP]{\@@_ln_significand:NNNNnnnN}
+% \begin{syntax}
+% \cs{@@_ln_significand:NNNNnnnN} \meta{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4} \meta{continuation}
+% \end{syntax}
+% 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 = - \log X$.
+% \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 < \c_four
+ \__int_eval:w \c_ten - #2
+ \else:
+ 6
+ \fi:
+ \or: 4
+ \or: 3
+ \or: 2
+ \or: 2
+ \or: 2
+ \else: 1
+ \fi:
+ ; { #1 #2 #3 #4 }
+ }
+% \end{macrocode}
+% We have thus found $c$. It is chosen 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_ \tex_romannumeral:D #1 _fixed_tl \exp_after:wN \cs_end:
+ \__int_value:w
+ \exp_after:wN \@@_ln_x_iv:nnnnnnnn
+ \tex_romannumeral:D -`0
+ \exp_after:wN \@@_ln_x_iii_var:NNNNNw
+ \int_use:N \__int_eval:w 9999 9999 + #1*#2#3 +
+ \exp_after:wN \@@_ln_x_iii:NNNNNw
+ \int_use:N \__int_eval:w 1 0000 0000 + #1*#4#5 ;
+ {20000} {0000} {0000} {0000}
+ } %^^A todo: reoptimize (a generalization attempt failed).
+\cs_new:Npn \@@_ln_x_iii:NNNNNw #1 #2#3#4#5 #6; { #1; {#2#3#4#5} {#6} }
+\cs_new:Npn \@@_ln_x_iii_var:NNNNNw #1 #2#3#4#5 #6; { {#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
+% occured 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).\footnote{Bruno: to be completed.}
+%
+% \begin{quote}
+% \cs{@@_ln_x_iv:NNNNNwnn}
+% \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:nnnnnnnn #1#2#3#4 #5#6#7#8
+ {
+ \exp_after:wN \@@_ln_div_i:w
+ \int_use:N \__int_eval:w #1 + \c_one ;
+ #5 #6 ; {#7} {#8}
+ {#1} {#2} {#3} {#4}
+ }
+\cs_new:Npn \@@_ln_div_i:w #1;
+ {
+ \exp_after:wN \@@_ln_div_ii:www
+ \__int_value:w #1 \exp_after:wN ;
+ \__int_value:w
+ \exp_after:wN \@@_div_mantissa_calc:Nwwnnnnnn
+ \int_use:N \__int_eval:w 999999 + 2 0000 0000 / #1 ; % Q1
+ }
+\cs_set_protected:Npn \@@_tmp:w #1#2
+ {
+ \cs_new:Npn #1 ##1; ##2; ##3; % y; Q1; B1B2; <- for k=1
+ {
+ \exp_after:wN \@@_div_mantissa_pack:NNN
+ \int_use:N \__int_eval:w ##2
+ \exp_after:wN #2
+ \__int_value:w ##1 \exp_after:wN ;
+ \__int_value:w
+ \exp_after:wN \@@_div_mantissa_calc:Nwwnnnnnn
+ \int_use:N \__int_eval:w 999999 + ##3 / ##1 ; % Q2
+ ##3 ;
+ }
+ }
+\@@_tmp:w \@@_ln_div_ii:www \@@_ln_div_iii:www
+\@@_tmp:w \@@_ln_div_iii:www \@@_ln_div_iv:www
+\@@_tmp:w \@@_ln_div_iv:www \@@_ln_div_v:www
+\@@_tmp:w \@@_ln_div_v:www \@@_ln_div_vi:www
+\cs_new:Npn \@@_ln_div_vi:www #1; #2; #3;#4#5 #6#7#8#9 %y;Q5;F1F2;F3F4x1x2x3x4
+ {
+ \exp_after:wN \@@_div_mantissa_pack:NNN
+ \int_use:N \__int_eval:w #2
+ \exp_after:wN \@@_div_mantissa_pack:NNN
+ \int_use:N \__int_eval:w 1000000 + #3 / #1 ; % Q6
+ }
+% \end{macrocode}
+% We now have essentially\footnote{Bruno: add a mention that
+% the error on $Q_{6}$ is bounded by $10$ (probably $6.7$),
+% and thus corresponds to an error of $10^{-23}$ on the final
+% result, small enough in all cases.}
+% \begin{quote}
+% \cs{@@_ln_div_after:Nw} \meta{fixed tl}
+% \cs{@@_div_mantissa_pack:NNN} $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^6 + Q_{4}$
+% \cs{@@_div_mantissa_pack:NNN} $10^6 + Q_{5}$
+% \cs{@@_div_mantissa_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_mantissa_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_use:N \__int_eval:w 9999 - #2 \exp_after:wN ;
+ \int_use:N \__int_eval:w 9999 - #3 \exp_after:wN ;
+ \int_use:N \__int_eval:w 9999 - #4 \exp_after:wN ;
+ \int_use:N \__int_eval:w 9999 - #5 \exp_after:wN ;
+ \int_use:N \__int_eval:w 9999 - #6 \exp_after:wN ;
+ \int_use:N \__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_use:N \__int_eval:w 9999 0000 + #3*#3
+ \exp_after:wN \@@_ln_square_t_pack:NNNNNw
+ \int_use:N \__int_eval:w 9999 0000 + 2*#3*#4
+ \exp_after:wN \@@_ln_square_t_pack:NNNNNw
+ \int_use:N \__int_eval:w 9999 0000 + 2*#3*#5 + #4*#4
+ \exp_after:wN \@@_ln_square_t_pack:NNNNNw
+ \int_use:N \__int_eval:w 9999 0000 + 2*#3*#6 + 2*#4*#5
+ \exp_after:wN \@@_ln_square_t_pack:NNNNNw
+ \int_use:N \__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_use:N \__int_eval:w -1 + 2*#3
+ \exp_after:wN \@@_ln_twice_t_pack:Nw
+ \int_use:N \__int_eval:w 9999 + 2*#4
+ \exp_after:wN \@@_ln_twice_t_pack:Nw
+ \int_use:N \__int_eval:w 9999 + 2*#5
+ \exp_after:wN \@@_ln_twice_t_pack:Nw
+ \int_use:N \__int_eval:w 9999 + 2*#6
+ \exp_after:wN \@@_ln_twice_t_pack:Nw
+ \int_use:N \__int_eval:w 9999 + 2*#7
+ \exp_after:wN \@@_ln_twice_t_pack:Nw
+ \int_use:N \__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:NwNn} \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}
+% \footnote{Bruno: add explanations.}
+%
+% 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 = \c_one
+ \@@_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_use:N \__int_eval:w #1 - \c_two ;
+ }
+ #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_use:N \__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$.\footnote{Bruno: that was wrong at some point, I must check.}
+% \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}
+% \footnote{Bruno: this \emph{\textbf{must}} be updated with correct values!}
+% \end{macro}
+%
+% \begin{macro}{\@@_ln_exponent:wn}
+% \begin{quote}
+% \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:
+ \c_zero \@@_case_return:nw { \@@_fixed_to_float:Nw 2 }
+ \or:
+ \exp_after:wN \@@_ln_exponent_one:ww \__int_value:w
+ \else:
+ \if_int_compare:w #2 > \c_zero
+ \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;
+ {
+ \c_zero
+ \exp_after:wN \@@_fixed_sub:wwN \c_@@_ln_x_fixed_tl ; #1;
+ \@@_fixed_to_float: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;
+ {
+ \c_four
+ \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:wN #1
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsection{Exponential}
+%
+% \subsubsection{Sign, exponent, and special numbers}
+%
+% \begin{macro}[int, EXP]{\@@_exp:w}
+% \begin{macrocode}
+\cs_new:Npn \@@_exp:w \s_@@ \@@_chk:w #1#2
+ {
+ \if_case:w #1 \exp_stop_f:
+ \@@_case_return_o:Nw \c_one_fp
+ \or:
+ \exp_after:wN \@@_exp_normal:w
+ \or:
+ \if_meaning:w 0 #2
+ \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 #1#2
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[aux, EXP]{\@@_exp_normal:w, \@@_exp_pos:Nnwnw}
+% \begin{macrocode}
+\cs_new:Npn \@@_exp_normal:w \s_@@ \@@_chk:w 1#1
+ {
+ \if_meaning:w 0 #1
+ \@@_exp_pos:NNwnw + \@@_fixed_to_float:wN
+ \else:
+ \@@_exp_pos:NNwnw - \@@_fixed_inv_to_float:wN
+ \fi:
+ }
+\cs_new:Npn \@@_exp_pos:NNwnw #1#2#3 \fi: #4#5;
+ {
+ \fi:
+ \exp_after:wN \@@_sanitize:Nw
+ \exp_after:wN 0
+ \__int_value:w #1 \__int_eval:w
+ \if_int_compare:w #4 < - \c_eight
+ \c_one
+ \exp_after:wN \@@_add_big_i:wNww
+ \int_use:N \__int_eval:w \c_one - #4 ;
+ 0 {1000}{0000}{0000}{0000} ; #5;
+ \tex_romannumeral:D
+ \else:
+ \if_int_compare:w #4 > \c_five % cf \c_@@_max_exponent_int
+ \exp_after:wN \@@_exp_overflow:
+ \tex_romannumeral:D
+ \else:
+ \if_int_compare:w #4 < \c_zero
+ \exp_after:wN \use_i:nn
+ \else:
+ \exp_after:wN \use_ii:nn
+ \fi:
+ {
+ \c_zero
+ \@@_decimate:nNnnnn { - #4 }
+ \@@_exp_Taylor:Nnnwn
+ }
+ {
+ \@@_decimate:nNnnnn { \c_sixteen - #4 }
+ \@@_exp_pos_large:NnnNwn
+ }
+ #5
+ {#4}
+ #2 0
+ \tex_romannumeral:D
+ \fi:
+ \fi:
+ \exp_after:wN \c_zero
+ }
+\cs_new:Npn \@@_exp_overflow:
+ { + \c_two * \c_@@_max_exponent_int ; {1000} {0000} {0000} {0000} ; }
+% \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})$. Our only task is to compute 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
+ {
+ \@@_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 = \c_one
+ \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_use:N \__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 lose 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_\tex_romannumeral:D #6:wN \exp_after:wN \cs_end:
+ \exp_after:wN \c_@@_one_fixed_tl
+ \exp_after:wN ;
+ \__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;
+ {
+ \@@_exp_Taylor:Nnnwn ? { } { } 0 #2; {}
+ \@@_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}
+% \begin{tabular}{>{$}c<{$}|*8{>{$}l<{$}}}
+% a^b &-\infty&-y &-n &\pm 0&+n &+y &+\infty&\nan \\
+% \hline
+% +\infty&+0 &+0 &+0 &+1&+\infty &+\infty&+\infty&\nan \\
+% 1<x &+0 &+x^{-y}&+x^{-n} &+1&+x^{n} &+x^{y} &+\infty&\nan \\
+% +1 &+1 &+1 &+1 &+1&+1 &+1 &+1 &+1 \\
+% 0<x<1 &+\infty&+x^{-y}&+x^{-n} &+1&+x^{n} &+x^{y} &+0 &\nan \\
+% +0 &+\infty&+\infty&+\infty &+1&+0 &+0 &+0 &\nan \\
+% -0 &\nan &\nan &\pm\infty &+1&\pm 0 &+0 &+0 &\nan \\
+% -1<-x<0&\nan &\nan &\pm x^{-n}&+1&\pm x^{n}&\nan &+0 &\nan \\
+% -1 &\nan &\nan &\pm 1 &+1&\pm 1 &\nan &\nan &\nan \\
+% -x<-1 &+0 &\nan &\pm x^{-n}&+1&\pm x^{n}&\nan &\nan &\nan \\
+% -\infty&+0 &+0 &\pm 0 &+1&\pm\infty&\nan &\nan &\nan \\
+% \nan &\nan &\nan &\nan &+1&\nan &\nan &\nan &\nan \\
+% \end{tabular}
+% \end{center}
+% 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 a 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: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: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
+ \tex_romannumeral:D -`0 \exp_after:wN \use:nn
+ \fi:
+ {
+ \if_meaning:w 1 #1
+ \exp_after:wN \@@_pow_normal: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. We
+% can thus know the result by comparing the type of $a$ with the sign
+% of $b$, since those conveniently take the same possible values, $0$
+% and~$2$.
+% \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
+ \else:
+ \@@_case_return_o:Nw \c_inf_fp
+ \fi:
+ \s_@@ \@@_chk:w #3#4
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[aux, EXP]{\@@_pow_normal: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:ww}.
+% \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:ww \s_@@ \@@_chk:w 1 #1#2#3; \s_@@ \@@_chk:w #4#5
+ {
+ \if_int_compare:w \pdftex_strcmp:D { #2 #3 }
+ { 1 {1000} {0000} {0000} {0000} } = \c_zero
+ \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 ~
+ \or:
+ \exp_after:wN \@@_pow_npos:Nww
+ \exp_after:wN #5
+ \or:
+ \if_meaning:w 2 #5 \exp_after:wN \reverse_if:N \fi:
+ \if_int_compare:w #2 > \c_zero
+ \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: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((\log\lvert x\rvert + n
+% \log 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: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: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 > \c_zero 0 \else: 2 \fi:
+ \exp_after:wN \@@_pow_npos_aux:Nnww
+ \exp_after:wN \@@_fixed_to_float:wN
+ \else:
+ -
+ \exp_after:wN \@@_pow_npos_aux:Nnww
+ \exp_after:wN \@@_fixed_inv_to_float:wN
+ \fi:
+ {#3}
+ }
+% \end{macrocode}
+% \end{macro}
+%
+%^^A begin[todo]
+% \begin{macro}[aux, EXP]{\@@_pow_npos_aux:Nnww}
+% 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 $-\log x$.
+% \begin{macrocode}
+\cs_new:Npn \@@_pow_npos_aux:Nnww #1#2#3#4; \s_@@ \@@_chk:w 1#5#6#7;
+ {
+ \__int_eval:w
+ \@@_ln_significand:NNNNnnnN #3#4
+ \@@_pow_exponent:wnN {#2}
+ \@@_fixed_mul:wwn #7 {0000}{0000} ;
+ \@@_pow_B:wwN #6;
+ #1 0 % fixed_to_float:wN
+ }
+\cs_new:Npn \@@_pow_exponent:wnN #1; #2
+ {
+ \if_int_compare:w #2 > \c_zero
+ \exp_after:wN \@@_pow_exponent:Nwnnnnnn % n\log 10 - (-\log x)
+ \exp_after:wN +
+ \else:
+ \exp_after:wN \@@_pow_exponent:Nwnnnnnn % -( |n|\log 10 + (-\log x) )
+ \exp_after:wN -
+ \fi:
+ #2; #1;
+ }
+\cs_new:Npn \@@_pow_exponent:Nwnnnnnn #1#2; #3#4#5#6#7#8;
+ { %^^A todo: use that in ln. %^^A todo: log(1.00...) too inaccurate?
+ \exp_after:wN \@@_fixed_mul_after:wwn
+ \int_use:N \__int_eval:w -5 0000
+ \exp_after:wN \@@_fixed_mul_pack:NNNNNw
+ \int_use:N \__int_eval:w 4 9995 0000 #1#2*23025 - #1 #3
+ \exp_after:wN \@@_fixed_mul_pack:NNNNNw
+ \int_use:N \__int_eval:w 4 9995 0000 #1 #2*8509 - #1 #4
+ \exp_after:wN \@@_fixed_mul_pack:NNNNNw
+ \int_use:N \__int_eval:w 4 9995 0000 #1 #2*2994 - #1 #5
+ \exp_after:wN \@@_fixed_mul_pack:NNNNNw
+ \int_use:N \__int_eval:w 4 9995 0000 #1 #2*0456 - #1 #6
+ \exp_after:wN \@@_fixed_mul_pack:NNNNNw
+ \int_use:N \__int_eval:w 5 0000 0000 #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 < \c_zero
+ \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_use:N \__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;
+ {
+ + \c_two * \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 - \c_one \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_neg:w}
+% 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_neg:w}. 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.6)| 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_neg:w
+ \or:
+ \if_int_compare:w \__int_eval:w #1 / \c_two = \c_one
+ \@@_invalid_operation:Nnww \c_nan_fp { ^ } #3; #4;
+ \tex_romannumeral:D -`0
+ \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_neg:w \@@_exp_after_o: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}[aux, rEXP]
+% {
+% \@@_pow_neg_case:w, \@@_pow_neg_case_aux:nnnnn,
+% \@@_pow_neg_case_aux:NNNNNNNNw
+% }
+% This function expects a floating point number, and \enquote{returns}
+% $-1$ if it is an even integer, $0$ if it is an odd integer, and $1$
+% if it is not an integer. Zeros are even, $\pm\infty$ and
+% \texttt{nan} are non-integers. The sign of normal numbers is
+% irrelevant to parity. If the exponent is greater than sixteen, then
+% the number is even. If the exponent is non-positive, the number
+% cannot be an integer. We also separate the ranges of exponent
+% $[1,8]$ and $[9,16]$. In the former case, check that the last $8$
+% digits are zero (otherwise we don't have an integer). In both
+% cases, consider the appropriate $8$ digits, either |#4#5| or |#2#3|,
+% remove the first few: we are then left with \meta{digit}
+% \meta{digits} |;| which would be the digits surrounding the decimal
+% period. If the \meta{digits} are non-zero, the number is not an
+% integer. Otherwise, check the parity of the \meta{digit} and return
+% \cs{c_zero} or \cs{c_minus_one}.
+% \begin{macrocode}
+\cs_new:Npn \@@_pow_neg_case:w \s_@@ \@@_chk:w #1#2#3;
+ {
+ \if_case:w #1 \exp_stop_f:
+ \c_minus_one
+ \or: \@@_pow_neg_case_aux:nnnnn #3
+ \else: \c_one
+ \fi:
+ }
+\cs_new:Npn \@@_pow_neg_case_aux:nnnnn #1#2#3#4#5
+ {
+ \if_int_compare:w #1 > \c_eight
+ \if_int_compare:w #1 > \c_sixteen
+ \c_minus_one
+ \else:
+ \exp_after:wN \exp_after:wN
+ \exp_after:wN \@@_pow_neg_case_aux:NNNNNNNNw
+ \prg_replicate:nn { \c_sixteen - #1 } { 0 } #4#5 ;
+ \fi:
+ \else:
+ \if_int_compare:w #1 > \c_zero
+ \if_int_compare:w #4#5 = \c_zero
+ \exp_after:wN \exp_after:wN
+ \exp_after:wN \@@_pow_neg_case_aux:NNNNNNNNw
+ \prg_replicate:nn { \c_eight - #1 } { 0 } #2#3 ;
+ \else:
+ \c_one
+ \fi:
+ \else:
+ \c_one
+ \fi:
+ \fi:
+ }
+\cs_new:Npn \@@_pow_neg_case_aux:NNNNNNNNw #1#2#3#4#5#6#7#8#9;
+ {
+ \if_int_compare:w 0 #9 = \c_zero
+ \if_int_odd:w #8 \exp_stop_f:
+ \c_zero
+ \else:
+ \c_minus_one
+ \fi:
+ \else:
+ \c_one
+ \fi:
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macrocode}
+%</initex|package>
+% \end{macrocode}
+%
+% \end{implementation}
+%
+% \PrintChanges
+%
+% \PrintIndex \ No newline at end of file