diff options
Diffstat (limited to 'Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx')
-rw-r--r-- | Master/texmf-dist/source/latex/l3kernel/l3fp-expo.dtx | 1252 |
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 |