EDIT (three years later) Answering this question has led to a complete reimplementation of floating point calculations in LaTeX3.
There are two types of (complex) macros in TeX: those ("expandable" macros) which work purely by expansion, turning their argument into something else without ever needing to store an intermediate result or otherwise affect the state of the TeX machine, and those ("non-expandable" macros) which alter the state.
TeX proper only computes with integers non-expandably, through \advance, \multiply and \divide, essentially +=, *=, and //=. These assignments are non-expandable. Instead of simply doing 2+3*4+5 or \add{\add{2}{\mul{3}{4}}}{5} one must do the equivalent of
x = 2
y = 3
y *= 4
x += y
x += 5
using intermediate variables, which is unpractical. The eTeX primitive \numexpr lets us evaluate expandably expressions which involve integers, (, ), +, -, *, /. As Joseph describes, it lets us define function-like macros
\def\add#1#2{\the\numexpr(#1)*(#2)\relax}
\def\abs#1{\the\numexpr\ifnum\numexpr#1<0-\fi(#1)\relax}
\add{\abs{-3}}{4} % = 7
Many packages (pgfmath, fp, ...) provide non-expandable macros to compute with floating point numbers, with varying levels of practicality. Since last year, l3fp evaluates floating point expressions expandably. A typical example:
\fp_eval:n { abs(-3)+4*ln(5.6e7) }
For historical purposes, I've left the code for a prototype of l3fp below.
EDIT (2011) I had written a parsing step for more user-friendliness, but it was buggy, so I removed it. I also rewrote the code completely a few times and implemented subtraction and division. The code:
\documentclass{article}
\usepackage{expl3}
\ExplSyntaxOn
% ====
% Floating point numbers which we'll call ``<ZBfp>'', of the form
% <sign> {<int>}{<4d>}{<4d>}{<4d>}T{<int expr>}
% where
% <sign> is "+" or "-" (mandatory)
% <int> is at most 4 digits, and never equal to zero (except cases below)
% <4d> represents 4 digits exactly
% <int expr> is something that "\number\numexpr...\relax" is happy with.
%
% Exceptions:
% - The <sign> can also be "X" for undefined.
% X{0000}{0000}{0000}{0000}T{0}
% This is converted to a sign of "+" if we use it further in the calculation.
% +{0000}{0000}{0000}{0000}T{0}
%
% Example:
% + {0234}{2345}{3456}{4567}T{0} means 234.234534564567
% - {1111}{2222}{3333}{4444}T{2} means -111122223333.4444
%
% ==== Helpers
\cs_new:Npn \ZBfp_brace_last:nw #1 #2 #{#1{#2}}
\cs_new:Npn \ZBfp_brace_last:nnnnnw #1 #2 #3 #4 #5 #6 # { #1 #2 #3 #4 #5 {#6} }
\cs_new:Npn \use_Bi_iiE:nn #1 #2 { {#1 #2} }
\cs_new:Npn \use_iv_delimit_by_q_stop:nnnnw #1#2#3#4#5\q_stop {#4}
\tl_const:Nn \c_ZBfp_braced_plus_tl { {+} }
\tl_const:Nn \c_ZBfp_braced_minus_tl { {-} }
% Constants
\tl_const:Nn \c_zero_ZBfp {+{0000}{0000}{0000}{0000}T{0}}
\tl_const:Nn \c_undefined_ZBfp {X{0000}{0000}{0000}{0000}T{0}}
% ==== Adding mantissas.
%
% Used as
% \ZBfp_add_mantissa:nnnnnnnn {<4d>}{<4d>}{<4d>}{<4d>}
% {<4d>}{<4d>}{<4d>}{<4d>} {<sign>} T {<exp>}
% to get
% <sign> {<4d>}{<4d>}{<4d>}{<4d>} T{<exp>}
%
% We split numbers efficiently by pre-carrying.
% Note that all the extra "-1", "+9999", "10000" sum up to zero.
% Also, it would be more efficient to replace those constants by registers?
\cs_new:Npn \ZBfp_add_mantissa:nnnnnnnn #1#2#3#4 #5#6#7#8{
\exp_after:wN \ZBfp_add_mantissa_aux:N
\tex_the:D \etex_numexpr:D -1 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D #1 + #5 + 9999 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D #2 + #6 + 9999 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D #3 + #7 + 9999 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D #4 + #8 + 10000
}
% Important: that "#8" above is followed by a brace group, \marg{sign}.
% We are left with
% "\ZBfp_add_mantissa_aux:N" <1d>{4d}{4d}{4d}{4d}{<sign>}T{<exp>}
% where <1d> is either "0" or "1", the overall carry.
\cs_new:Npn \ZBfp_add_mantissa_aux:N #1 {
\if_meaning:w 0 #1
\exp_after:wN \ZBfp_add_mantissa_aux_z:nnnnnnw
\else:
\exp_after:wN \ZBfp_add_mantissa_aux_i:nnnnnnw
\fi:
#1
}
% If the first digit is "0", we discard it and keep four brace groups.
% Otherwise, we keep it in a brace group, together with three more.
% In both cases, we place the sign and the exponent at the right place.
% Note that the exponent is not evaluated.
\cs_new:Npn \ZBfp_add_mantissa_aux_z:nnnnnnw #1 #2 #3 #4 #5 #6 T #7 {
#6 {#2}{#3}{#4}{#5} T{#7} }
\cs_new:Npn \ZBfp_add_mantissa_aux_i:nnnnnnw #1 #2 #3 #4 #5 #6 T #7 {
#6 {000 #1}{#2}{#3}{#4} T{#7+1} }
% ===== Subtracting mantissa
% Used as
% \ZBfp_sub_back_mantissa:nnnnnnnn <A> <B> {<sign>} T {<exp>}
% where <A> and <B> have the form {<4d>}{<4d>}{<4d>}{<4d>}. This
% calculates $<sign>*( - <A> + <B> )$ (hence the name "back") and
% f-expands to the result in the form
% <sign> {<4d>}{<4d>}{<4d>}{<4d>} T{<exp>}
%
% Again, pre-carrying. We check after the calculation if the
% subtraction should have been done in the other direction.
% Thus, doing subtraction as - <big mantissa> + <small mantissa>
% should be twice slower than the other way around, but that's
% the price to pay to get a more optimized subtraction
% - <small> + <big>, used more frequently in the ``full''
% subtraction and in sines.
%
% We use bigger shifts that we had for addition, because we want
% to ensure that every intermediate numexpr we use has five digits
% for "\ZBfp_brace_last:nw" to work.
\cs_new:Npn \ZBfp_sub_back_mantissa:nnnnnnnn #1#2#3#4 #5#6#7#8 {
\exp_after:wN \ZBfp_sub_back_mantissa_aux:w
\tex_the:D \etex_numexpr:D #5 - #1 - 2 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D #6 - #2 + 19998 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D #7 - #3 + 19998 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D #8 - #4 + 20000
}
% Again, following "#8" is a brace group, \marg{sign}.
% We are left with
% "\ZBfp_sub_back_mantissa_aux:w" <w>{4d}{4d}{4d}{<sign>}T{<exp>},
% where <w> is any integer between -9999 and +9999. Here we need to take
% care of the fact that e.g. -11{2222}{2222}{2222} means
% "-11 + .222222222222", not "-11.222222222222"! Convince yourself
% that if <w> is zero, then we are in the ``positive'' case.
\cs_new:Npn \ZBfp_sub_back_mantissa_aux:w #1 #{
\if_num:w #1 < 0 \exp_stop_f:
\exp_after:wN \ZBfp_sub_back_mantissa_neg:nnnnn
\else:
\exp_after:wN \ZBfp_sub_back_mantissa_pos:nnnn
\fi:
{#1}
}
% \ZBfp_sub_back_mantissa_neg:nnnnn {-<4d>}{<4d>}{<4d>}{<4d>} {<sign>} T {<exp>}
% If the result of the subtraction was negative, then we need
% to correct it and correct the sign.
%
% > "#1", "#2", "#3", "#4" form the body of the number
% > "#5" is the <sign>
% > "#6" would be the exponent expression
%
% Here, we know that "#1" is negative. In fact, our number is
% equal to <sign>"(#1 + .#2#3#4)". We calculate "|#1|-.#2#3#4",
% which is positive, and place the correct sign (namely, $-$<sign>)
% in a brace group after it, then pass the whole thing to
% "\ZBfp_add_mantissa_pos:nnnnnw".
%
% Once more, pre-carry, to make sure everything has 5 digits when passed
% to "\ZBfp_brace_last:nw".
%
% On the first piece below, "-2 - #1 + ...". From the level below
% (i.e. the "..."), we get "1", or sometimes "2", and we know that
% "10000<#1<0", so the whole thing is "0<=...<10000", at most 4 digits.
%
\cs_new:Npn \ZBfp_sub_back_mantissa_neg:nnnnn #1 #2 #3 #4 #5 {
\exp_after:wN \ZBfp_brace_last:nw
\exp_after:wN \ZBfp_sub_back_mantissa_pos:nnnn
\tex_the:D \etex_numexpr:D -2 - #1 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D 19998 - #2 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D 19998 - #3 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D 20000 - #4
\if:w -#5
\exp_after:wN \c_ZBfp_braced_plus_tl
\else:
\exp_after:wN \c_ZBfp_braced_minus_tl
\fi:
}
% Note that the final line gets fully expanded to "{+}" or "{-}"
% when \TeX\ is trying to finish the number "#4". Specifically, the sign
% we get is the opposite of "#5". Any faster test (braces are important)?
% "\ZBfp_sub_back_mantissa_pos:nnnn" {<4d>}{<4d>}{<4d>}{<4d>} {<sign>} T {<exp>}
% Once we have taken care of the signs, we need to check whether
% the first, second, etc. pieces is/are zero, and shift exponents
% accordingly.
% > "#1", "#2", "#3", "#4" form the body of the number,
% > "#5" is the sign
% > "#7" is the exponent
\cs_new:Npn \ZBfp_sub_back_mantissa_pos:nnnn #1 #2 #3 #4 {
\cs:w ZBfp_sub_maux_
\if_num:w #1 = 0 \exp_stop_f: i
\if_num:w #2 = 0 \exp_stop_f: i
\if_num:w #3 = 0 \exp_stop_f: i
\if_num:w #4 = 0 \exp_stop_f: i
\fi:
\fi:
\fi:
\fi:
:w
\cs_end:
{#1}
{#2}
{#3}
{#4}
}
% Just reformat digits in the right way.
\cs_new:Npx\ZBfp_sub_maux_iiii:w #1#2#3#4 #5T#6 {\c_zero_ZBfp}
\cs_new:Npn\ZBfp_sub_maux_iii:w#1#2#3#4 #5T#6 {#5 {#4}{0000}{0000}{0000}T{#6-3}}
\cs_new:Npn\ZBfp_sub_maux_ii:w #1#2#3#4 #5T#6 {#5 {#3} {#4} {0000}{0000}T{#6-2}}
\cs_new:Npn\ZBfp_sub_maux_i:w #1#2#3#4 #5T#6 {#5 {#2} {#3} {#4} {0000}T{#6-1}}
\cs_new:Npn\ZBfp_sub_maux_:w #1#2#3#4 #5T#6 {#5 {#1} {#2} {#3} {#4} T {#6} }
% ======= Full addition and subtraction.
% Syntax:
% "\ZBfp_add:nn" \marg{ZBfp expr.} \marg{ZBfp expr.}
% "\ZBfp_sub:nn" \marg{ZBfp expr.} \marg{ZBfp expr.}
% where <ZBfp> are floating point expressions, which will be
% fully f-expanded before the calculation.
% Plan:
% > expand the left <ZBfp expr.>
% > swap, and expand the right <ZBfp expr.>
% > check signs
% > shift digits according to exponent
\cs_new:Npn \ZBfp_add:nn #1 {
\exp_after:wN \ZBfp_add_aux:NwNn
\tex_romannumeral:D -`\0 #1 \q_stop +
}
\cs_new:Npn \ZBfp_sub:nn #1 {
\exp_after:wN \ZBfp_add_aux:NwNn
\tex_romannumeral:D -`\0 #1 \q_stop -
}
% We now have
% "\ZBfp_add_aux:NwNn" <sign A> <A> T {<eA>} \q_stop <operation> {<expr B>}
% Here <sign A> <A> T {<eA>} is our first floating point number:
% <A> is the body, {<4d>}{<4d>}{<4d>}{<4d>}, and <eA> is the exponent.
% Also, <operation> is "+" or "-".
% We then move <expr B> to the front and f-expand it.
\cs_new:Npn \ZBfp_add_aux:NwNn #1 #2 \q_stop #3 #4 {
\exp_after:wN \ZBfp_add_signs:NNN \exp_after:wN #1 \exp_after:wN #3
\tex_romannumeral:D -`\0 #4 #2
}
% We now have
% "\ZBfp_add_aux:NwNn" <sign A> <operation> <sign B> <B>T{<eB>} <A>T{<eA>}
% If the product <sign A><operation><sign B> is +,
% use "\ZBfp_add_exponent:w", else "\ZBfp_sub_exponent:w".
\cs_new:Npn \ZBfp_add_signs:NNN #1 #2 #3 {
\cs:w ZBfp_
\if:w #1 \if:w #2 #3 + \else: - \fi: add \else: sub \fi:
_exponent:w \cs_end:
#1
}
% We now have "\ZBfp_..._exponent:w" <sign A> <B>T{<eB>} <A>T{<eA>}
% Take care of exponents. If they are the same, do nothing.
% If they differ, shift in the relevant direction using
% "\ZBfp_decimate_do". We first check the equality case
% because we will arrange to be in that case for our calculations
% of sine and cosine (which we want to optimize).
\cs_new:Npn \ZBfp_add_exponent:w #1 #2 T #3 #4 T #5 {
\if_num:w #3 > #5 \exp_stop_f:
\exp_after:wN \use_i:nn
\else:
\exp_after:wN \use_ii:nn
\fi:
{
\exp_after:wN \ZBfp_decimate_do:nNnnnn \exp_after:wN {
\tex_the:D \etex_numexpr:D #3 - #5 }
\ZBfp_add_mantissa:nnnnnnnn #4 #2 {#1} T {#3}
}
{
\exp_after:wN \ZBfp_decimate_do:nNnnnn \exp_after:wN {
\tex_the:D \etex_numexpr:D #5 - #3 }
\ZBfp_add_mantissa:nnnnnnnn #2 #4 {#1} T {#5}
}
}
% "\ZBfp_decimate_do:nNnnnn {<shift>} <func> {<4d>}{<4d>}{<4d>}{<4d>}
% will shift the blocks of digits by <shift> blocks to the right,
% adding "{0000}" blocks at the start when shifting, then place
% <func> in front of 4 blocks of 4 digits. Requires $<shift> \geq 0$.
% When shifting results in a zero number, instead of performing the
% computation, we simply remove <func> and the 4 blocks of digits,
% presumably leaving only the other operand on the input stream.
% In other words, if the difference in exponent is ${}>3$, the small
% number is to small ($10000^{-4}$), and we throw it away.
\cs_new:Npn \ZBfp_decimate_do:nNnnnn #1 {
\if_case:w #1 \exp_stop_f:
% zero shift: do nothing.
\or: \exp_after:wN \ZBfp_decimate_i:Nnnnn
\or: \exp_after:wN \ZBfp_decimate_ii:Nnnnn
\or: \exp_after:wN \ZBfp_decimate_iii:Nnnnn
\else: \exp_after:wN \use_none:nnnnn
\fi:
}
% "\ZBfp_decimate_1:w" <func> {<4d>} {<4d>} {<4d>} {<4d>}
\cs_new:Npn \ZBfp_decimate_i:Nnnnn #1 #2#3#4#5 {#1 {0000} {#2} {#3} {#4}}
\cs_new:Npn \ZBfp_decimate_ii:Nnnnn #1 #2#3#4#5 {#1 {0000}{0000} {#2} {#3}}
\cs_new:Npn \ZBfp_decimate_iii:Nnnnn #1 #2#3#4#5 {#1 {0000}{0000}{0000}{#2}}
% "\ZBfp_sub_exponent:w" <sign A> <B>T{<eB>} <A>T{<eA>}
% For subtraction, signs are tricky. Given the above, we want to
% calculate $<sign A>*(<A>*10000^{<eA>} - <B>*10000^{<eB>})$.
% \begin{itemize}
% \item
% If $ <eB> \leq <eA> $, then we will ``decimate'' <B> (possibly
% with a zero <shift> when there is equality), and essentially
% calculate "sub_back" <B'> <A> <sign A> T {<eA>}, where <B'> is
% the result of decimating. That gives
% $ <sign A> * (- <B> + <A>) * 10000^{<eA>} $.
% \item
% Otherwise, we need to change the sign. Unfortunately, that sign
% has to be sent quite deep (just before the "T" of the exponent).
% Maybe that could be improved, but right now, we first calculate
% the correct sign and place it at the right spot using
% "\ZBfp_sub_exponent_aux:Nw". Some will complain that the
% conditional is not terminated when placing the sign, but no worry,
% it gets terminated since we just continue executing things after
% placing the sign.
%
% Then decimate <A>. And calculate "sub_back" <A'> <B> {- <sign A>} T {<eB>},
% which is equal to $ - <sign A> *(-<A>+<B>) * 10000^{<eB>}$.
% \end{itemize}
\cs_new:Npn \ZBfp_sub_exponent:w #1 #2 T #3 #4 T #5 {
\if_num:w #3 > #5 \exp_stop_f:
\exp_after:wN \use_i:nn
\else:
\exp_after:wN \use_ii:nn
\fi:
{
\exp_after:wN \ZBfp_sub_exponent_aux:Nw \if:w -#1 + \else: - \fi:
\exp_after:wN \ZBfp_decimate_do:nNnnnn \exp_after:wN {
\tex_the:D \etex_numexpr:D #3 - #5 }
\ZBfp_sub_back_mantissa:nnnnnnnn #4 #2 T {#3}
}
{
\exp_after:wN \ZBfp_decimate_do:nNnnnn \exp_after:wN {
\tex_the:D \etex_numexpr:D #5 - #3 }
\ZBfp_sub_back_mantissa:nnnnnnnn #2 #4 {#1} T {#5}
}
}
\cs_new:Npn \ZBfp_sub_exponent_aux:Nw #1 #2 T {#2 {#1} T}
% ===== Multiplication
% Syntax:
% "\ZBfp_mul:nn" \marg{ZBfp expr.} \marg{ZBfp expr.}
% First expand the first expression, then swap, and expand the
% second expression. Then treat exponents and place the sign at
% the right place.
\cs_new:Npn \ZBfp_mul:nn #1 {
\exp_after:wN \ZBfp_mul_aux:Nwn
\tex_romannumeral:D -`\0 #1 \q_stop
}
\cs_new:Npn \ZBfp_mul_aux:Nwn #1 \q_stop #2 {
\exp_after:wN \ZBfp_mul_signs_expo:NwnNwn
\tex_romannumeral:D -`\0 #2 #1
}
% The signs are "#1" and "#4", "#3" and "#6" are the exponents, and
% each "#2" and "#5" are 16 digits split in blocks of 4.
\cs_new:Npn \ZBfp_mul_signs_expo:NwnNwn #1 #2 T#3 #4 #5 T#6 {
\exp_after:wN \ZBfp_mul_mantissa:wnnnnnnnn
\exp_after:wN #1
\exp_after:wN #4
\exp_after:wN T
\exp_after:wN {\tex_the:D \etex_numexpr:D #3 + #6}
\q_stop
#2 #5
}
% Now we have
% "\ZBfp_mul_mantissa:wnnnnnnnn" <sign><sign> T {<exp>} \q_stop
% {<4d>}{<4d>}{<4d>}{<4d>} {<4d>}{<4d>}{<4d>}{<4d>}
% TODO: Round instead of truncating, by adapting one of the "99990000" shifts
% to something like "50000000" and review "\ZBfp_mul_mantissa_aux_signs"
% once that is done.
\cs_new:Npn \ZBfp_mul_mantissa:wnnnnnnnn #1\q_stop #2#3#4#5 #6#7#8#9 {
\exp_after:wN \ZBfp_mul_mantissa_aux:w
\tex_the:D \etex_numexpr:D -10000 +
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 99990000 + #2*#6 +
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 99990000 + #2*#7+#3*#6+
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 99990000 + #2*#8+#3*#7+#4*#6+
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 99990000 + #2*#9+#3*#8+#4*#7+#5*#6 +
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 99990000 + #3*#9+#4*#8+#5*#7+
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 99990000 + #4*#9+#5*#8 +
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 100000000 + #5*#9
\ZBfp_mul_mantissa_aux_signs:NN #1
}
\cs_new:Npn \ZBfp_mul_mantissa_aux_signs:NN #1 #2 {
\if:w #1 #2
\exp_after:wN \use_i:nn
\else:
\exp_after:wN \use_ii:nn
\fi:
{ {} \q_stop {+} }
{ {} \q_stop {-} }
}
% We are left with
% "\ZBfp_mul_mantissa_aux:w
% <int>{4d}{4d}{4d}{4d}...{4d}{}\q_stop{<sign>}T{<exp>}"
% where <int> (positive) may range from 0 to 9999. If it is zero,
% we need to remove it and shift the exponent. But we know that
% <int> is in its simplest form, since it comes from a "numexpr"
% calculation, so <int>=0 is equivalent to having its first character
% equal to 0 (I guess "\if_meaning:w" is quicker than "\if_num:w"
% because it involves no expansion).
\cs_new:Npn \ZBfp_mul_mantissa_aux:w #1 #{
\if_meaning:w 0 #1
\exp_after:wN \ZBfp_mul_mantissa_aux_z:nnnnnnww
\else:
\exp_after:wN \ZBfp_mul_mantissa_aux_i:nnnnnww
\fi:
{#1}
}
\cs_new:Npn \ZBfp_mul_mantissa_aux_z:nnnnnnww #1 #2 #3 #4 #5 #6 \q_stop #7 T#8 {
#7 {#2}{#3}{#4}{#5}T{#8}}
\cs_new:Npn \ZBfp_mul_mantissa_aux_i:nnnnnww #1 #2 #3 #4 #5 \q_stop #6 T#7 {
#6 {#1}{#2}{#3}{#4}T{#7+1}}
% ===== Division [not tested properly!]
% Syntax:
% "\ZBfp_div:nn" \marg{ZBfp expr.} \marg{ZBfp expr.}
% First expand the first expression, then swap, and expand the
% second expression. Then treat exponents and place the sign at
% the right place.
\cs_new:Npn \ZBfp_div:nn #1 {
\exp_after:wN \ZBfp_div_aux:Nwn
\tex_romannumeral:D -`\0 #1 \q_stop
}
\cs_new:Npn \ZBfp_div_aux:Nwn #1 \q_stop #2 {
\exp_after:wN \ZBfp_div_check_zero:ww
\tex_romannumeral:D -`\0 #2 \q_stop #1 \q_stop
}
% First let's check if the numerator or denominator is zero.
% Note that "#1" is the denominator!
\cs_new:Npn \ZBfp_div_check_zero:ww #1 \q_stop #2 \q_stop {
\ZBfp_if_zero_aux:NnnnnwNN #1 \use_i_delimit_by_q_stop:nw \use_none:n
{\c_undefined_ZBfp}
\ZBfp_if_zero_aux:NnnnnwNN #2 \use_i_delimit_by_q_stop:nw \use_none:n
{\c_zero_ZBfp}
\ZBfp_div_i:Nn #1 #2
\q_stop
}
% We strip leading zeros from the denominator (specifically <B0>).
\cs_new:Npn \ZBfp_div_i:Nn #1 #2 {
\exp_after:wN \ZBfp_div_ii:NnnnnwNww
\exp_after:wN #1
\exp_after:wN {\tex_number:D #2}
}
% "\ZBfp_div_ii:NnnnnwNww"
% <Bsign> {<trimmed B0>}{<B1>}{<B2>}{<B3>} T {<Bexp>}
% <Asign> {<A0>}{<A1>}{<A2>}{<A3>} T {<Aexp>} \q_stop
% We grab the whole body of <A> as one argument.
\cs_new:Npn \ZBfp_div_ii:NnnnnwNww #1 #2#3#4#5T#6 #7 #8T#9 \q_stop{
\exp_after:wN \ZBfp_div_iii:nnnnnnnnn \tex_number:D #2 \exp_stop_f:
#3 #4 #5 0000\q_stop % <B>, unpacked
#8 % <A>, in blocks of 4
#1 #7 % signs
T {#9 - (#6)} % exponent
Z {\use_none:nnnn #2 000} % extra zeros compensating the power.
}
% Repack <B>.
\cs_new:Npn \ZBfp_div_iii:nnnnnnnnn #1#2#3#4#5#6#7#8#9 {
\ZBfp_div_iiii:nnnnnw {{#1#2#3#4#5}{#6#7#8#9}}
}
\cs_new:Npn \ZBfp_div_iiii:nnnnnw #1 #2#3#4#5 #6 \q_stop {
\ZBfp_div_v:nwnn #1 {#2#3#4#5} \q_stop
}
% Now we have
% "\ZBfp_div_v:nwnn" {<B'0>}{<B'1>}{<B'2>} \q_stop
% {<A0>}{<A1>}{<A2>}{<A3>} <Bsign><Asign> T{<exp>} Z{<zeros>}
%
% where each of <Ai> and <B'i> have 4 digits, except <B'0>, which
% has exactly 5 digits.
%
% We will now use as a quotient
% \[
% <Q0> = \operatorname{Round} \frac{<A0><A1>}{<B'0>+1} - 1.
% \]
% The extra $+1$ and $-1$ are needed to make sure that everything
% stays positive.\footnote{Better ideas are welcome!} The numerator
% of <Q0> is between $10^4$ and $10^8$. The denominator is
% between $10^4$ and $10^5$. So $-1\leq <Q0> \leq 10^4$.
%
\cs_new:Npn \ZBfp_div_v:nwnn #1 #2 \q_stop #3 #4 {
\exp_after:wN \ZBfp_div_vi:nnnnnnn \exp_after:wN
{
\tex_the:D \etex_numexpr:D #3#4 / (#1+1) - 1
}
{#1} #2 {#3 #4}
}
% We calculate $<A> - <Q0>*<B'>$.
%
% "\ZBfp_div_vi:nnnnnnn" {<Q0>} {<B'0>}{<B'1>}{<B'2>}
% {<A0><A1>}{<A2>}{<A3>} <Bsign><Asign> T{<exp>} Z{<zeros>}
%
\cs_new:Npn \ZBfp_div_vi:nnnnnnn #1 #2#3#4 #5#6#7 {
\exp_after:wN \ZBfp_brace_last:nw
\exp_after:wN \ZBfp_div_vii:nnn
\tex_the:D \etex_numexpr:D (#5 - #1*#2)*10000 + #6 - #1*#3 - 20000 +
\exp_after:wN
\ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 200000000 + #7 - #1*#4
{#2}{#3}{#4} {#1}
}
% The term in parentheses is
% \begin{align*}
% <A0><A1> - <Q0>*<B'0>
% &= <A0><A1> - \frac{<B'0>}{<B'0>+1} <A0><A1> + \alpha <B'0>
% \\
% &= \frac{<A0><A1>}{<B'0>+1} + \alpha <B'0>
% \end{align*}
% where $1/2 \leq \alpha \leq 3/2$ (because of e\TeX's rounding
% and our $-1$). This is at most $10^4 + 1.5*<B'0> \leq 1.5\cdot 10^5$.
% Multiplying by $10^4$ won't overflow (but it can be quite close).
%
% All in all, the result should be some long number <s1s2> obeying
% $<s1s2> \leq 15000\cdot<B'0>$ (more or less).
% "\ZBfp_div_vii:nnn" {<s1s2>}{<s3>} {<B'0>}{<B'1>}{<B'2>} {<Q0>}
% <Bsign><Asign> T{<exp>} Z{<zeros>}
%
% We need to be very careful, because <s1s2> can be up to $1.5*10^9$,
% close to overflow. Calculate the next quotient
% \[
% <Q1> = \operatorname{Round} \frac{<s1s2>}{<B'0> + 1} - 1.
% \]
% A few lines above, we get $<s1s2> \leq 15000\cdot <B'0> + 10^8$,
% so $<Q1> \leq 15000 + 10^8/<B'0> \leq 25000$ or so.
%
\cs_new:Npn \ZBfp_div_vii:nnn #1 #2 #3 {
\exp_after:wN \ZBfp_div_viii:nnnnnnn \exp_after:wN
{
\tex_the:D \etex_numexpr:D #1/(#3 + 1) - 1
}
{#1} {#2} {#3}
}
% "\ZBfp_div_viii:nnnnnnn" {<Q1>} {<s1s2>}{<s3>} {<B'0>}{<B'1>}{<B'2>}
% {<Q0>} <Bsign><Asign> T{<exp>} Z{<zeros>}
%
% We now add one level down, to keep enough precision.
% Not sure I'm doing that rigt.
%
\cs_new:Npn \ZBfp_div_viii:nnnnnnn #1 #2#3 #4#5#6 #7 {
\exp_after:wN \ZBfp_brace_last:nw \exp_after:wN \ZBfp_div_ix:nnn
\tex_the:D \etex_numexpr:D (#2 - #1*#4)*10000 + #3 - #1*#5 - 30000 +
\exp_after:wN \ZBfp_brace_last:nnnnnw
\tex_the:D \etex_numexpr:D 300000000 - #1*#6
{#4} {#5} {#7} {#1}
}
% we won't need <B'2> anymore, so we drop "#6".
% "\ZBfp_div_ix:nnn" {<t1t2>}{<t3>} {<B'0>}{<B'1>} {<Q0>}{<Q1>}
% <Bsign><Asign> T{<exp>} Z{<zeros>}
\cs_new:Npn \ZBfp_div_ix:nnn #1#2 #3 {
\exp_after:wN \ZBfp_div_x:nnnnnnn \exp_after:wN
{
\tex_the:D \etex_numexpr:D #1/(#3 + 1) - 1
}
{#1} {#2} {#3}
}
% "\ZBfp_div_x:nnnnnnn" {<Q2>} {<t1t2>}{<t3>} {<B'0>}{<B'1>} {<Q0>}{<Q1>}
% <Bsign><Asign> T{<exp>} Z{<zeros>}
\cs_new:Npn \ZBfp_div_x:nnnnnnn #1 #2#3 #4#5 #6#7 {
\exp_after:wN \ZBfp_brace_last:nw \exp_after:wN \ZBfp_div_xi:nnnnnww
\tex_the:D \etex_numexpr:D (#2 - #1*#4)*10000 + #3 - #1*#5
{#4} {#6}{#7}{#1}
}
% "\ZBfp_div_xi:nnnnn" {<u>} {<B'0>} {<Q0>}{<Q1>}{<Q2>}
% <Bsign><Asign> T{<exp>} Z{<zeros>}
\cs_new:Npn \ZBfp_div_xi:nnnnnww #1 #2 #3 #4 #5 #6 T#7 {
\exp_after:wN \ZBfp_div_xii:nnnnw
\tex_romannumeral:D \etex_numexpr:D -1 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D 9999 + #3 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D 9999 + #4 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D 9999 + #5 + \exp_after:wN \ZBfp_brace_last:nw
\tex_the:D \etex_numexpr:D 10000 + #1 / #2
\exp_after:wN {\tex_the:D \etex_numexpr:D #7\exp_after:wN} #6
}
% "\ZBfp_div_xii:nnnnw" {<Q0>}{<Q1>}{<Q2>}{<Q3>} {<exp>}
% <Bsign><Asign> Z{<zeros>}
% now the <zeros> are explicitly "000" or "00" or "0" or " ",
% and each <Qi> is 4 digits.
\cs_new:Npn \ZBfp_div_xii:nnnnw #1#2#3#4 #5 Z#6 {
\exp_after:wN \ZBfp_div_repack_i:nnnnnnnn #6 #1#2#3#4 0000\q_stop #5
}
% "\ZBfp_div_repack_i:nnnnnnnn" <20-23 digits> \q_stop
\cs_new:Npn \ZBfp_div_repack_i:nnnnnnnn #1#2#3#4 #5#6#7#8 {
\if_num:w #1#2#3#4 = \c_zero
\exp_after:wN \ZBfp_div_repack_ii_small:nnNNNNNNN
\else:
\exp_after:wN \ZBfp_div_repack_ii_large:nnNNNNNNN
\fi:
{#1#2#3#4} {#5#6#7#8}
}
\cs_new:Npn \ZBfp_div_repack_ii_large:nnNNNNNNN #1 #2 #3#4#5#6 #7#8#9 {
\ZBfp_div_repack_iii_large:nnnnNw {#1} {#2} {#3#4#5#6} {#7#8#9}
}
\cs_new:Npn \ZBfp_div_repack_iii_large:nnnnNw #1 #2 #3 #4#5 #6\q_stop #7 #8#9{
\if:w #8#9
\exp_after:wN +
\else:
\exp_after:wN -
\fi:
{#1}{#2}{#3}{#4#5} T{#7}
}
\cs_new:Npn \ZBfp_div_repack_ii_small:nnNNNNNNN #1 #2 #3#4#5#6 #7#8#9 {
\ZBfp_div_repack_iii_small:nnnNNNNNw {#2} {#3#4#5#6} {#7#8#9}
}
\cs_new:Npn \ZBfp_div_repack_iii_small:nnnNNNNNw #1 #2 #3#4 #5#6#7#8 #9\q_stop{
\ZBfp_div_repack_iv_small:nnNN { {#1} {#2} {#3#4} {#5#6#7#8} }
}
% "\ZBfp_div_repack_iv_small:nnNN" {<digits>} {<exp>} <sign><sign>
\cs_new:Npn \ZBfp_div_repack_iv_small:nnNN #1 #2 #3#4 {
\if:w #3#4
\exp_after:wN +
\else:
\exp_after:wN -
\fi:
#1 T{#2-1}
}
% ===== Absolute value
\cs_new:Npn \ZBfp_abs:n #1 {
\exp_after:wN \ZBfp_abs_aux:N \tex_romannumeral:D -`\0 #1
}
\cs_new:Npn \ZBfp_abs_aux:N #1 {+}
% ===== Negate
\cs_new:Npn \ZBfp_neg:n #1 {
\exp_after:wN \ZBfp_neg_aux:N \tex_romannumeral:D -`\0 #1
}
\cs_new:Npn \ZBfp_neg_aux:N #1 {
\if:w - #1
\exp_after:wN +
\else:
\exp_after:wN -
\fi:
}
% ===== Test for zero
% All of the functions below are based on one auxiliary function.
\cs_new:Npn \ZBfp_if_zero_aux:NnnnnwNN #1 #2#3#4#5 T#6 #7#8 {
\if_num:w \etex_numexpr:D #2#3+#4#5 = \c_zero
\exp_after:wN #7
\else:
\exp_after:wN #8
\fi:
}
\prg_new_conditional:Npnn \ZBfp_if_zero:n #1 {p,T,F,TF} {
\ZBfp_if_zero_aux:NnnnnwNN #1
\prg_return_true: \prg_return_false:
}
\prg_new_conditional:Npnn \ZBfp_if_zero:f #1 {p,T,F,TF} {
\exp_after:wN \ZBfp_if_zero_aux:NNNnnnnw
\tex_romannumeral:D -`\0 #1
\prg_return_true: \prg_return_false:
}
\begin{document}
\cs_generate_variant:Nn \iow_term:n {f}
\iow_term:f { \ZBfp_mul:nn
{\ZBfp_add:nn {+{1234}{0000}{0000}{0000}T{2}} {+{0034}{5678}{0000}{0000}T{3}}}
{\ZBfp_add:nn {+{1111}{1111}{1111}{1111}T{0}} {+{2222}{2222}{2222}{2222}T{1}}}
}
\iow_term:f {
\ZBfp_div:nn
{-{0001}{2309}{8124}{3772}T{2}}
{+{2384}{7623}{4098}{1230}T{-2}}
}
\end{document}
\addetc. And I think that an optimized expandable version would beat an optimized assignment-full version, since in the end the calculations are the same. (However, in terms of maintainability, the story is different.) – Bruno Le Floch Feb 20 '11 at 11:19\expandafteruse inl3fp. – Joseph Wright Feb 20 '11 at 11:34\numexpris happy until2e9. – Bruno Le Floch Feb 20 '11 at 11:47fppackage, with a few optimisations to gain some speed. – Joseph Wright Feb 20 '11 at 12:29l3fpand I'll add it to the SVN :-) – Joseph Wright Feb 20 '11 at 12:30xparseif you don't get time before me :). Waiting to see what happens to the code for\::eon LaTeX-L as well. – Bruno Le Floch Feb 20 '11 at 12:49.1, which might reasonably be regarded as valid input. – Joseph Wright Feb 20 '11 at 13:13l3fpspeed testing was done using loops (I think I did 10000 repetitions for sine, deliberated deleting the stored values as part of the loop), with about 100 000 repetitions for division and perhaps 500 000 for addition. These were all of the order of 10-20 s on my PC usingtime. I am quite serious about re-implementing things forl3fp, by the way - as long as it gives the right output, the back-end can be done in a number of ways, andl3fpis still very much new code. – Joseph Wright Feb 21 '11 at 07:20\number\numexpr #1-#5-1+\expandafter\carry:w ~ ~ \number\numexpr 1#2-#6-1+\expandafter\carry:w ~ ~ \number\numexpr 1#3-#7-1+\expandafter\carry:w ~ ~ \number\numexpr 1#4-#8\relaxwhere each argument is 4 digits (and~are just for clarity here).\carry:wputs the last four digits in braces, leaving the carry to be added with the previous level. – Bruno Le Floch Feb 21 '11 at 14:18\fp_add:NNNNNNNNN(internal addition), you only use\c_one_thousand_million. But you just said you used 1+9+ 3 digits? – Bruno Le Floch Feb 21 '11 at 20:21