summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/source/latex-dev/l3kernel/l3fp-expo.dtx
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/source/latex-dev/l3kernel/l3fp-expo.dtx')
-rw-r--r--Master/texmf-dist/source/latex-dev/l3kernel/l3fp-expo.dtx1380
1 files changed, 1380 insertions, 0 deletions
diff --git a/Master/texmf-dist/source/latex-dev/l3kernel/l3fp-expo.dtx b/Master/texmf-dist/source/latex-dev/l3kernel/l3fp-expo.dtx
new file mode 100644
index 00000000000..6fbc9de52c7
--- /dev/null
+++ b/Master/texmf-dist/source/latex-dev/l3kernel/l3fp-expo.dtx
@@ -0,0 +1,1380 @@
+% \iffalse meta-comment
+%
+%% File: l3fp-expo.dtx
+%
+% Copyright (C) 2011-2024 The LaTeX Project
+%
+% It may be distributed and/or modified under the conditions of the
+% LaTeX Project Public License (LPPL), either version 1.3c of this
+% license or (at your option) any later version. The latest version
+% of this license is in the file
+%
+% https://www.latex-project.org/lppl.txt
+%
+% This file is part of the "l3kernel bundle" (The Work in LPPL)
+% and all files in that bundle must be distributed together.
+%
+% -----------------------------------------------------------------------
+%
+% The development version of the bundle can be found at
+%
+% https://github.com/latex3/latex3
+%
+% for those people who are interested.
+%
+%<*driver>
+\documentclass[full,kernel]{l3doc}
+\begin{document}
+ \DocInput{\jobname.dtx}
+\end{document}
+%</driver>
+% \fi
+%
+% \title{^^A
+% The \pkg{l3fp-expo} module\\
+% Floating point exponential-related functions^^A
+% }
+% \author{^^A
+% The \LaTeX{} Project\thanks
+% {^^A
+% E-mail:
+% \href{mailto:latex-team@latex-project.org}
+% {latex-team@latex-project.org}^^A
+% }^^A
+% }
+% \date{Released 2024-04-11}
+%
+% \maketitle
+%
+% \begin{documentation}
+%
+% \end{documentation}
+%
+% \begin{implementation}
+%
+% \section{\pkg{l3fp-expo} implementation}
+%
+% \begin{macrocode}
+%<*package>
+% \end{macrocode}
+%
+% \begin{macrocode}
+%<@@=fp>
+% \end{macrocode}
+%
+% \begin{macro}[EXP]
+% {
+% \@@_parse_word_exp:N ,
+% \@@_parse_word_ln:N ,
+% \@@_parse_word_fact: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 ? }
+\cs_new:Npn \@@_parse_word_fact:N
+ { \@@_parse_unary_function:NNN \@@_fact_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]$ is 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}
+% {
+% \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]{\@@_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}[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 < \c_one_int
+ 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}[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{syntax}
+% \meta{continuation} \Arg{Y_1} \Arg{Y_2} \Arg{Y_3} \Arg{Y_4} \Arg{Y_5} \Arg{Y_6} |;|
+% \end{syntax}
+% 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}[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 to be used is 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 \eTeX{} 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{syntax}
+% \cs{@@_ln_x_iv:wnnnnnnnn} \meta{1 or 2} \meta{8d} |;| \Arg{4d} \Arg{4d} \meta{fixed-tl}
+% \end{syntax}
+% 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{syntax}
+% \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{syntax}
+% 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{syntax}
+% \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{syntax}
+% ^^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{syntax}
+% \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{syntax}
+% 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{syntax}
+% \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{syntax}
+% 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 = \c_one_int
+ \@@_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{syntax}
+% \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{syntax}
+% 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 is 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{syntax}
+% \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{syntax}
+% Compute \meta{exponent} times $\ln(10)$. Apart from the cases where
+% \meta{exponent} is $0$ or $1$, the result is necessarily 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 > \c_zero_int
+ \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}[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}[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 < \c_zero_int
+ \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}[EXP]{\@@_exp_Taylor:Nnnwn}
+% \begin{macro}[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 = \c_one_int
+ \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{variable}{\c_@@_exp_intarray}
+% The integer array has $6\times 9\times 4=216$ items encoding the
+% values of $\exp(j\times 10^i)$ for $j=1,\dots,9$ and $i=-1,\dots,4$.
+% Each value is expressed as $\simeq 10^p \times 0.m_1m_2m_3$ with
+% three $8$-digit blocks $m_1$, $m_2$, $m_3$ and an integer
+% exponent~$p$ (one more than the scientific exponent), and these are
+% stored in the integer array as four items: $p$, $10^8+m_1$,
+% $10^8+m_2$, $10^8+m_3$. The various exponentials are stored in
+% increasing order of $j\times 10^i$.
+%
+% Storing this data in an integer array makes it slightly harder to
+% access (slower, too), but uses $16$ bytes of memory per exponential
+% stored, while storing as tokens used around $40$ tokens; tokens have
+% an especially large footprint in Unicode-aware engines.
+% \begin{macrocode}
+\intarray_const_from_clist:Nn \c_@@_exp_intarray
+ {
+ 1 , 1 1105 1709 , 1 1807 5647 , 1 6248 1171 ,
+ 1 , 1 1221 4027 , 1 5816 0169 , 1 8339 2107 ,
+ 1 , 1 1349 8588 , 1 0757 6003 , 1 1039 8374 ,
+ 1 , 1 1491 8246 , 1 9764 1270 , 1 3178 2485 ,
+ 1 , 1 1648 7212 , 1 7070 0128 , 1 1468 4865 ,
+ 1 , 1 1822 1188 , 1 0039 0508 , 1 9748 7537 ,
+ 1 , 1 2013 7527 , 1 0747 0476 , 1 5216 2455 ,
+ 1 , 1 2225 5409 , 1 2849 2467 , 1 6045 7954 ,
+ 1 , 1 2459 6031 , 1 1115 6949 , 1 6638 0013 ,
+ 1 , 1 2718 2818 , 1 2845 9045 , 1 2353 6029 ,
+ 1 , 1 7389 0560 , 1 9893 0650 , 1 2272 3043 ,
+ 2 , 1 2008 5536 , 1 9231 8766 , 1 7740 9285 ,
+ 2 , 1 5459 8150 , 1 0331 4423 , 1 9078 1103 ,
+ 3 , 1 1484 1315 , 1 9102 5766 , 1 0342 1116 ,
+ 3 , 1 4034 2879 , 1 3492 7351 , 1 2260 8387 ,
+ 4 , 1 1096 6331 , 1 5842 8458 , 1 5992 6372 ,
+ 4 , 1 2980 9579 , 1 8704 1728 , 1 2747 4359 ,
+ 4 , 1 8103 0839 , 1 2757 5384 , 1 0077 1000 ,
+ 5 , 1 2202 6465 , 1 7948 0671 , 1 6516 9579 ,
+ 9 , 1 4851 6519 , 1 5409 7902 , 1 7796 9107 ,
+ 14 , 1 1068 6474 , 1 5815 2446 , 1 2146 9905 ,
+ 18 , 1 2353 8526 , 1 6837 0199 , 1 8540 7900 ,
+ 22 , 1 5184 7055 , 1 2858 7072 , 1 4640 8745 ,
+ 27 , 1 1142 0073 , 1 8981 5684 , 1 2836 6296 ,
+ 31 , 1 2515 4386 , 1 7091 9167 , 1 0062 6578 ,
+ 35 , 1 5540 6223 , 1 8439 3510 , 1 0525 7117 ,
+ 40 , 1 1220 4032 , 1 9431 7840 , 1 8020 0271 ,
+ 44 , 1 2688 1171 , 1 4181 6135 , 1 4484 1263 ,
+ 87 , 1 7225 9737 , 1 6812 5749 , 1 2581 7748 ,
+ 131 , 1 1942 4263 , 1 9524 1255 , 1 9365 8421 ,
+ 174 , 1 5221 4696 , 1 8976 4143 , 1 9505 8876 ,
+ 218 , 1 1403 5922 , 1 1785 2837 , 1 4107 3977 ,
+ 261 , 1 3773 0203 , 1 0092 9939 , 1 8234 0143 ,
+ 305 , 1 1014 2320 , 1 5473 5004 , 1 5094 5533 ,
+ 348 , 1 2726 3745 , 1 7211 2566 , 1 5673 6478 ,
+ 391 , 1 7328 8142 , 1 2230 7421 , 1 7051 8866 ,
+ 435 , 1 1970 0711 , 1 1401 7046 , 1 9938 8888 ,
+ 869 , 1 3881 1801 , 1 9428 4368 , 1 5764 8232 ,
+ 1303 , 1 7646 2009 , 1 8905 4704 , 1 8893 1073 ,
+ 1738 , 1 1506 3559 , 1 7005 0524 , 1 9009 7592 ,
+ 2172 , 1 2967 6283 , 1 8402 3667 , 1 0689 6630 ,
+ 2606 , 1 5846 4389 , 1 5650 2114 , 1 7278 5046 ,
+ 3041 , 1 1151 7900 , 1 5080 6878 , 1 2914 4154 ,
+ 3475 , 1 2269 1083 , 1 0850 6857 , 1 8724 4002 ,
+ 3909 , 1 4470 3047 , 1 3316 5442 , 1 6408 6591 ,
+ 4343 , 1 8806 8182 , 1 2566 2921 , 1 5872 6150 ,
+ 8686 , 1 7756 0047 , 1 2598 6861 , 1 0458 3204 ,
+ 13029 , 1 6830 5723 , 1 7791 4884 , 1 1932 7351 ,
+ 17372 , 1 6015 5609 , 1 3095 3052 , 1 3494 7574 ,
+ 21715 , 1 5297 7951 , 1 6443 0315 , 1 3251 3576 ,
+ 26058 , 1 4665 6719 , 1 0099 3379 , 1 5527 2929 ,
+ 30401 , 1 4108 9724 , 1 3326 3186 , 1 5271 5665 ,
+ 34744 , 1 3618 6973 , 1 3140 0875 , 1 3856 4102 ,
+ 39087 , 1 3186 9209 , 1 6113 3900 , 1 6705 9685 ,
+ }
+% \end{macrocode}
+% \end{variable}
+%
+% \begin{macro}[rEXP]
+% {
+% \@@_exp_pos_large:NnnNwn ,
+% \@@_exp_large_after:wwn ,
+% \@@_exp_large:NwN ,
+% \@@_exp_intarray:w ,
+% \@@_exp_intarray_aux:w ,
+% }
+% 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 \cs{@@_exp_large:NwN}, whose
+% |#1| is the \meta{exponent}, |#2| is the current mantissa, and |#3|
+% is the \meta{digit}. At the end, \cs{@@_exp_large_after:wwn} moves
+% on to the Taylor series, eventually multiplied with the mantissa
+% that we have just computed.
+% \begin{macrocode}
+\cs_new:Npn \@@_exp_pos_large:NnnNwn #1#2#3 #4#5; #6
+ {
+ \exp_after:wN \exp_after:wN \exp_after:wN \@@_exp_large:NwN
+ \exp_after:wN \exp_after:wN \exp_after:wN #6
+ \exp_after:wN \c_@@_one_fixed_tl
+ \int_value:w #3 #4 \exp_stop_f:
+ #5 00000 ;
+ }
+\cs_new:Npn \@@_exp_large:NwN #1#2; #3
+ {
+ \if_case:w #3 ~
+ \exp_after:wN \@@_fixed_continue:wn
+ \else:
+ \exp_after:wN \@@_exp_intarray:w
+ \int_value:w \@@_int_eval:w 36 * #1 + 4 * #3 \exp_after:wN ;
+ \fi:
+ #2;
+ {
+ \if_meaning:w 0 #1
+ \exp_after:wN \@@_exp_large_after:wwn
+ \else:
+ \exp_after:wN \@@_exp_large:NwN
+ \int_value:w \@@_int_eval:w #1 - 1 \exp_after:wN \scan_stop:
+ \fi:
+ }
+ }
+\cs_new:Npn \@@_exp_intarray:w #1 ;
+ {
+ +
+ \__kernel_intarray_item:Nn \c_@@_exp_intarray
+ { \@@_int_eval:w #1 - 3 \scan_stop: }
+ \exp_after:wN \use_i:nnn
+ \exp_after:wN \@@_fixed_mul:wwn
+ \int_value:w 0
+ \exp_after:wN \@@_exp_intarray_aux:w
+ \int_value:w \__kernel_intarray_item:Nn
+ \c_@@_exp_intarray { \@@_int_eval:w #1 - 2 }
+ \exp_after:wN \@@_exp_intarray_aux:w
+ \int_value:w \__kernel_intarray_item:Nn
+ \c_@@_exp_intarray { \@@_int_eval:w #1 - 1 }
+ \exp_after:wN \@@_exp_intarray_aux:w
+ \int_value:w \__kernel_intarray_item:Nn \c_@@_exp_intarray {#1} ; ;
+ }
+\cs_new:Npn \@@_exp_intarray_aux:w 1 #1#2#3#4#5 ; { ; {#1#2#3#4} {#5} }
+\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) &-\text{integer} &\pm 0 &+\text{integer} &(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)^b\infty &+1 &(-1)^b 0 &+0 &+0 &\nan \\
+% (-1,0) &+\infty &\nan &(-1)^b\abs{a}^{b} &+1 &(-1)^b\abs{a}^{b} &\nan &+0 &\nan \\
+% -1 &+1 &\nan &(-1)^b &+1 &(-1)^b &\nan &+1 &\nan \\
+% (-\infty,-1) &+0 &\nan &(-1)^b\abs{a}^{b} &+1 &(-1)^b\abs{a}^{b} &\nan &+\infty &\nan \\
+% -\infty &+0 &+0 &(-1)^b 0 &+1 &(-1)^b\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) integer exponents, as $(-1)^b$ 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}[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}[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}[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:w 0 \@@_str_if_eq:nn { #2 #3 } { 1 {1000} {0000} {0000} {0000} }
+ \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 > \c_zero_int
+ \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}[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 > \c_zero_int 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}[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 > \c_zero_int
+ \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 < \c_zero_int
+ \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:NwN
+ \exp_after:wN 5
+ \c_@@_one_fixed_tl
+ }
+% \end{macrocode}
+% \end{macro}
+%^^A end[todo]
+%
+% \begin{macro}[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)| gives $+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 = \c_one_int
+ \@@_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}[rEXP]
+% {
+% \@@_pow_neg_case:w, \@@_pow_neg_case_aux:nnnnn,
+% \@@_pow_neg_case_aux:Nnnw
+% }
+% 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,
+% $0$~if the number is an odd integer, 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. After \cs{@@_decimate:nNnnnn} the
+% argument |#1| of \cs{@@_pow_neg_case_aux:Nnnw} is a rounding digit,
+% |0|~if and only if the number was an integer, and |#3| is the $8$
+% least significant digits of that integer.
+% \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
+ {
+ \if_int_compare:w #1 > \c_@@_prec_int
+ -1
+ \else:
+ \@@_decimate:nNnnnn { \c_@@_prec_int - #1 }
+ \@@_pow_neg_case_aux:Nnnw
+ {#2} {#3} {#4} {#5}
+ \fi:
+ }
+\cs_new:Npn \@@_pow_neg_case_aux:Nnnw #1#2#3#4 ;
+ {
+ \if_meaning:w 0 #1
+ \if_int_odd:w #3 \exp_stop_f:
+ 0
+ \else:
+ -1
+ \fi:
+ \else:
+ 1
+ \fi:
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \subsection{Factorial}
+%
+% \begin{variable}{\c_@@_fact_max_arg_int}
+% The maximum integer whose factorial fits in the exponent range is
+% $3248$, as $3249!\sim 10^{10000.8}$
+% \begin{macrocode}
+\int_const:Nn \c_@@_fact_max_arg_int { 3248 }
+% \end{macrocode}
+% \end{variable}
+%
+% \begin{macro}[EXP]{\@@_fact_o:w}
+% First detect $\pm 0$ and $+\infty$ and \texttt{nan}. Then note that
+% factorial of anything with a negative sign (except $-0$) is
+% undefined. Then call \cs{@@_small_int:wTF} to get an integer as the
+% argument, and start a loop. This is not the most efficient way of
+% computing the factorial, but it works all right. Of course we work
+% with $24$ digits instead of~$16$. It is easy to check that
+% computing factorials with this precision is enough.
+% \begin{macrocode}
+\cs_new:Npn \@@_fact_o:w #1 \s_@@ \@@_chk:w #2#3#4; @
+ {
+ \if_case:w #2 \exp_stop_f:
+ \@@_case_return_o:Nw \c_one_fp
+ \or:
+ \or:
+ \if_meaning:w 0 #3
+ \exp_after:wN \@@_case_return_same_o:w
+ \fi:
+ \or:
+ \@@_case_return_same_o:w
+ \fi:
+ \if_meaning:w 2 #3
+ \@@_case_use:nw { \@@_invalid_operation_o:fw { fact } }
+ \fi:
+ \@@_fact_pos_o:w
+ \s_@@ \@@_chk:w #2 #3 #4 ;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_fact_pos_o:w, \@@_fact_int_o:w}
+% Then check the input is an integer, and call
+% \cs{@@_facorial_int_o:n} with that \texttt{int} as an argument. If
+% it's too big the factorial overflows. Otherwise call
+% \cs{@@_sanitize:Nw} with a positive sign marker~|0| and an integer
+% expression that will mop up any exponent in the calculation.
+% \begin{macrocode}
+\cs_new:Npn \@@_fact_pos_o:w #1;
+ {
+ \@@_small_int:wTF #1;
+ { \@@_fact_int_o:n }
+ { \@@_invalid_operation_o:fw { fact } #1; }
+ }
+\cs_new:Npn \@@_fact_int_o:n #1
+ {
+ \if_int_compare:w #1 > \c_@@_fact_max_arg_int
+ \@@_case_return:nw
+ {
+ \exp_after:wN \exp_after:wN \exp_after:wN \@@_overflow:w
+ \exp_after:wN \c_inf_fp
+ }
+ \fi:
+ \exp_after:wN \@@_sanitize:Nw
+ \exp_after:wN 0
+ \int_value:w \@@_int_eval:w
+ \@@_fact_loop_o:w #1 . 4 , { 1 } { } { } { } { } { } ;
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macro}[EXP]{\@@_fact_loop_o:w}
+% The loop receives an integer |#1| whose factorial we want to
+% compute, which we progressively decrement, and the result so far as
+% an extended-precision number |#2| in the form
+% \meta{exponent}|,|\meta{mantissa}|;|. The loop goes in steps of two
+% because we compute |#1*#1-1| as an integer expression (it must fit
+% since |#1| is at most $3248$), then multiply with the result so far.
+% We don't need to fill in most of the mantissa with zeros because
+% \cs{@@_ep_mul:wwwwn} first normalizes the extended precision number
+% to avoid loss of precision. When reaching a small enough number
+% simply use a table of factorials less than $10^8$. This limit is
+% chosen because the normalization step cannot deal with larger
+% integers.
+% \begin{macrocode}
+\cs_new:Npn \@@_fact_loop_o:w #1 . #2 ;
+ {
+ \if_int_compare:w #1 < 12 \exp_stop_f:
+ \@@_fact_small_o:w #1
+ \fi:
+ \exp_after:wN \@@_ep_mul:wwwwn
+ \exp_after:wN 4 \exp_after:wN ,
+ \exp_after:wN { \int_value:w \@@_int_eval:w #1 * (#1 - 1) }
+ { } { } { } { } { } ;
+ #2 ;
+ {
+ \exp_after:wN \@@_fact_loop_o:w
+ \int_value:w \@@_int_eval:w #1 - 2 .
+ }
+ }
+\cs_new:Npn \@@_fact_small_o:w #1 \fi: #2 ; #3 ; #4
+ {
+ \fi:
+ \exp_after:wN \@@_ep_mul:wwwwn
+ \exp_after:wN 4 \exp_after:wN ,
+ \exp_after:wN
+ {
+ \int_value:w
+ \if_case:w #1 \exp_stop_f:
+ 1 \or: 1 \or: 2 \or: 6 \or: 24 \or: 120 \or: 720 \or: 5040
+ \or: 40320 \or: 362880 \or: 3628800 \or: 39916800
+ \fi:
+ } { } { } { } { } { } ;
+ #3 ;
+ \@@_ep_to_float_o:wwN 0
+ }
+% \end{macrocode}
+% \end{macro}
+%
+% \begin{macrocode}
+%</package>
+% \end{macrocode}
+%
+% \end{implementation}
+%
+% \PrintChanges
+%
+% \PrintIndex