% \iffalse meta-comment % %% File: l3fp-extended.dtx Copyright (C) 2011-2012 The LaTeX3 Project %% %% It may be distributed and/or modified under the conditions of the %% LaTeX Project Public License (LPPL), either version 1.3c of this %% license or (at your option) any later version. The latest version %% of this license is in the file %% %% http://www.latex-project.org/lppl.txt %% %% This file is part of the "l3kernel bundle" (The Work in LPPL) %% and all files in that bundle must be distributed together. %% %% The released version of this bundle is available from CTAN. %% %% ----------------------------------------------------------------------- %% %% The development version of the bundle can be found at %% %% http://www.latex-project.org/svnroot/experimental/trunk/ %% %% for those people who are interested. %% %%%%%%%%%%% %% NOTE: %% %%%%%%%%%%% %% %% Snapshots taken from the repository represent work in progress and may %% not work or may contain conflicting material! We therefore ask %% people _not_ to put them into distributions, archives, etc. without %% prior consultation with the LaTeX Project Team. %% %% ----------------------------------------------------------------------- %% % %<*driver> \RequirePackage{l3bootstrap} \GetIdInfo$Id: l3fp-extended.dtx 4339 2012-11-24 19:16:43Z joseph $ {L3 Floating-point extended precision fixed-points} \documentclass[full]{l3doc} \begin{document} \DocInput{\jobname.dtx} \end{document} % % \fi % % \title{The \textsf{l3fp-extended} package\thanks{This file % has version number \ExplFileVersion, last % revised \ExplFileDate.}\\ % Fixed points with extended precision for internal use} % \author{^^A % The \LaTeX3 Project\thanks % {^^A % E-mail: % \href{mailto:latex-team@latex-project.org} % {latex-team@latex-project.org}^^A % }^^A % } % \date{Released \ExplFileDate} % % \maketitle % % \begin{documentation} % % \end{documentation} % % \begin{implementation} % % \section{\pkg{l3fp-extended} implementation} % % \begin{macrocode} %<*initex|package> % \end{macrocode} % % \begin{macrocode} %<@@=fp> % \end{macrocode} % % \subsection{Description of extended fixed points} % % In this module, we work on (almost) fixed-point numbers with % extended ($24$ digits) precision. This is used in the computation of % Taylor series for the logarithm, exponential, and trigonometric % functions. Since we eventually only care about the $16$ first digits % of the final result, some of the calculations are not performed with % the full $24$-digit precision. In other words, the last two blocks of % each fixed point number may be wrong as long as the error is small % enough to be rounded away when converting back to a floating point % number. The fixed point numbers are expressed as % \begin{quote} % \Arg{a_1} \Arg{a_2} \Arg{a_3} \Arg{a_4} \Arg{a_5} \Arg{a_6} |;| % \end{quote} % where each \meta{a_i} is exactly $4$ digits (ranging from |0000| to % |9999|), except \meta{a_1}, which may be any \enquote{not-too-large} % non-negative integer, with or without trailing zeros. Here, % \enquote{not-too-large} depends on the specific function (see the % corresponding comments for details). Checking for overflow is the % responsibility of the code calling those functions. The fixed point % number $a$ corresponding to the representation above is $a = % \sum_{i=1}^{6} \meta{a_i} \cdot 10^{-4i}$. % % Most functions we define here have the form % \begin{syntax} % \cs{@@_fixed_\meta{calculation}:wwn} \meta{operand_1} |;| \meta{operand_2} |;| \Arg{continuation} % \end{syntax} % They perform the \meta{calculation} on the two \meta{operands}, then % feed the result ($6$ brace groups followed by a semicolon) to the % \meta{continuation}, responsible for the next step of the calculation. % Some functions only accept an \texttt{N}-type \meta{continuation}. % This allows constructions such as % \begin{quote} % \cs{@@_fixed_add:wwn} \meta{X_1} |;| \meta{X_2} |;| \\ % \cs{@@_fixed_mul:wwn} \meta{X_3} |;| \\ % \cs{@@_fixed_add:wwn} \meta{X_4} |;| \\ % \end{quote} % to compute $(X_1+X_2)\cdot X_3 + X_4$. This turns out to be very % appropriate for computing continued fractions and Taylor series. % % At the end of the calculation, the result is turned back to a floating % point number using \cs{@@_fixed_to_float:Nw}. This function has to % change the exponent of the floating point number: it must be used % after starting an integer expression for the overall exponent of the % result. % % \subsection{Helpers for extended fixed points} % % \begin{variable}[int]{\c_@@_one_fixed_tl} % The extended fixed-point number~$1$, used in \pkg{l3fp-expo}. % \begin{macrocode} \tl_const:Nn \c_@@_one_fixed_tl { {10000} {0000} {0000} {0000} {0000} {0000} } % \end{macrocode} % \end{variable} % % \begin{macro}[int, EXP]{\@@_fixed_continue:wn} % This function does nothing. Of course, there is no bound on % $a_1$ (except \TeX{}'s own $2^{31}-1$). % \begin{macrocode} \cs_new:Npn \@@_fixed_continue:wn #1; #2 { #2 #1; } % \end{macrocode} % \end{macro} % % \begin{macro}[int, EXP]{\@@_fixed_add_one:wN} % \begin{syntax} % \cs{@@_fixed_add_one:wN} \meta{a} |;| \meta{continuation} % \end{syntax} % This function adds $1$ to the fixed point \meta{a}, by changing % $a_1$ to $10000+a_1$, then calls the \meta{continuation}. This % requires $a_1 \leq 2^{31} - 10001$. % \begin{macrocode} \cs_new:Npn \@@_fixed_add_one:wN #1#2; #3 { \exp_after:wN #3 \exp_after:wN { \int_use:N \__int_eval:w \c_ten_thousand + #1 } #2 ; } % \end{macrocode} % \end{macro} % % \begin{macro}[aux, EXP]{\@@_fixed_mul_after:wn} % The fixed point operations which involve multiplication end by % calling this auxiliary. It braces the last block of digits, and % places the \meta{continuation} |#2| in front. The % \meta{continuation} was brought up through the expansions by % the packing functions. % \begin{macrocode} \cs_new:Npn \@@_fixed_mul_after:wn #1; #2 { #2 {#1} } % \end{macrocode} % \end{macro} % % \subsection{Dividing a fixed point number by a small integer} % % \begin{macro}[int, EXP]{\@@_fixed_div_int:wwN} % \begin{macro}[aux, EXP] % { % \@@_fixed_div_int:wnN, \@@_fixed_div_int_auxi:wnn, % \@@_fixed_div_int_auxii:wnn, \@@_fixed_div_int_pack:Nw, % \@@_fixed_div_int_after:Nw % } % \begin{syntax} % \cs{@@_fixed_div_int:wwN} \meta{a} |;| \meta{n} |;| \meta{continuation} % \end{syntax} % Divides the fixed point number \meta{a} by the (small) integer % $0<\meta{n}<10^4$ and feeds the result to the \meta{continuation}. % There is no bound on $a_1$. % % The arguments of the \texttt{i} auxiliary are 1: one of the $a_{i}$, % 2: $n$, 3: the \texttt{ii} or the \texttt{iii} auxiliary. It % computes a (somewhat tight) lower bound $Q_{i}$ for the ratio % $a_{i}/n$. % % The \texttt{ii} auxiliary receives $Q_{i}$, $n$, and $a_{i}$ as % arguments. It adds $Q_{i}$ to a surrounding integer expression, and % starts a new one with the initial value $9999$, which ensures that % the result of this expression will have $5$ digits. The auxiliary % also computes $a_{i}-n\cdot Q_{i}$, placing the result in front of % the $4$ digits of $a_{i+1}$. The resulting $a'_{i+1} = 10^{4} % (a_{i} - n \cdot Q_{i}) + a_{i+1}$ serves as the first argument for % a new call to the \texttt{i} auxiliary. % % When the \texttt{iii} auxiliary is called, the situation looks like % this: % \begin{quote} % \cs{@@_fixed_div_int_after:Nw} \meta{continuation} \\ % $-1 + Q_{1}$ \\ % \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{2}$ \\ % \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{3}$ \\ % \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{4}$ \\ % \cs{@@_fixed_div_int_pack:Nw} $9999 + Q_{5}$ \\ % \cs{@@_fixed_div_int_pack:Nw} $9999$ \\ % \cs{@@_fixed_div_int_auxii:wnn} $Q_{6}$ |;| \Arg{n} \Arg{a_{6}} % \end{quote} % where expansion is happening from the last line up. The % \texttt{iii} auxiliary adds $Q_{6} + 2 \simeq a_{6}/n + 1$ to the % last $9999$, giving the integer closest to $10000 + a_{6}/n$. % % Each \texttt{pack} auxiliary receives $5$ digits followed by a % semicolon. The first digit is added as a carry to the integer % expression above, and the $4$ other digits are braced. Each call to % the \texttt{pack} auxiliary thus produces one brace group. The last % brace group is produced by the \texttt{after} auxiliary, which % places the \meta{continuation} as appropriate. % \begin{macrocode} \cs_new:Npn \@@_fixed_div_int:wwN #1#2#3#4#5#6 ; #7 ; #8 { \exp_after:wN \@@_fixed_div_int_after:Nw \exp_after:wN #8 \int_use:N \__int_eval:w \c_minus_one \@@_fixed_div_int:wnN #1; {#7} \@@_fixed_div_int_auxi:wnn #2; {#7} \@@_fixed_div_int_auxi:wnn #3; {#7} \@@_fixed_div_int_auxi:wnn #4; {#7} \@@_fixed_div_int_auxi:wnn #5; {#7} \@@_fixed_div_int_auxi:wnn #6; {#7} \@@_fixed_div_int_auxii:wnn ; } \cs_new:Npn \@@_fixed_div_int:wnN #1; #2 #3 { \exp_after:wN #3 \int_use:N \__int_eval:w #1 / #2 - \c_one ; {#2} {#1} } \cs_new:Npn \@@_fixed_div_int_auxi:wnn #1; #2 #3 { + #1 \exp_after:wN \@@_fixed_div_int_pack:Nw \int_use:N \__int_eval:w 9999 \exp_after:wN \@@_fixed_div_int:wnN \int_use:N \__int_eval:w #3 - #1*#2 \__int_eval_end: } \cs_new:Npn \@@_fixed_div_int_auxii:wnn #1; #2 #3 { + #1 + \c_two ; } \cs_new:Npn \@@_fixed_div_int_pack:Nw #1 #2; { + #1; {#2} } \cs_new:Npn \@@_fixed_div_int_after:Nw #1 #2; { #1 {#2} } % \end{macrocode} % \end{macro} % \end{macro} % % \subsection{Adding and subtracting fixed points} % % \begin{macro}[int, EXP]{\@@_fixed_add:wwn, \@@_fixed_sub:wwn} % \begin{macro}[aux, EXP] % { % \@@_fixed_add:Nnnnnwnn, % \@@_fixed_add:nnNnnnwn, % \@@_fixed_add_pack:NNNNNwn, % \@@_fixed_add_after:NNNNNwn % } % \begin{syntax} % \cs{@@_fixed_add:wwn} \meta{a} |;| \meta{b} |;| \Arg{continuation} % \end{syntax} % Computes $a+b$ (resp.\ $a-b$) and feeds the result to the % \meta{continuation}. This function requires $0\leq % a_{1},b_{1}<50000$, and requires the result to be positive (this % happens automatically for addition). The two functions only differ % a sign, hence use a common auxiliary. It would be nice to grab the % $12$ brace groups in one go; only $9$ parameters are allowed. Start % by grabbing the two signs, $a_{1}, \ldots, a_{4}$, the rest of $a$, % and $b_{1}$ and $b_{2}$. The second auxiliary receives the rest of % $a$, the sign multiplying $b$, the rest of $b$, and the % \meta{continuation} as arguments. After going down through the % various level, we go back up, packing digits and bringing the % \meta{continuation} (|#8|, then |#7|) from the end of the argument % list to its start. % \begin{macrocode} \cs_new_nopar:Npn \@@_fixed_add:wwn { \@@_fixed_add:Nnnnnwnn + } \cs_new_nopar:Npn \@@_fixed_sub:wwn { \@@_fixed_add:Nnnnnwnn - } \cs_new:Npn \@@_fixed_add:Nnnnnwnn #1 #2#3#4#5 #6; #7#8 { \exp_after:wN \@@_fixed_add_after:NNNNNwn \int_use:N \__int_eval:w 9 9999 9998 + #2#3 #1 #7#8 \exp_after:wN \@@_fixed_add_pack:NNNNNwn \int_use:N \__int_eval:w 1 9999 9998 + #4#5 \@@_fixed_add:nnNnnnwn #6 #1 } \cs_new:Npn \@@_fixed_add:nnNnnnwn #1#2 #3 #4#5 #6#7 ; #8 { #3 #4#5 \exp_after:wN \@@_fixed_add_pack:NNNNNwn \int_use:N \__int_eval:w 2 0000 0000 #3 #6#7 + #1#2 ; {#8} ; } \cs_new:Npn \@@_fixed_add_pack:NNNNNwn #1 #2#3#4#5 #6; #7 { + #1 ; {#7} {#2#3#4#5} {#6} } \cs_new:Npn \@@_fixed_add_after:NNNNNwn 1 #1 #2#3#4#5 #6; #7 { #7 {#1#2#3#4#5} {#6} } % \end{macrocode} % \end{macro} % \end{macro} % % \subsection{Multiplying fixed points} % % \begin{macro}[int, EXP]{\@@_fixed_mul:wwn} % \begin{macro}[aux, EXP]{\@@_fixed_mul:nnnnnnnwn} % \begin{syntax} % \cs{@@_fixed_mul:wwn} \meta{a} |;| \meta{b} |;| \Arg{continuation} % \end{syntax} % Computes $a\times b$ and feeds the result to \meta{continuation}. % This function requires $0\leq a_{1}, b_{1} < 10000$. Once more, we % need to play around the limit of $9$ arguments for \TeX{} macros. % Note that we don't need to obtain an exact rounding, contrarily to % the |*| operator, so things could be harder. We wish to perform % carries in % \begin{align*} % a \times b = % & a_{1} \cdot b_{1} \cdot 10^{-8} \\ % & + (a_{1} \cdot b_{2} + a_{2} \cdot b_{1}) \cdot 10^{-12} \\ % & + (a_{1} \cdot b_{3} + a_{2} \cdot b_{2} % + a_{3} \cdot b_{1}) \cdot 10^{-16} \\ % & + (a_{1} \cdot b_{4} + a_{2} \cdot b_{3} % + a_{3} \cdot b_{2} + a_{4} \cdot b_{1}) \cdot 10^{-20} \\ % & + \left(a_{2} \cdot b_{4} + a_{3} \cdot b_{3} + a_{4} \cdot b_{2} % + \frac{a_{3} \cdot b_{4} + a_{4} \cdot b_{3} % + a_{1} \cdot b_{6} + a_{2} \cdot b_{5} % + a_{5} \cdot b_{2} + a_{6} \cdot b_{1}}{10^{4}} % + a_{1} \cdot b_{5} + a_{5} \cdot b_{1}\right) \cdot 10^{-24} % + O(10^{-24}), % \end{align*} % where the $O(10^{-24})$ stands for terms which are at most $5\cdot % 10^{-24}$; ignoring those leads to an error of at most % $5$~\texttt{ulp}. Note how the first $15$~terms only depend on % $a_{1},\ldots{},a_{4}$ and $b_{1},\ldots,b_{4}$, while the last % $6$~terms only depend on $a_{1},a_{2},a_{5},a_{6}$, and the % corresponding parts of~$b$. Hence, the first function grabs % $a_{1},\ldots,a_{4}$, the rest of $a$, and $b_{1},\ldots,b_{4}$, and % writes the $15$ first terms of the expression, including a left % parenthesis for the fraction. The \texttt{i} auxiliary receives % $a_{5}$, $a_{6}$, $b_{1}$, $b_{2}$, $a_{1}$, $a_{2}$, $b_{5}$, % $b_{6}$ and finally the \meta{continuation} as arguments. It writes % the end of the expression, including the right parenthesis and the % denominator of the fraction. The packing auxiliaries bring the % \meta{continuation} up through the expansion chain, as |#7|, and it % is finally placed in front of the $6$ brace groups by % \cs{@@_fixed_mul_after:wn}. % \begin{macrocode} \cs_new:Npn \@@_fixed_mul:wwn #1#2#3#4 #5; #6#7#8#9 { \exp_after:wN \@@_fixed_mul_after:wn \int_use:N \__int_eval:w \c_@@_leading_shift_int \exp_after:wN \@@_pack:NNNNNwn \int_use:N \__int_eval:w \c_@@_middle_shift_int + #1*#6 \exp_after:wN \@@_pack:NNNNNwn \int_use:N \__int_eval:w \c_@@_middle_shift_int + #1*#7 + #2*#6 \exp_after:wN \@@_pack:NNNNNwn \int_use:N \__int_eval:w \c_@@_middle_shift_int + #1*#8 + #2*#7 + #3*#6 \exp_after:wN \@@_pack:NNNNNwn \int_use:N \__int_eval:w \c_@@_middle_shift_int + #1*#9 + #2*#8 + #3*#7 + #4*#6 \exp_after:wN \@@_pack:NNNNNwn \int_use:N \__int_eval:w \c_@@_trailing_shift_int + #2*#9 + #3*#8 + #4*#7 + ( #3*#9 + #4*#8 + \@@_fixed_mul:nnnnnnnwn #5 {#6}{#7} {#1}{#2} } \cs_new:Npn \@@_fixed_mul:nnnnnnnwn #1#2 #3#4 #5#6 #7#8 ; #9 { #1*#4 + #2*#3 + #5*#8 + #6*#7 ) / \c_ten_thousand + #1*#3 + #5*#7 ; {#9} ; } % \end{macrocode} % \end{macro} % \end{macro} % % \subsection{Combining product and sum of fixed points} % % \begin{macro}[int, EXP] % { % \@@_fixed_mul_add:wwwn, % \@@_fixed_mul_sub_back:wwwn, % \@@_fixed_mul_one_minus_mul:wwn, % } % \begin{syntax} % \cs{@@_fixed_mul_add:wwwn} \meta{a} |;| \meta{b} |;| \meta{c} |;| \Arg{continuation} % \cs{@@_fixed_mul_sub_back:wwwn} \meta{a} |;| \meta{b} |;| \meta{c} |;| \Arg{continuation} % \cs{@@_fixed_one_minus_mul:wwn} \meta{a} |;| \meta{b} |;| \Arg{continuation} % \end{syntax} % Compute $a\times b + c$, $c - a\times b$, and $1 - a\times b$ and % feed the result to the \meta{continuation}. Those functions require % $0\leq a_{1}, b_{1}, c_{1} \leq 10000$. Since those functions are % at the heart of the computation of Taylor expansions, we % over-optimize them a bit, and in particular we do not factor out the % common parts of the three functions. % % For definiteness, consider the task of computing $a\times b + c$. % We will perform carries in % \begin{align*} % a \times b + c = % & (a_{1} \cdot b_{1} + c_{1} c_{2})\cdot 10^{-8} \\ % & + (a_{1} \cdot b_{2} + a_{2} \cdot b_{1}) \cdot 10^{-12} \\ % & + (a_{1} \cdot b_{3} + a_{2} \cdot b_{2} + a_{3} \cdot b_{1} % + c_{3} c_{4}) \cdot 10^{-16} \\ % & + (a_{1} \cdot b_{4} + a_{2} \cdot b_{3} + a_{3} \cdot b_{2} % + a_{4} \cdot b_{1}) \cdot 10^{-20} \\ % & + \Big(a_{2} \cdot b_{4} + a_{3} \cdot b_{3} + a_{4} \cdot b_{2} % + \frac{a_{3} \cdot b_{4} + a_{4} \cdot b_{3} % + a_{1} \cdot b_{6} + a_{2} \cdot b_{5} % + a_{5} \cdot b_{2} + a_{6} \cdot b_{1}}{10^{4}} % + a_{1} \cdot b_{5} + a_{5} \cdot b_{1} % + c_{5} c_{6} \Big) \cdot 10^{-24} % + O(10^{-24}), % \end{align*} % where $c_{1} c_{2}$, $c_{3} c_{4}$, $c_{5} c_{6}$ denote the % $8$-digit number obtained by juxtaposing the two blocks of digits of % $c$, and $\cdot$ denotes multiplication. The task is obviously % tough because we have $18$ brace groups in front of us. % % Each of the three function starts the first two levels (the first, % corresponding to $10^{-4}$, is empty), with $c_{1} c_{2}$ in the % first level, calls the \texttt{i} auxiliary with arguments described % later, and adds a trailing ${} + c_{5}c_{6}$ |;| % \Arg{continuation}~|;|. The ${} + c_{5}c_{6}$ piece, which is % omitted for \cs{@@_fixed_one_minus_mul:wwn}, will be taken in the % integer expression for the $10^{-24}$ level. The % \meta{continuation} is placed correctly to be taken upstream by % packing auxiliaries. % \begin{macrocode} \cs_new:Npn \@@_fixed_mul_add:wwwn #1; #2; #3#4#5#6#7#8; #9 { \exp_after:wN \@@_fixed_mul_after:wn \int_use:N \__int_eval:w \c_@@_big_leading_shift_int \exp_after:wN \@@_pack_big:NNNNNNwn \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #3 #4 \@@_fixed_mul_add:Nwnnnwnnn + + #5 #6 ; #2 ; #1 ; #2 ; + + #7 #8 ; {#9} ; } \cs_new:Npn \@@_fixed_mul_sub_back:wwwn #1; #2; #3#4#5#6#7#8; #9 { \exp_after:wN \@@_fixed_mul_after:wn \int_use:N \__int_eval:w \c_@@_big_leading_shift_int \exp_after:wN \@@_pack_big:NNNNNNwn \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #3 #4 \@@_fixed_mul_add:Nwnnnwnnn - + #5 #6 ; #2 ; #1 ; #2 ; - + #7 #8 ; {#9} ; } \cs_new:Npn \@@_fixed_one_minus_mul:wwn #1; #2; #3 { \exp_after:wN \@@_fixed_mul_after:wn \int_use:N \__int_eval:w \c_@@_big_leading_shift_int \exp_after:wN \@@_pack_big:NNNNNNwn \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + 1 0000 0000 \@@_fixed_mul_add:Nwnnnwnnn - ; #2 ; #1 ; #2 ; - ; {#3} ; } % \end{macrocode} % \end{macro} % % \begin{macro}[aux, EXP]{\@@_fixed_mul_add:Nwnnnwnnn} % \begin{syntax} % \cs{@@_fixed_mul_add:Nwnnnwnnn} \meta{op} |+| \meta{c_3} \meta{c_4} |;| % ~~\meta{b} |;| \meta{a} |;| \meta{b} |;| \meta{op} % ~~|+| \meta{c_5} \meta{c_6} |;| % \end{syntax} % Here, \meta{op} is either |+| or |-|. Arguments |#3|, |#4|, |#5| % are \meta{b_1}, \meta{b_2}, \meta{b_3}; arguments |#7|, |#8|, |#9| % are \meta{a_1}, \meta{a_2}, \meta{a_3}. We can build three levels: % $a_{1} \cdot b_{1}$ for $10^{-8}$, $(a_{1} \cdot b_{2} + a_{2} \cdot % b_{1})$ for $10^{-12}$, and $(a_{1} \cdot b_{3} + a_{2} \cdot b_{2} % + a_{3} \cdot b_{1} + c_{3} c_{4})$ for $10^{-16}$. The $a$--$b$ % products huse the sign |#1|. Note that |#2| is empty for % \cs{@@_fixed_one_minus_mul:wwn}. We call the \texttt{ii} auxiliary % for levels $10^{-20}$ and $10^{-24}$, keeping the pieces of \meta{a} % we've read, but not \meta{b}, since there is another copy later in % the input stream. % \begin{macrocode} \cs_new:Npn \@@_fixed_mul_add:Nwnnnwnnn #1 #2; #3#4#5#6; #7#8#9 { #1 #7*#3 \exp_after:wN \@@_pack_big:NNNNNNwn \int_use:N \__int_eval:w \c_@@_big_middle_shift_int #1 #7*#4 #1 #8*#3 \exp_after:wN \@@_pack_big:NNNNNNwn \int_use:N \__int_eval:w \c_@@_big_middle_shift_int #1 #7*#5 #1 #8*#4 #1 #9*#3 #2 \exp_after:wN \@@_pack_big:NNNNNNwn \int_use:N \__int_eval:w \c_@@_big_middle_shift_int #1 \@@_fixed_mul_add:nnnnwnnnn {#7}{#8}{#9} } % \end{macrocode} % \end{macro} % % \begin{macro}[aux, EXP]{\@@_fixed_mul_add:nnnnwnnnn} % \begin{syntax} % \cs{@@_fixed_mul_add:nnnnwnnnn} \meta{a} |;| \meta{b} |;| \meta{op} % ~~|+| \meta{c_5} \meta{c_6} |;| % \end{syntax} % Level $10^{-20}$ is $(a_{1} \cdot b_{4} + a_{2} \cdot b_{3} + a_{3} % \cdot b_{2} + a_{4} \cdot b_{1})$, multiplied by the sign, which was % inserted by the \texttt{i} auxiliary. Then we prepare level % $10^{-24}$. We don't have access to all parts of \meta{a} and % \meta{b} needed to make all products. Instead, we prepare the % partial expressions % \begin{align*} % & b_{1} + a_{4} \cdot b_{2} + a_{3} \cdot b_{3} + a_{2} \cdot b_{4} + a_{1} \\ % & b_{2} + a_{4} \cdot b_{3} + a_{3} \cdot b_{4} + a_{2} . % \end{align*} % Obviously, those expressions make no mathematical sense: we will % complete them with $a_{5} \cdot {}$ and ${} \cdot b_{5}$, and with % $a_{6} \cdot b_{1} + a_{5} \cdot {}$ and ${} \cdot b_{5} + a_{1} % \cdot b_{6}$, and of course with the trailing ${} + c_{5} c_{6}$. % To do all this, we keep $a_{1}$, $a_{5}$, $a_{6}$, and the % corresponding pieces of \meta{b}. % \begin{macrocode} \cs_new:Npn \@@_fixed_mul_add:nnnnwnnnn #1#2#3#4#5; #6#7#8#9 { ( #1*#9 + #2*#8 + #3*#7 + #4*#6 ) \exp_after:wN \@@_pack_big:NNNNNNwn \int_use:N \__int_eval:w \c_@@_big_trailing_shift_int \@@_fixed_mul_add:nnnnwnnwN { #6 + #4*#7 + #3*#8 + #2*#9 + #1 } { #7 + #4*#8 + #3*#9 + #2 } {#1} #5; {#6} } % \end{macrocode} % \end{macro} % % \begin{macro}[aux, EXP]{\@@_fixed_mul_add:nnnnwnnwN} % \begin{syntax} % \cs{@@_fixed_mul_add:nnnnwnnwN} \Arg{partial_1} \Arg{partial_2} % ~~\Arg{a_1} \Arg{a_5} \Arg{a_6} |;| \Arg{b_1} \Arg{b_5} \Arg{b_6} |;| % ~~\meta{op} |+| \meta{c_5} \meta{c_6} |;| % \end{syntax} % Complete the \meta{partial_1} and \meta{partial_2} expressions as % explained for the \texttt{ii} auxiliary. The second one is divided % by $10000$: this is the carry from level $10^{-28}$. The trailing % ${} + c_{5} c_{6}$ is taken into the expression for level % $10^{-24}$. Note that the total of level $10^{-24}$ is in the % interval $[-5\cdot 10^{8}, 6\cdot 10^{8}$ (give or take a couple of % $10000$), hence adding it to the shift gives a $10$-digit number, as % expected by the packing auxiliaries. See \pkg{l3fp-aux} for the % definition of the shifts and packing auxiliaries. % \begin{macrocode} \cs_new:Npn \@@_fixed_mul_add:nnnnwnnwN #1#2 #3#4#5; #6#7#8; #9 { #9 (#4* #1 *#7) #9 (#5*#6+#4* #2 *#7+#3*#8) / \c_ten_thousand } % \end{macrocode} % \end{macro} % % \subsection{Converting from fixed point to floating point} % % \begin{macro}[int, rEXP] % {\@@_fixed_to_float:wN, \@@_fixed_to_float:Nw} % \begin{syntax} % \ldots{} \cs{__int_eval:w} \meta{exponent} \cs{@@_fixed_to_float:wN} \Arg{a_1} \Arg{a_2} \Arg{a_3} \Arg{a_4} \Arg{a_5} \Arg{a_6} |;| \meta{sign} % \end{syntax} % yields % \begin{quote} % \meta{exponent'} |;| \Arg{a'_1} \Arg{a'_2} \Arg{a'_3} \Arg{a'_4} |;| % \end{quote} % And the \texttt{to_fixed} version gives six brace groups instead of % $4$, ensuring that $1000\leq\meta{a'_1}\leq 9999$. At this stage, we % know that \meta{a_1} is positive (otherwise, it is sign of an error % before), and we assume that it is less than $10^8$.\footnote{Bruno: % I must double check this assumption.} % %^^A todo: round properly when rounding to infinity: I need to know the sign. % \begin{macrocode} \cs_new:Npn \@@_fixed_to_float:Nw #1#2; { \@@_fixed_to_float:wN #2; #1 } \cs_new:Npn \@@_fixed_to_float:wN #1#2#3#4#5#6; #7 { + \c_four % for the 8-digit-at-the-start thing. \exp_after:wN \exp_after:wN \exp_after:wN \@@_fixed_to_loop:N \exp_after:wN \use_none:n \int_use:N \__int_eval:w 1 0000 0000 + #1 \exp_after:wN \@@_use_none_stop_f:n \__int_value:w 1#2 \exp_after:wN \@@_use_none_stop_f:n \__int_value:w 1#3#4 \exp_after:wN \@@_use_none_stop_f:n \__int_value:w 1#5#6 \exp_after:wN ; \exp_after:wN ; } \cs_new:Npn \@@_fixed_to_loop:N #1 { \if_meaning:w 0 #1 - \c_one \exp_after:wN \@@_fixed_to_loop:N \else: \exp_after:wN \@@_fixed_to_loop_end:w \exp_after:wN #1 \fi: } \cs_new:Npn \@@_fixed_to_loop_end:w #1 #2 ; { \if_meaning:w ; #1 \exp_after:wN \@@_fixed_to_float_zero:w \else: \exp_after:wN \@@_pack_twice_four:wNNNNNNNN \exp_after:wN \@@_pack_twice_four:wNNNNNNNN \exp_after:wN \@@_fixed_to_float_pack:ww \exp_after:wN ; \fi: #1 #2 0000 0000 0000 0000 ; } \cs_new:Npn \@@_fixed_to_float_zero:w ; 0000 0000 0000 0000 ; { - \c_two * \c_@@_max_exponent_int ; {0000} {0000} {0000} {0000} ; } \cs_new:Npn \@@_fixed_to_float_pack:ww #1 ; #2#3 ; ; { \if_int_compare:w #2 > \c_four \exp_after:wN \@@_fixed_to_float_round_up:wnnnnw \fi: ; #1 ; } \cs_new:Npn \@@_fixed_to_float_round_up:wnnnnw ; #1#2#3#4 ; { \exp_after:wN \@@_basics_pack_high:NNNNNw \int_use:N \__int_eval:w 1 #1#2 \exp_after:wN \@@_basics_pack_low:NNNNNw \int_use:N \__int_eval:w 1 #3#4 + \c_one ; } % \end{macrocode} % \end{macro} % % \begin{macro}[rEXP, int]{\@@_fixed_inv_to_float:wN, \@@_fixed_div_to_float:ww} % Starting from \texttt{fixed_dtf} $A$ |;| $B$ |;| we want to compute % $A/B$, and express it as a floating point number. Normalize both % numbers by removing leading brace groups of zeros and leaving the % appropriate exponent shift in the input stream. % \begin{macrocode} \cs_new:Npn \@@_fixed_inv_to_float:wN #1#2; #3 { + \__int_eval:w % ^^A todo: remove the +? \if_int_compare:w #1 < \c_one_thousand \@@_fixed_dtf_zeros:wNnnnnnn \fi: \@@_fixed_dtf_no_zero:Nwn + {#1} #2 \s_@@ \@@_fixed_dtf_approx:n {10000} {0000} {0000} {0000} {0000} {0000} ; } \cs_new:Npn \@@_fixed_div_to_float:ww #1#2; #3#4; { \if_int_compare:w #1 < \c_one_thousand \@@_fixed_dtf_zeros:wNnnnnnn \fi: \@@_fixed_dtf_no_zero:Nwn - {#1} #2 \s_@@ { \if_int_compare:w #3 < \c_one_thousand \@@_fixed_dtf_zeros:wNnnnnnn \fi: \@@_fixed_dtf_no_zero:Nwn + {#3} #4 \s_@@ \@@_fixed_dtf_approx:n } } \cs_new:Npn \@@_fixed_dtf_no_zero:Nwn #1#2 \s_@@ #3 { #3 #2; } \cs_new:Npn \@@_fixed_dtf_zeros:wNnnnnnn \fi: \@@_fixed_dtf_no_zero:Nwn #1#2#3#4#5#6#7 { \fi: #1 \c_minus_one \exp_after:wN \use_i_ii:nnn \exp_after:wN \@@_fixed_dtf_zeros:NN \exp_after:wN #1 \int_use:N \__int_eval:w 10 0000 + #2 \__int_eval_end: #3#4#5#6#7 ; 1 ; } \cs_new:Npn \@@_fixed_dtf_zeros:NN #1#2 { \if_meaning:w 0 #2 #1 \c_one \else: \@@_fixed_dtf_zeros_end:wNww #2 \fi: \@@_fixed_dtf_zeros:NN #1 } \cs_new:Npn \@@_fixed_dtf_zeros_end:wNww #1 \fi: \@@_fixed_dtf_zeros:NN #2 #3; #4 \s_@@ { \fi: \if_meaning:w ; #1 #2 \c_two * \c_@@_max_exponent_int \use_i_ii:nnn \fi: \@@_fixed_dtf_zeros_auxi:ww #1#3 0000 0000 0000 0000 0000 0000 ; } \cs_new:Npn \@@_fixed_dtf_zeros_auxi:ww { \@@_pack_twice_four:wNNNNNNNN \@@_pack_twice_four:wNNNNNNNN \@@_pack_twice_four:wNNNNNNNN \@@_fixed_dtf_zeros_auxii:ww ; } \cs_new:Npn \@@_fixed_dtf_zeros_auxii:ww #1; #2; #3 { #3 #1; } % \end{macrocode} % \newcommand{\eTeXfrac}[2]{\left[\frac{#1}{#2}\right]} % We get % \begin{quote} % \cs{@@_fixed_dtf_approx:n} \meta{B'} |;| \meta{A'} |;| % \end{quote} % where \meta{B'} and \meta{A'} are each $6$ brace groups, % representing fixed point numbers in the range $[0.1,1)$. Denote by % $x\in[1000,9999]$ and $y\in[0,9999]$ the first two groups of % \meta{B'}. We first find an estimate $a$ for the inverse of $B'$ by % computing % \begin{align*} % \alpha &= \eTeXfrac{10^{9}}{x+1} \\ % \beta &= \eTeXfrac{10^{9}}{x} \\ % a &= 10^{3} \alpha + (\beta-\alpha) \cdot % \left(10^{3}-\eTeXfrac{y}{10}\right) - 1750, % \end{align*} % where $\eTeXfrac{\bullet}{\bullet}$ denotes \eTeX{}'s rounding % division. The idea is to interpolate between $\alpha$ and $\beta$ % with a parameter $y/10^{4}$. The shift by $1750$ helps to ensure % that $a$ is an underestimate of the correct value. We will prove % that % \[ % 1 - 2.255\cdot 10^{-5} < \frac{B'a}{10^{8}} < 1 . % \] % We can then compute the inverse $B'a/10^{8}$ using $1/(1-\epsilon) % \simeq (1+\epsilon)(1+\epsilon^{2})$, which is correct up to a % relative error of $\epsilon^4 < 2.6\cdot 10^{-19}$. Since we target % a $16$-digit value, this is small enough. % % Let us prove the upper bound first. % \begin{align}\label{l3fp-fixed-eTeXfrac} % 10^{7} B'a % & < \left(10^{3} x + \eTeXfrac{y}{10} + \frac{3}{2}\right) % \left(\left(10^{3}-\eTeXfrac{y}{10}\right) \beta % + \eTeXfrac{y}{10} \alpha - 1750\right) % \\& < \left(10^{3} x + \eTeXfrac{y}{10} + \frac{3}{2}\right) % \left(\left(10^{3}-\eTeXfrac{y}{10}\right) % \left(\frac{10^{9}}{x} + \frac{1}{2} \right) % + \eTeXfrac{y}{10} \left(\frac{10^{9}}{x+1} + \frac{1}{2} \right) % - 1750\right) % \\& < \left(10^{3} x + \eTeXfrac{y}{10} + \frac{3}{2}\right) % \left(\frac{10^{12}}{x} % - \eTeXfrac{y}{10} \frac{10^{9}}{x(x+1)} % - 1250\right) % \end{align} % We recognize a quadratic polynomial in $[y/10]$ with a negative % leading coefficient, $([y/10]+a)(b-c[y/10]) \leq (b+ca)^2/(4c)$. % Hence, % \[ % 10^{7} B'a % < \frac{10^{15}}{x(x+1)} \left( % x + \frac{1}{2} + \frac{3}{4} 10^{-3} % - 6.25\cdot 10^{-10} x(x+1) \right)^2 % \] % We want to prove that the squared expression is less than $x(x+1)$, % which we do by simplifying the difference, and checking its sign, % \[ % x(x+1) - \left(x + \frac{1}{2} + \frac{3}{4} 10^{-3} % - 6.25\cdot 10^{-10} x(x+1) \right)^2 % > - \frac{1}{4} (1+1.5\cdot 10^{-3})^2 - 10^{-3} x % + 1.25\cdot 10^{-9} x(x+1)(x+0.5) % > 0. % \] % % Now, the lower bound. The same computation as % \eqref{l3fp-fixed-eTeXfrac} imply % \[ % 10^{7} B'a % > \left(10^{3} x + \eTeXfrac{y}{10} - \frac{1}{2}\right) % \left(\frac{10^{12}}{x} - \eTeXfrac{y}{10} \frac{10^{9}}{x(x+1)} % - 2250\right) % \] % This time, we want to find the minimum of this quadratic polynomial. % Since the leading coefficient is still negative, the minimum is % reached for one of the extreme values $y=0$ or $y=9999$, and we % easily check the bound for those values. % % We have proven that the algorithm will give us a precise enough % answer. Incidentally, the upper bound that we derived tells us that % $a < 10^{8}/B \leq 10^{9}$, hence we can compute $a$ safely as a % \TeX{} integer, and even add $10^{9}$ to it to ease grabbing of all % the digits. % \begin{macrocode} \cs_new:Npn \@@_fixed_dtf_approx:n #1 { \exp_after:wN \@@_fixed_dtf_approx:wnn \int_use:N \__int_eval:w 10 0000 0000 / ( #1 + \c_one ) ; {#1} } \cs_new:Npn \@@_fixed_dtf_approx:wnn #1; #2#3 { % \assert:n { \tl_count:n {#1} = 6 } \exp_after:wN \@@_fixed_dtf_approx:NNNNNw \int_use:N \__int_eval:w 10 0000 0000 - 1750 + #1000 + (10 0000 0000/#2-#1) * (1000-#3/10) ; {#2}{#3} } \cs_new:Npn \@@_fixed_dtf_approx:NNNNNw 1#1#2#3#4#5#6; #7; #8; { + \c_four % because of the line below "dtf_epsilon" here. \@@_fixed_mul:wwn {000#1}{#2#3#4#5}{#6}{0000}{0000}{0000} ; #7; \@@_fixed_dtf_epsilon:wN \@@_fixed_mul:wwn {000#1}{#2#3#4#5}{#6}{0000}{0000}{0000} ; \@@_fixed_mul:wwn #8; \@@_fixed_to_float:wN ? } \cs_new:Npn \@@_fixed_dtf_epsilon:wN #1#2#3#4#5#6; { % \assert:n { #1 = 0000 } % \assert:n { #2 = 9999 } \exp_after:wN \@@_fixed_dtf_epsilon:NNNNNww \int_use:N \__int_eval:w 1 9999 9998 - #3#4 + \exp_after:wN \@@_fixed_dtf_epsilon_pack:NNNNNw \int_use:N \__int_eval:w 2 0000 0000 - #5#6 ; {0000} ; } \cs_new:Npn \@@_fixed_dtf_epsilon_pack:NNNNNw #1#2#3#4#5#6; { #1 ; {#2#3#4#5} {#6} } \cs_new:Npn \@@_fixed_dtf_epsilon:NNNNNww #1#2#3#4#5#6; #7; { \@@_fixed_mul:wwn %^^A todo: optimize to use \@@_mul_significand. {0000} {#2#3#4#5} {#6} #7 ; {0000} {#2#3#4#5} {#6} #7 ; \@@_fixed_add_one:wN \@@_fixed_mul:wwn {10000} {#2#3#4#5} {#6} #7 ; } % \end{macrocode} % \end{macro} % % \begin{macrocode} % % \end{macrocode} % % \end{implementation} % % \PrintChanges % % \PrintIndex