diff options
author | Karl Berry <karl@freefriends.org> | 2012-07-23 17:15:00 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2012-07-23 17:15:00 +0000 |
commit | 134349701bddf7cbbacf6030c6b9f9838aff96fa (patch) | |
tree | ec1140c46e1c0347a671a6fa3cf8af5a79e95f93 /Master/texmf-dist/source/latex/l3kernel/l3fp-extended.dtx | |
parent | f7855c12c18bb97b7b9e49ab685ee558d8c0b47b (diff) |
l3kernel 3990 (17jul12)
git-svn-id: svn://tug.org/texlive/trunk@27108 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/source/latex/l3kernel/l3fp-extended.dtx')
-rw-r--r-- | Master/texmf-dist/source/latex/l3kernel/l3fp-extended.dtx | 674 |
1 files changed, 674 insertions, 0 deletions
diff --git a/Master/texmf-dist/source/latex/l3kernel/l3fp-extended.dtx b/Master/texmf-dist/source/latex/l3kernel/l3fp-extended.dtx new file mode 100644 index 00000000000..d38c6de8e17 --- /dev/null +++ b/Master/texmf-dist/source/latex/l3kernel/l3fp-extended.dtx @@ -0,0 +1,674 @@ +% \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{l3names} +\GetIdInfo$Id: l3fp-extended.dtx 2474 2011-06-17 12:54:02Z bruno $ + {L3 Floating-point extended precision fixed-points} +\documentclass[full]{l3doc} +\begin{document} + \DocInput{\jobname.dtx} +\end{document} +%</driver> +% \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} +% +% 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, except +% \meta{a_1}, which may be any positive \TeX{} integer. 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} |;| \meta{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. +% 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. +% +% \begin{variable}{\c_@@_one_fixed_tl} +% \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. +% \begin{macrocode} +\cs_new:Npn \@@_fixed_continue:wn #1; #2 { #2 #1; } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[int, EXP]{\@@_fixed_div_int:wwN} +% \begin{macro}[aux, EXP] +% { +% \@@_fixed_div_int_i:wnN, \@@_fixed_div_int_ii:wnn, +% \@@_fixed_div_int_end:wnn, \@@_fixed_div_int_pack:Nw, +% \@@_fixed_div_int_after:Nw +% } +% \begin{syntax} +% \cs{@@_fixed_div_int:wwN} \meta{a} |;| \meta{n} |;| \meta{function} +% \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{function}. The +% \texttt{wnN} auxiliary receives $a_{i}$, $n$, and a continuation +% function as arguments, and computes a (rather tight) lower bound +% $Q_{i}$ for the quotient. The \texttt{wnn} auxiliary receives +% $Q_{i}$, $n$, and $a_{i}$. It adds $Q_{i}$ to a surrounding integer +% expression, and starts a new one. It also computes $a_{i}-n\cdot +% Q_{i}$, putting the result in front of $a_{i+1}$ to serve as the +% first argument for a new call to the \texttt{wnN} auxiliary. At the +% end, the path we took to the lowest levels rewinds: the +% \texttt{pack} auxiliary receives $5$ digits, braces the last $4$, +% and carries the leading digit to the level above. The offsets used +% to ensure a given number of digits are as follows: we first subtract +% $1$ from the top-level, then add $9999$ at every subsequent level, +% and add $2$ to the last level. This last number is not $1$, because +% it compensates for the |- \c_one| in the \texttt{wnN} auxiliary. +% \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_i:wnN + #1; {#7} \@@_fixed_div_int_ii:wnn + #2; {#7} \@@_fixed_div_int_ii:wnn + #3; {#7} \@@_fixed_div_int_ii:wnn + #4; {#7} \@@_fixed_div_int_ii:wnn + #5; {#7} \@@_fixed_div_int_ii:wnn + #6; {#7} \@@_fixed_div_int_end:wnn ; + } +\cs_new:Npn \@@_fixed_div_int_i: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_ii: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_i:wnN + \int_use:N \__int_eval:w #3 - #1*#2 \__int_eval_end: + } +\cs_new:Npn \@@_fixed_div_int_end: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} +% +% \begin{macro}[int, EXP]{\@@_fixed_add_one:wN} +% \begin{syntax} +% \cs{@@_fixed_add_one:wN} \Arg{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4} \Arg{X_5} \Arg{X_6} |;| \meta{function} +% \end{syntax} +% \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 10000 + #1 } #2 ; + } +% \end{macrocode} +% \end{macro} +% +% \begin{macro}[int, EXP] +% {\@@_fixed_add:wwN, \@@_fixed_sub:wwN, \@@_fixed_sub_back:wwN} +%^^A todo: remove sub_back. +% \begin{macro}[aux, EXP] +% { +% \@@_fixed_add_i:NNnnnnwnn, +% \@@_fixed_add_ii:NnnNnnnnw, +% \@@_fixed_add_pack:NNNNNwN, +% \@@_fixed_add_after:NNNNNwN +% } +% \begin{syntax} +% \cs{@@_fixed_add:wwN} \Arg{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4} \Arg{X_5} \Arg{X_6} |;| \Arg{Y_1} \Arg{Y_2} \Arg{Y_3} \Arg{Y_4} \Arg{Y_5} \Arg{Y_6} |;| \meta{function} +% \end{syntax} +% Computes $X+Y$ (resp.\ $X-Y$ and $Y-X$) and feeds the result to +% \meta{function}. The three functions only differ by some signs and +% use a common auxiliary. It would be nice to grab the $12$ brace +% groups in one go, only $9$ arguments are allowed. Start by grabbing +% the two signs, $X_{1}, \ldots, X_{4}$, the rest of $X$, and $Y_{1}$ +% and $Y_{2}$. The second auxiliary receives the sign of $X$, the +% rest of $X$, the sign of $Y$, the rest of $Y$, and the +% \meta{function}. After going down through the various level, we go +% back up, packing digits and bringing the \meta{function} (|#9|, then +% |#7|) from the end of the argument list to its start. +% \begin{macrocode} +\cs_new_nopar:Npn \@@_fixed_add:wwN { \@@_fixed_add_i:NNnnnnwnn + + } +\cs_new_nopar:Npn \@@_fixed_sub:wwN { \@@_fixed_add_i:NNnnnnwnn + - } +\cs_new_nopar:Npn \@@_fixed_sub_back:wwN { \@@_fixed_add_i:NNnnnnwnn - + } +\cs_new:Npn \@@_fixed_add_i:NNnnnnwnn #1#2 #3#4#5#6 #7; #8#9 + { + \exp_after:wN \@@_fixed_add_after:NNNNNwN + \int_use:N \__int_eval:w 1 9999 9998 #1 #3#4 #2 #8#9 + \exp_after:wN \@@_fixed_add_pack:NNNNNwN + \int_use:N \__int_eval:w 1 9999 9998 #1 #5#6 + \@@_fixed_add_ii:NnnNnnnnw #1 #7 #2 + } +\cs_new:Npn \@@_fixed_add_ii:NnnNnnnnw #1 #2#3 #4 #5#6 #7#8 ; #9 + { + #4 #5#6 + \exp_after:wN \@@_fixed_add_pack:NNNNNwN + \int_use:N \__int_eval:w 2 0000 0000 #4 #7#8 #1 #2#3 ; #9 ; + } +\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 #2#3#4#5 #6; #7 + { + \exp_after:wN #7 + \exp_after:wN { \int_use:N \__int_eval:w - 2 0000 + #1#2#3#4#5 } + {#6} + } +% \end{macrocode} +% \end{macro} +% \end{macro} +% +% \begin{macro}[int, EXP]{\@@_fixed_mul:wwn} +% \begin{macro}[aux, EXP] +% { +% \@@_fixed_mul_i:nnnnnnnn , +% \@@_fixed_mul_pack:NNNNNw , +% \@@_fixed_mul_after:wwn +% } +% \begin{syntax} +% \cs{@@_fixed_mul:wwn} \Arg{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4} \Arg{X_5} \Arg{X_6} |;| \Arg{Y_1} \Arg{Y_2} \Arg{Y_3} \Arg{Y_4} \Arg{Y_5} \Arg{Y_6} |;| \Arg{tokens} +% \end{syntax} +% Computes $X\times Y$ and feeds the result to \meta{function}. It +% would be nice to grab the $12$ brace groups in one go, but that's +% not possible. On the other hand, we don't need to obtain an exact +% rounding, contrarily to the case in \cs{@@_*_o:ww}, so things are +% not quite as bad as they may seem. The parenthesis computing the +% seventh group of digits (computed because we need to know its +% potentially large carry) is closed by +% \cs{@@_fixed_mul_i:nnnnnnnn}, once we access the last two brace +% groups, which were not read before. Also, in +% \cs{@@_fixed_mul_after:wwn}, |#3| is the continuation +% tokens.\footnote{Bruno: insist on the difference compared to +% \cs{@@_fixed_add:wwN}.} +% \begin{macrocode} +\cs_new:Npn \@@_fixed_mul:wwn #1#2#3#4 #5; #6#7#8#9 + { + \exp_after:wN \@@_fixed_mul_after:wwn + \int_use:N \__int_eval:w \c_@@_leading_shift_int + \exp_after:wN \@@_pack:NNNNNw + \int_use:N \__int_eval:w \c_@@_middle_shift_int + + #1*#6 + \exp_after:wN \@@_pack:NNNNNw + \int_use:N \__int_eval:w \c_@@_middle_shift_int + + #1*#7 + #2*#6 + \exp_after:wN \@@_pack:NNNNNw + \int_use:N \__int_eval:w \c_@@_middle_shift_int + + #1*#8 + #2*#7 + #3*#6 + \exp_after:wN \@@_pack:NNNNNw + \int_use:N \__int_eval:w \c_@@_middle_shift_int + + #1*#9 + #2*#8 + #3*#7 + #4*#6 + \exp_after:wN \@@_pack:NNNNNw + \int_use:N \__int_eval:w \c_@@_trailing_shift_int + + #2*#9 + #3*#8 + #4*#7 + + ( #3*#9 + #4*#8 + + \@@_fixed_mul_i:nnnnnnnn #5 {#6}{#7} {#1}{#2} + } +\cs_new:Npn \@@_fixed_mul_i:nnnnnnnn #1#2 #3#4 #5#6 #7#8 + { #1*#4 + #2*#3 + #5*#8 + #6*#7 )/10000 + #1*#3 + #5*#7 ; } +\cs_new:Npn \@@_fixed_mul_pack:NNNNNw + #1 #2#3#4#5 #6; { + #1#2#3#4#5 ; {#6} } +\cs_new:Npn \@@_fixed_mul_after:wwn #1; #2; #3 { #3 {#1} #2 ; } +% \end{macrocode} +% \end{macro} +% \end{macro} +% +% \begin{macro}[int, EXP]{\@@_fixed_mul_add:wwwn, \@@_fixed_mul_sub_back:wwwn} +% \begin{syntax} +% \cs{@@_fixed_mul_add:wwn} \Arg{X_1} \Arg{X_2} \Arg{X_3} \Arg{X_4} \Arg{X_5} \Arg{X_6} |;| \Arg{Y_1} \Arg{Y_2} \Arg{Y_3} \Arg{Y_4} \Arg{Y_5} \Arg{Y_6} |;| \Arg{Z_1} \Arg{Z_2} \Arg{Z_3} \Arg{Z_4} \Arg{Z_5} \Arg{Z_6} |;| \Arg{tokens} +% \end{syntax} +% These functions compute $X\times Y + Z$ or $Z-X\times Y$ and feed +% the result to the \meta{tokens}. This is tough because we have $18$ +% brace groups in front of us. +% \begin{macrocode} +\cs_new:Npn \@@_fixed_one_minus_mul:wwn #1; #2#3#4#5; + { + \exp_after:wN \@@_fixed_mul_after:wwn + \int_use:N \__int_eval:w \c_@@_big_leading_shift_int + \c_ten_thousand + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + \@@_fixed_mul_add_i:Nnwnnwnnn + - 00; {#2}{#3}{#4}; #1; {#2}{#3}{#4}#5; - 00 ; + } +\cs_new:Npn \@@_fixed_mul_add:wwwn #1; #2#3#4#5; #6#7#8#9 + { + \exp_after:wN \@@_fixed_mul_after:wwn + \int_use:N \__int_eval:w \c_@@_big_leading_shift_int + #6 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #7 + \@@_fixed_mul_add_i:Nnwnnwnnn + + {#8}{#9}; {#2}{#3}{#4}; #1; {#2}{#3}{#4}#5; + + } +\cs_new:Npn \@@_fixed_mul_sub_back:wwwn #1; #2#3#4#5; #6#7#8#9 + { + \exp_after:wN \@@_fixed_mul_after:wwn + \int_use:N \__int_eval:w \c_@@_big_leading_shift_int + #6 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #7 + \@@_fixed_mul_add_i:Nnwnnwnnn + - {#8}{#9}; {#2}{#3}{#4}; #1; {#2}{#3}{#4}#5; - + } +\cs_new:Npn \@@_fixed_mul_add_i:Nnwnnwnnn #1 #2#3; #4#5#6; #7#8#9 + { % sg z3z4; y1y2y3; x1x2x3 x4x5x6; y1y2y3y4y5y6; sg z5z6; + #1 #7*#4 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #2 + #1 #7*#5 #1 #8*#4 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #3 + #1 #7*#6 #1 #8*#5 #1 #9*#4 + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_middle_shift_int + #1 \@@_fixed_mul_add_ii:nnnnwnnnn {#7}{#8}{#9} + } +\cs_new:Npn \@@_fixed_mul_add_ii:nnnnwnnnn #1#2#3#4#5; #6#7#8#9 + { % x1x2x3x4 x5x6; y1y2y3y4 y5y6; sg z5z6; + ( #1*#9 + #2*#8 + #3*#7 + #4*#6 ) + \exp_after:wN \@@_pack_big:NNNNNNw + \int_use:N \__int_eval:w \c_@@_big_trailing_shift_int + \@@_fixed_mul_add_iii:nnnnwnnwN + { #6 + #4*#7 + #3*#8 + #2*#9 + #1 } + { #7 + #4*#8 + #3*#9 + #2 } + {#1} #5; + {#6} + } +\cs_new:Npn \@@_fixed_mul_add_iii:nnnnwnnwN #1#2 #3#4#5; #6#7#8; #9 + { % {y1+x4*y2+x3*y3+x2*y4+x1} {y2+x4*y3+x3*y4+x2} + % x1x5x6; y1y5y6; sg z5z6; + % => + % sg (x5*y1+x4*y2+x3*y3+x2*y4+x1*y5) + % sg (x6*y1+x5*y2+x4*y3+x3*y4+x2*y5+x1*y6)/10000 + % + z5z6; + #9 (#4* #1 *#7) + #9 (#5*#6+#4* #2 *#7+#3*#8) / \c_ten_thousand + + \@@_use_s:nn + } +% \end{macrocode} +% \end{macro} +% +% \begin{macrocode} +\cs_new:Npn \@@_fixed_to_float:Nw #1#2; { \@@_fixed_to_float:wN #2; #1 } +% \end{macrocode} +% +% \begin{macro}[int, rEXP]{\@@_fixed_to_float:wN} +% \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: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 + \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_ii:ww + #1#3 0000 0000 0000 0000 0000 0000 ; + } +\cs_new:Npn \@@_fixed_dtf_zeros_ii:ww + { + \@@_pack_twice_four:wNNNNNNNN + \@@_pack_twice_four:wNNNNNNNN + \@@_pack_twice_four:wNNNNNNNN + \@@_fixed_dtf_zeros_iii:ww + ; + } +\cs_new:Npn \@@_fixed_dtf_zeros_iii: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_ii:wnn + \int_use:N \__int_eval:w 10 0000 0000 / ( #1 + \c_one ) ; + {#1} + } +\cs_new:Npn \@@_fixed_dtf_approx_ii:wnn #1; #2#3 + { +%<assert> \assert:n { \tl_count:n {#1} = 6 } + \exp_after:wN \@@_fixed_dtf_approx_iii: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_iii: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> \assert:n { #1 = 0000 } +%<assert> \assert:n { #2 = 9999 } + \exp_after:wN \@@_fixed_dtf_epsilon_ii: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_ii:NNNNNww #1#2#3#4#5#6; #7; + { + \@@_fixed_mul:wwn %^^A todo: optimize to use \@@_mul_mantissa. + {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} +%</initex|package> +% \end{macrocode} +% +% \end{implementation} +% +% \PrintChanges +% +% \PrintIndex
\ No newline at end of file |