%
% div16.tex -- perform 16.16 bit fixed point division
% by pts@fazekas.hu at Fri Dec 27 22:40:52 CET 2002 --- ...
%
% This .tex file contains TeX macros for double precision multiplication,
% division of integer and real numbers and image size scaling. The most
% remarkable macro in this file is \@@muldivletposdim (see the comments in
% front of it for detailed documentation).
%
% If you process this file with
% Plain TeX, the unit regression testsuite will be run, and the results of
% the operations issued at the end of this file will be printed to the
% terminal and into the .log file.
%
% The documentation of this macro package is scattered through this file as
% embedded TeX comments.
%
% History
% ~~~~~~~
% started at Fri Dec 27 22:40:52 CET 2002
% --- Sat Dec 28 01:26:52 CET 2002
% 31.31 bit arithmetic: Tue Jan  7 18:12:53 CET 2003
%   skippable, partial \@@muldivpos at Tue Jan  7 21:00:17 CET 2003
%   skippable, finished, doesn't work at Tue Jan  7 21:43:25 CET 2003
%   first successful test case at Tue Jan  7 22:34:32 CET 2003
%   all test cases OK at Tue Jan  7 23:30:32 CET 2003
% Fri Jan 24 12:56:23 CET 2003
%   uses laemu.sty
%

\expandafter\ifx\csname ifLaTeX\endcsname\relax\input laemu.sty\relax\fi%
\NeedsTeXFormat{LaTeX2e}
\ProvidesPackage{div16b}[2003/01/24 v0.0 high precision div and mul]

%
% SUXX: `\count0=2\count1' and `\advance\count0 by 2\count1' disallowed
% Dat:  `\count0=-\count1' and `\advance\count0 by -\count1' OK
%
% Imp: @@ namespace
%

% \def\@empty{} % from laemu.sty
% \catcode`\@=11 % \makeatletter

%\def\loop#1\repeat{%
%  \def\iterate{#1\relax\expandafter\iterate\fi}%
%  \iterate \let\iterate\relax%
%}%
%\let\repeat\fi%

% vvv so we can nest a \@@loop inside \loop or vica versa
\def\@@loop#1\@@repeat{%
  \def\@@iterate{#1\relax\expandafter\@@iterate\fi}%
  \@@iterate \let\@@iterate\relax%
}%
\let\@@repeat\fi%

%\repeat
%alma
%\iftrue
%korte
%\loop

\def\@@gobbleend#1\end{}%

%** `\@@aftergroupbegin ... \end' types `...' token-by-token \aftergroup
%** Argument of \@@aftergroupbegin must not contain \end, {, or }
\def\@@aftergroupbegin#1{%
  \ifx#1\end%
    %\message{done}
  \else%
    %\message{\string#1}%
    \aftergroup#1%
    \expandafter\@@aftergroupbegin%
  \fi%
}%

%** Doesn't make a \global assignment to #1, but keeps its value after the
%** the next `}'.
%** @param #1 name of a \count, \dimen or \skip register (i.e \count0 or
%**   \jot)
\def\@@noresetaftergroup#1{%
  \@@aftergroupbegin#1=\end%
  \expandafter\@@aftergroupbegin\the#1 \end%
  % ^^^ The space before \end is required, because in LaTeX, \end is expandable
  \aftergroup\relax% finish assignment
}%
%** Doesn't work if the definition of #1 contains strange tokens
%** @param #1 a macro \cs name
\def\@@keepdefaftergroup#1{%
  \aftergroup\def%
  \aftergroup#1%
  \aftergroup{%
  \expandafter\@@aftergroupbegin#1\end%
  \aftergroup}%
}%

%** Performs the `#1 := #2' assignment, avoids `Dimension too large' error if
%** both #1 and #2 are \counts.
%** @param #1 a \count or \dimen register specification
%** @param #2 a \count or \dimen register specification, an integer, a real
%**   number or a dimension. Numbers are assumed to have unit `sp'.
\def\@@letdimcount#1#2{{%
  \afterassignment\@@letdimcount@a\count0=#2 \end%
  % Dat: now \@@letdimcount@c contains any of "", "sp", ".45", ".45pt"
  \afterassignment\@@gobbleend\dimen0=1\@@letdimcount@c sp \end% "1", "1sp" "1.45sp" "1.45pt"
  % Dat: 1.99999999sp == 1sp
  \ifdim\dimen0=1sp%
    %\showthe\count0 \message{\the\count0-sp}%
    \afterassignment\@@gobbleend#1=\count0 sp \end%
  \else%
    \afterassignment\@@gobbleend\dimen0=#2 sp \end%
    \afterassignment\@@gobbleend#1=\dimen0 sp \end%
  \fi%
  \@@noresetaftergroup{#1}% return(\T)
}}%
\def\@@letdimcount@a#1\end{\def\@@letdimcount@c{#1}}%
%** Performs the `#1 := #2' assignment, doesn't avoid `Dimension too large'
%** error if both #1 and #2 are \counts.
%** `\@@letdimcount#1#2' is the same as `\@@letdimcountany#1#2{sp}', but
%** \@@letdimcount avoids the error.
%** @param #1 a \count or \dimen register specification
%** @param #2 a \count or \dimen register specification, an integer, a real
%**   number or a dimension. Numbers are assumed to have unit `#3'.
%** @param #3 default metric unit
\def\@@letdimcountany#1#2#3{{%
  \afterassignment\@@letdimcount@a\dimen0=#2 #3 bp\end%
  #1=\dimen0%
  \@@noresetaftergroup{#1}% return(\T)
}}%

\iffalse
% Regression test for \@@letdimcount:
\@@letdimcount{\count0}{2147483647} \message{\the\count0}
\@@letdimcount{\count0}{2147483647sp} \message{\the\count0}
\@@letdimcount{\count0}{2147483647.0sp} \message{\the\count0}
\@@letdimcount{\count0}{2147483647.999999sp} \message{\the\count0}
\@@letdimcount{\count0}{16383.99999pt} \message{\the\count0}
\@@letdimcount{\dimen0}{16383.99999pt} \message{\the\dimen0}
%\@@letdimcount{\count0}{16384pt} % Dimension too large.
%\@@letdimcount{\dimen0}{2000000000} % Dimension too large.
\@@letdimcount{\dimen0}{\count0} \message{\the\dimen0}
\fi

% ---

% !! why no \space in \write?
%** Performs the assignment `#1 := \dimen0 / \dimen1',
%** @param #1 is a name of a TeX dimen/count register (i.e \jot or \count7).
%**   If dimen, then it contains quotient in dimension `pt' (i.e 1in / 1cm
%**   is 2.54pt). If count, then its value divided by 65536 is the quotient
%**   (so this is 16.16 bit fixed-point notation).
%** \dimen0 must not be negative. \dimen1 must be positive.
%** Overwrites \dimen2.
\def\@@divletpos#1{%
  \dimen2=1pt\relax%
  \multiply#1 0% `#1 := 0' works for both \dimen and \count registers
  \loop\ifdim\dimen0>0pt%
    \immediate\write16{0 a=\the\dimen0:b=\the\dimen1:m=\the\dimen2:q=\the#1}
    \@@loop\ifdim\dimen1<0.5\dimen0%
      \advance\dimen1\dimen1%
      \advance\dimen2\dimen2%
      % ^^^ this is the only statement where an arithmetic overflow can occur.
      %     an arithmetic overflow means that the divisor is zero is very
      %     small positive, so the result will be a ``dimen too large''
    \@@repeat%
    \immediate\write16{1 a=\the\dimen0\space b=\the\dimen1\space m=\the\dimen2\space q=\the#1}
    \@@loop\ifdim\dimen1>\dimen0%
      \divide\dimen1 2%
      \divide\dimen2 2%
      \ifdim\dimen2=0pt%
        \dimen0\dimen2%
        \dimen1\dimen2%
      \fi%
    \@@repeat%
    \immediate\write16{2 a=\the\dimen0\space b=\the\dimen1\space m=\the\dimen2\space q=\the#1}
    \ifdim\dimen0>0pt%
      \@@loop\ifdim\dimen0<8192pt%
        \multiply\dimen0 2%
        \multiply\dimen1 2%
      \@@repeat%
    \fi%
    \immediate\write16{3 a=\the\dimen0\space b=\the\dimen1\space m=\the\dimen2\space q=\the#1}
    \advance\dimen0-\dimen1%
    \advance#1\dimen2%
    % ^^^ Dat: because of this, \dimen2 must be a dimen (and not a count)
    \immediate\write16{4 a=\the\dimen0\space b=\the\dimen1\space m=\the\dimen2\space q=\the#1}
  \repeat%
  %\ifdim\dimen0>0pt \message{SUXX}\fi%
}%

%** The same as \@@divletpos, except for #2 and #3 may be either positive or
%** negative, and they might be either dimen or count registers. (A count is
%** measured in sp.)
%** Please don't pass \dimen0, \dimen1 or \dimen2 as argument.
%** @param #1 is a name of a TeX dimen/count register (i.e \jot or \count7).
%**   If dimen, then it contains quotient in dimension `pt' (i.e 1in / 1cm
%**   is 2.54pt). If count, then its value divided by 65536 is the quotient
%**   (so this is 16.16 bit fixed-point notation).
%** @param #2 and #3 are token lists
%**   containing a _positive_ integer or floating point value
%**   (measured in sp, must be specified explicitly), or a dimen such as `2.54cm'.
\def\@@divletdim#1#2#3{{%
  \afterassignment\@@gobbleend\dimen0=#2 sp\end% \dimen0 := #1, convert count -> sp
  \afterassignment\@@gobbleend\dimen1=#3 sp\end%
  \ifdim\dimen1=0pt%
    \divide\dimen0 by0% trigger a TeX division-by-zero error
  \else%
    \ifdim\dimen1<0pt%
      \dimen1-\dimen1%
      \ifdim\dimen0<0pt%
        \dimen0-\dimen0%
        \@@divletpos{#1}%
      \else
        \@@divletpos{#1}%
        #1-#1%
      \fi%
    \else%
      \ifdim\dimen0<0pt%
        \dimen0-\dimen0%
        \@@divletpos{#1}%
        #1-#1%
      \else
        \@@divletpos{#1}%
      \fi%
    \fi%
  \fi%
  % Imp: round to nearest integer
  \@@noresetaftergroup{#1}%
}}%

% --- Dobule (high) precision multiplication and division.

% The type co31.31 denotes two \count registers, representing a nonnegative
% value in 31.31 bit fixed point notation. The sign bits of the registers
% are unused.

%** Subtracts #3#4 from #1#2, modifies #1#2 in place
%** @param #1 \control-sequence on less
%** @param #2 \control-sequence on greater-or-equal
%** @param #3#4 co31.31
%** @param #5#6 co31.31
\def\@@dolt#1#2#3#4#5#6{%
  \ifnum#3<#5%
    \expandafter#1%
  \else%
    \ifnum#3>#5%
      \expandafter#2%
    \else%
      \ifnum#4<#6%
        \expandafter#1%
      \else%
        \expandafter#2%
      \fi%
    \fi%
  \fi%
}%

%** Subtracts #3#4 from #1#2, modifies #1#2 in place
%** @param #1#2 co31.31
%** @param #3#4 co31.31
\def\@@letsub#1#2#3#4#5{%
  \ifnum#2<#4%
    \advance#1-1%
    \advance#2-#4%
    \advance#2-#5% #2 += 0x80000000
  \else%
    \advance#2-#4%
  \fi%
  \advance#1-#3%
}%

%** Adds #3#4 to #1#2, modifies #1#2 in place
%** @param #1#2 co31.31
%** @param #3#4 co31.31
%** @param #5 \countx, whose value is -0x80000000 (-2147483648)
\def\@@letadd#1#2#3#4#5{%
  % \message{add:\the#3:\the#4}%
  \advance#2#5% #2 -= 0x80000000
  \ifnum#2<-#4% orig#2-0x80000000<-#4 <=> orig#2+#4 < 0x80000000
    \advance#2-#5% #2 += 0x80000000
  \else%
    \advance#1 1%
  \fi%
  \advance#2#4%
  \advance#1#3%
}%

%** Multiplies #1#2 by two, modifies #1#2 in place
%** @param #1#2 co31.31
\def\@@lettwice#1#2{%
  \advance#1#1%
  \ifnum#2>1073741823%
    \advance#1 1%
    \advance#2-1073741824%
  \fi%
  \advance#2#2%
}%

%** Divides #1#2 by two, discards LSB, modifies #1#2 in place
%** @param #1#2 co31.31
\def\@@lethalf#1#2{%
  \divide#2 2%
  \ifodd#1%
    \advance#2 1073741824%
  \fi%
  \divide#1 2%
}%

% vvv SUXX: \newif is `\outer'
\newif\ifLOOP
% \def\x{ \expandafter\newif\csname ifLOOP\endcsname}%

%** `\@@muldivletposdim#1#2#3#4' performs the assignment `#1 := #2 * #3 / #4'.
%** `\@@muldivletposdim#1{1pt}#3#4' performs the assignment `#1 := #3 / #4'.
%** `\@@muldivletposdim#1#2#3{1pt}' performs the assignment `#1 := #2 * #3'.
%**
%** All source arguments (#2, #3, #4) can be valid TeX dimen/count values, or
%** \dimen / \count register specifications. Input values may be negative or
%** positive or zero.
%**
%** -- Before each multiplication and division, the arguments are converted
%**    to `pt' (if the argument is a bare number or a \count register, it is
%**    treated as being a dimension in `sp'), the operation is performed on
%**    the numbers, and `pt' is added
%**    to the resulting number, so the result will be a dimension.
%** -- The default unit of measure is `sp'. That is, all bare numbers and
%**    \count registers are treated as being a dimension in `sp'. This
%**    applies to both source and target argument. (You may notice that 1pt
%**    == 65536sp. So the operations are performed in a fixed-point
%**    arithmetic with 16 bits in the fraction.)
%** -- The multiplication is performed exactly (i.e without loss of
%**    precision). However, the division and the 31.31 -> 15.16 rounding
%**    after the division may degrade precision. Internal calculations are
%**    performed in fixed-point arithmetic with 31 bits in the fraction.
%** -- The following constants are valid source arguments, and represent the
%**    same number (one and a half): `1.5pt', `98304sp', `98304 sp',
%**    `98304', `98304.1sp', `98304.9sp'.
%** -- To specify the real number `1.5', please specify `1.5 pt'.
%** -- To specify the integer `42', please specify `42 pt'.
%** -- If you do only division, set #2 := `1pt'. In that special case you may
%**    specify bare integers (or `... sp') for _both_ #3 and #4.
%** -- The following are valid arguments: `\count0',
%**    `\dimen7', `\jot', `\interdisplaylinepenalty'
%** -- Due to a limitation of TeX's fixed-precision arithmetic, it is
%**    impossible to specify exactly `1.5pt' in inches. `0.02075958251953in'
%**    is `1.49974pt', and `0.02075958251954in' is `1.50084pt'. (Of course,
%**    `1.5pt' is `1.5pt'.) That's because TeX converts the real constant
%**    0.02075958251953 to the fraction `1360/65536' and
%**    0.02075958251954 to the fraction `1361/65536', before multiplying.
%** -- If you specify >=2 of \dimen0, \dimen1 and \dimen2 as input, please
%**    specify them in this order.
%**
%** @param #1 target, must be a \dimen or \count register
%**   specification (including `\count0' and `\jot'). If a \dimen register,
%**   its maximum will be 1073741823sp. If a \count register, its maximum
%**   will be 2147483647.
%** @param #2 must be >=0
%** @param #3 must be >=0
%** @param #4 must be >0
%**
\def\@@muldivletposdim#1#2#3#4{{%
  %\afterassignment\@@gobbleend\dimen0=#2 sp\end% \dimen0 := #1, convert count -> sp
  %\afterassignment\@@gobbleend\dimen1=#3 sp\end%
  %\afterassignment\@@gobbleend\dimen2=#4 sp\end%
  \countdef\Y 2 \@@letdimcount\Y{#4}%
  \countdef\F 14 \ifnum\Y<0 \Y-\Y \F-1 \else \F1 \fi%
%  \showthe\Y%
  \ifnum0=\Y % SUXX: \ifnum\Y=0 doesn't work
    \divide\Y0% trigger a TeX division-by-zero error (`arithmetic overflow')
  \else%
    % vvv declare my input arguments as local variables
    \countdef\XA 0 \@@letdimcount\XA{#2}% treat dimen as in `sp'
    \ifnum\XA<0 \XA-\XA \F-\F\fi%
    \countdef\XB 1 \@@letdimcount\XB{#3}%
    \ifnum\XB<0 \XB-\XB \F-\F\fi%
    %\countdef\XA 0 \XA\dimen0% treat dimen as in `sp'
    %\countdef\XB 1 \XB\dimen1%
    %\countdef\Y  2 \Y \dimen2%
    % vvv declare local variables
    \countdef\C  3%
    \countdef\B  4%
    \countdef\R  5%
    \countdef\S  6%
    \countdef\T  7%
    \countdef\U  8%
    \countdef\D  9%
    \countdef\A 10%
    \countdef\L 11%
    \countdef\M 12%
    % vvv ++++ comment out the line below to get rid of the error ++++
    \countdef\G 13 \relax \G-1073741824 \advance\G\G% G := -0x80000000
    %
    % Imp: support negative numbers here
    \ifnum\XB=\Y% shortcut
      \T\XA%
    \else%
      \ifnum\XA=\Y% shortcut
        \T\XB%
      \else% really multiply and divide
        \@@muldivpos%
        %
        % vvv convert (round) \T\U: co31.31 -> \T: dimen15.16
        %     return (\U>=0x7fffc000) ? ((\T+1)<<16) : (\T<<16)+((\U+0x4000)>>15)
        \multiply\T65536% result too large
        \ifnum\U<2147467264% 0x80000000-0x4000
          \advance\U16384%
          \divide\U32768%
          \advance\T\U%
        \else%
          \advance\T65536%
        \fi%
      \fi%
    \fi%
    \ifnum\F<0 \T-\T\fi%
    %\afterassignment\@@gobbleend#1=\T sp\end% works if #1 is a dimen/count
    \@@letdimcount{#1}\T%
    \@@noresetaftergroup{#1}% return(\T)
  \fi% else of division-by-zero
}}%

%** Performs the assignment `\T\U := \XA * \XB / \Y'. Destroys contents of
%** \XA \XB \Y \C \B \R \S \T \U \D \A \L and \M during operations. All
%** control sequences mentioned in this paragraph must be predefined \count
%** registers.
%**   Assumes \newif\LOOP, sets \LOOPfalse upon return
%**   Inputs are: \XA \XB \Y and \G. \XA \XB \Y must be in 15.16 bit
%** fixed-point notation (i.e dimension in `pt' converted to a \count).
%** \G must be -0x80000000.
%**   Outputs are \T and \U. \T\U is a real number in 31.31 bit fixed-point
%** notation (the sign bits are unused).
%**   This is helper routine called by \@@muldivletposdim.
\def\@@muldivpos{%
  %\message{A}%
  \R\XA \divide\R65536 \S-\R \multiply\S65536 \advance\S\XA% my(\R,\S)=(( \XA >>16), \XA &0xffff);
  %\message{A.S=\the\S,R=\the\XA}%
  \T\XB \divide\T65536 \U-\T \multiply\U65536 \advance\U\XB% my(\T,\U)=(( \XB >>16), \XB &0xffff);
  \D0 \A\Y \divide\A65536 \multiply\A32768 % my (\D,\A)=(0, ((( \Y >>16)&0xffff)<<15)); % round() final result
  \L\R \advance\L\L \multiply\L\T \M0 \@@letadd\D\A\L\M\G% \@@letadd( (\D,\A), ((\R*\T)<<1, 0)
  \L\R \multiply\L\U \M\L \divide\L32768 \multiply\L32768 \advance\M-\L \divide\L32768 \multiply\M65536 \@@letadd\D\A\L\M\G% \@@letadd( (\D,\A), ((\R*\U)>>15, ((\R*\U)&0x7fff)<<16)
  \L\S \multiply\L\T \M\L \divide\L32768 \multiply\L32768 \advance\M-\L \divide\L32768 \multiply\M65536 \@@letadd\D\A\L\M\G% \@@letadd( (\D,\A), ((\S*\T)>>15, ((\S*\T)&0x7fff)<<16)
  \ifnum\S>32767 %
    \L0 \M\U \multiply\M32768 \@@letadd\D\A\L\M\G% \@@letadd( (\D,\A), (0, \U<<15)
    \advance\S-32768%
  \fi%
  %\message{B.S=\the\S,U=\the\U}%
  \L0 \M\S \multiply\M\U \@@letadd\D\A\L\M\G% \@@letadd( (\D,\A), (0, \S*\U)
  \C\Y \divide\C 32768 \B-\C \multiply\B32768 \advance\B\Y \multiply\B65536 % my(\C,\B)=(( \Y >>15),(( \Y &0x7fff)<<16));
  \L1 \M0 \T0 \U0% my (\L,\M)=(1,0); (\T,\U)=(0,0);
  %\showthe\D \showthe\A \showthe\C \showthe\B
  %\message{C}%
  \loop% vvv while (\D!=0 or \A!=0) % DA != 0 % \D!=0 <=> \D>0
    %\message{if}%
    \LOOPfalse%
    \ifnum\D>0 \LOOPtrue \fi% Dat: must have a space after `>0'
    \ifnum\A>0 \LOOPtrue \fi%
    \ifLOOP% loop condition
    \R\D \S\A \@@lethalf\R\S% (\R,\S)=(\D,\A); _Half(\R,\S);
    \@@loop% vvv while (_lt((\C,\B),(\R,\S))) % while (CB<DA/2)
      %\message{ina}%
      \@@dolt\LOOPtrue\LOOPfalse\C\B\R\S%
      \ifLOOP% loop condition
      \@@lettwice\C\B%
      \ifnum\L=1073741824% Dat: \L\M is a power of 2
        \divide\L0% die "overflow or division by 0" if \L<0
      \fi%
      \@@lettwice\L\M%
    \@@repeat%
    \@@loop% vvv while (_lt((\D,\A),(\C,\B)))
      %\immediate\write16{D=\the\D:A=\the\A:C=\the\C:B=\the\B}%
      %\message{inb}%
      \@@dolt\LOOPtrue\LOOPfalse\D\A\C\B%
      \ifLOOP% loop condition
      \ifnum\D<1073741824 \@@lettwice\D\A \else \@@lethalf\C\B \fi%
      \@@lethalf\L\M%
      \ifnum0=\L \ifnum0=\M %
        \D0 \A0% (\D,\A)=(0,0) if \L==0 and \M==0;
        \B0 \C0%
      \fi\fi%
    \@@repeat%
    \@@letsub\D\A\C\B\G%
    \@@letadd\T\U\L\M\G%
  \repeat%
  \ifnum\T>32767\divide\T0\fi% die "overflow2 or division by 0" if \T>=0x8000
  % ^^^ !! Imp: more meaningful overflow error message (here and at other places)
  % return value: \T\U: co31.31
  % to convert it to \dimen, do: (\U>=0x7fffc000) ? ((\T+1)<<16) : (\T<<16)+((\U+0x4000)>>15)
}%

\endinput%