OpenCores
URL https://opencores.org/ocsvn/eco32/eco32/trunk

Subversion Repositories eco32

[/] [eco32/] [trunk/] [fp/] [implementation/] [arith/] [mmix-arith.tex] - Blame information for rev 177

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 15 hellwig 
\input cwebmac
2 
% This file is part of the MMIXware package (c) Donald E Knuth 1999
3 
% This material goes at the beginning of all MMIXware CWEB files
4 
 
5 
\def\topofcontents{
6 
  \leftline{\sc\today\ at \hours}\bigskip\bigskip
7 
  \centerline{\titlefont\title}}
8 
 
9 
\font\ninett=cmtt9
10 
\def\botofcontents{\vskip 0pt plus 1filll
11 
    \ninerm\baselineskip10pt
12 
    \noindent\copyright\ 1999 Donald E. Knuth
13 
    \bigskip\noindent
14 
    This file may be freely copied and distributed, provided that
15 
    no changes whatsoever are made. All users are asked to help keep
16 
    the {\ninett MMIX}ware files consistent and ``uncorrupted,''
17 
    identical everywhere in the world. Changes are permissible only
18 
    if the modified file is given a new name, different from the names of
19 
    existing files in the {\ninett MMIX}ware package,
20 
    and only if the modified file is clearly identified
21 
    as not being part of that package.
22 
    (The {\ninett CWEB} system has a ``change file'' facility by
23 
    which users can easily make minor alterations without modifying
24 
    the master source files in any way. Everybody is supposed to use
25 
    change files instead of changing the files.)
26 
    The author has tried his best to produce correct and useful programs,
27 
    in order to help promote computer science research,
28 
    but no warranty of any kind should be assumed.}
29 
 
30 
 
31 
\def\title{MMIX-ARITH}
32 
 
33 
\def\MMIX{\.{MMIX}}
34 
\def\MMIXAL{\.{MMIXAL}}
35 
\def\Hex#1{\hbox{$^{\scriptscriptstyle\#}$\tt#1}} % experimental hex constant
36 
\def\dts{\mathinner{\ldotp\ldotp}}
37 
\def\<#1>{\hbox{$\langle\,$#1$\,\rangle$}}\let\is=\longrightarrow
38 
\def\ff{\\{ff\kern-.05em}}
39 
 
40 
 
41 
 
42 
 
43 
\N{1}{1}Introduction. The subroutines below are used to simulate 64-bit \MMIX\
44 
arithmetic on an old-fashioned 32-bit computer---like the one the author
45 
had when he wrote \MMIXAL\ and the first \MMIX\ simulators in 1998 and 1999.
46 
All operations are fabricated from 32-bit arithmetic, including
47 
a full implementation of the IEEE floating point standard,
48 
assuming only that the \CEE/ compiler has a 32-bit unsigned integer type.
49 
 
50 
Some day 64-bit machines will be commonplace and the awkward manipulations of
51 
the present program will look quite archaic. Interested readers who have such
52 
computers will be able to convert the code to a pure 64-bit form without
53 
difficulty, thereby obtaining much faster and simpler routines. Meanwhile,
54 
however, we can simulate the future and hope for continued progress.
55 
 
56 
This program module has a simple structure, intended to make it
57 
suitable for loading with \MMIX\ simulators and assemblers.
58 
 
59 
\Y\B\8\#\&{include} \.{<stdio.h>}\6
60 
\8\#\&{include} \.{<string.h>}\6
61 
\8\#\&{include} \.{<ctype.h>}\6
62 
\X2:Stuff for \CEE/ preprocessor\X\7
63 
\&{typedef} \&{enum} ${}\{{}$\5
64 
\1${}\\{false},\39{}$\\{true}\5
65 
\2${}\}{}$ \&{bool};\7
66 
\X3:Tetrabyte and octabyte type definitions\X\6
67 
\X36:Other type definitions\X\6
68 
\X4:Global variables\X\6
69 
\X5:Subroutines\X\par
70 
\fi
71 
 
72 
\M{2}Subroutines of this program are declared first with a prototype,
73 
as in {\mc ANSI C}, then with an old-style \CEE/ function definition.
74 
Here are some preprocessor commands that make this work correctly with both
75 
new-style and old-style compilers.
76 
 
77 
\Y\B\4\X2:Stuff for \CEE/ preprocessor\X${}\E{}$\6
78 
\8\#\&{ifdef} \.{\_\_STDC\_\_}\6
79 
\8\#\&{define} \.{ARGS}(\\{list}) \5\\{list}\6
80 
\8\#\&{else}\6
81 
\8\#\&{define} \.{ARGS}(\\{list}) \5(\,)\6
82 
\8\#\&{endif}\par
83 
\U1.\fi
84 
 
85 
\M{3}The definition of type \&{tetra} should be changed, if necessary, so that
86 
it represents an unsigned 32-bit integer.
87 
 
88 
\Y\B\4\X3:Tetrabyte and octabyte type definitions\X${}\E{}$\6
89 
\&{typedef} \&{unsigned} \&{int} \&{tetra};\C{ for systems conforming to the
90 
LP-64 data model }\6
91 
\&{typedef} \&{struct} ${}\{{}$\1\6
92 
\&{tetra} \|h${},{}$ \|l;\2\6
93 
${}\}{}$ \&{octa};\C{ two tetrabytes make one octabyte }\par
94 
\U1.\fi
95 
 
96 
\M{4}\B\D$\\{sign\_bit}$ \5
97 
((\&{unsigned}) \T{\^80000000})\par
98 
\Y\B\4\X4:Global variables\X${}\E{}$\6
99 
\&{octa} \\{zero\_octa};\C{ \PB{$\\{zero\_octa}.\|h\K\\{zero\_octa}.\|l\K%
100 
\T{0}$} }\6
101 
\&{octa} \\{neg\_one}${}\K\{{-}\T{1},\39{-}\T{1}\}{}$;\C{ \PB{$\\{neg\_one}.\|h%
102 
\K\\{neg\_one}.\|l\K{-}\T{1}$} }\6
103 
\&{octa} \\{inf\_octa}${}\K\{\T{\^7ff00000},\39\T{0}\}{}$;\C{ floating point $+%
104 
\infty$ }\6
105 
\&{octa} \\{standard\_NaN}${}\K\{\T{\^7ff80000},\39\T{0}\}{}$;\C{ floating
106 
point NaN(.5) }\6
107 
\&{octa} \\{aux};\C{ auxiliary output of a subroutine }\6
108 
\&{bool} \\{overflow};\C{ set by certain subroutines for signed arithmetic }\par
109 
\As9, 30, 32, 69\ETs75.
110 
\U1.\fi
111 
 
112 
\M{5}It's easy to add and subtract octabytes, if we aren't terribly
113 
worried about speed.
114 
 
115 
\Y\B\4\X5:Subroutines\X${}\E{}$\6
116 
\&{octa} \\{oplus}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
117 
\hbox{}\6{}\&{octa} ${}\\{oplus}(\|y,\39\|z{}$)\C{ compute $y+z$ }\1\1\6
118 
\&{octa} \|y${},{}$ \|z;\2\2\6
119 
${}\{{}$\5
120 
\1\&{octa} \|x;\7
121 
${}\|x.\|h\K\|y.\|h+\|z.\|h{}$;\5
122 
${}\|x.\|l\K\|y.\|l+\|z.\|l;{}$\6
123 
\&{if} ${}(\|x.\|l<\|y.\|l){}$\1\5
124 
${}\|x.\|h\PP;{}$\2\6
125 
\&{return} \|x;\6
126 
\4${}\}{}$\2\7
127 
\&{octa} \\{ominus}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
128 
\hbox{}\6{}\&{octa} ${}\\{ominus}(\|y,\39\|z{}$)\C{ compute $y-z$ }\1\1\6
129 
\&{octa} \|y${},{}$ \|z;\2\2\6
130 
${}\{{}$\5
131 
\1\&{octa} \|x;\7
132 
${}\|x.\|h\K\|y.\|h-\|z.\|h{}$;\5
133 
${}\|x.\|l\K\|y.\|l-\|z.\|l;{}$\6
134 
\&{if} ${}(\|x.\|l>\|y.\|l){}$\1\5
135 
${}\|x.\|h\MM;{}$\2\6
136 
\&{return} \|x;\6
137 
\4${}\}{}$\2\par
138 
\As6, 7, 8, 12, 13, 24, 25, 26, 27, 28, 29, 31, 34, 37, 38, 39, 40, 41, 44, 46,
139 
50, 54, 60, 61, 62, 68, 82, 85, 86, 88, 89, 91\ETs93.
140 
\U1.\fi
141 
 
142 
\M{6}In the following subroutine, \PB{\\{delta}} is a signed quantity that is
143 
assumed to fit in a signed tetrabyte.
144 
 
145 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
146 
\&{octa} \\{incr}\,\,${}\.{ARGS}((\&{octa},\39\&{int})){}$;\5
147 
\hbox{}\6{}\&{octa} ${}\\{incr}(\|y,\39\\{delta}{}$)\C{ compute $y+\delta$ }\1%
148 
\1\6
149 
\&{octa} \|y;\6
150 
\&{int} \\{delta};\2\2\6
151 
${}\{{}$\5
152 
\1\&{octa} \|x;\7
153 
${}\|x.\|h\K\|y.\|h{}$;\5
154 
${}\|x.\|l\K\|y.\|l+\\{delta};{}$\6
155 
\&{if} ${}(\\{delta}\G\T{0}){}$\5
156 
${}\{{}$\1\6
157 
\&{if} ${}(\|x.\|l<\|y.\|l){}$\1\5
158 
${}\|x.\|h\PP;{}$\2\6
159 
\4${}\}{}$\5
160 
\2\&{else} \&{if} ${}(\|x.\|l>\|y.\|l){}$\1\5
161 
${}\|x.\|h\MM;{}$\2\6
162 
\&{return} \|x;\6
163 
\4${}\}{}$\2\par
164 
\fi
165 
 
166 
\M{7}Left and right shifts are only a bit more difficult.
167 
 
168 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
169 
\&{octa} \\{shift\_left}\,\,${}\.{ARGS}((\&{octa},\39\&{int})){}$;\5
170 
\hbox{}\6{}\&{octa} ${}\\{shift\_left}(\|y,\39\|s{}$)\C{ shift left by $s$
171 
bits, where $0\le s\le64$ }\1\1\6
172 
\&{octa} \|y;\6
173 
\&{int} \|s;\2\2\6
174 
${}\{{}$\1\6
175 
\&{while} ${}(\|s\G\T{32}){}$\1\5
176 
${}\|y.\|h\K\|y.\|l,\39\|y.\|l\K\T{0},\39\|s\MRL{-{\K}}\T{32};{}$\2\6
177 
\&{if} (\|s)\5
178 
${}\{{}$\5
179 
\1\&{register} \&{tetra} \\{yhl}${}\K\|y.\|h\LL\|s,{}$ \\{ylh}${}\K\|y.\|l\GG(%
180 
\T{32}-\|s);{}$\7
181 
${}\|y.\|h\K\\{yhl}+\\{ylh}{}$;\5
182 
${}\|y.\|l\MRL{{\LL}{\K}}\|s;{}$\6
183 
\4${}\}{}$\2\6
184 
\&{return} \|y;\6
185 
\4${}\}{}$\2\7
186 
\&{octa} \\{shift\_right}\,\,${}\.{ARGS}((\&{octa},\39\&{int},\39\&{int})){}$;\5
187 
\hbox{}\6{}\&{octa} ${}\\{shift\_right}(\|y,\39\|s,\39\|u{}$)\C{ shift right,
188 
arithmetically if $u=0$ }\1\1\6
189 
\&{octa} \|y;\6
190 
\&{int} \|s${},{}$ \|u;\2\2\6
191 
${}\{{}$\1\6
192 
\&{while} ${}(\|s\G\T{32}){}$\1\5
193 
${}\|y.\|l\K\|y.\|h,\39\|y.\|h\K(\|u\?\T{0}:{-}(\|y.\|h\GG\T{31})),\39\|s%
194 
\MRL{-{\K}}\T{32};{}$\2\6
195 
\&{if} (\|s)\5
196 
${}\{{}$\5
197 
\1\&{register} \&{tetra} \\{yhl}${}\K\|y.\|h\LL(\T{32}-\|s),{}$ \\{ylh}${}\K%
198 
\|y.\|l\GG\|s;{}$\7
199 
${}\|y.\|h\K(\|u\?\T{0}:({-}(\|y.\|h\GG\T{31}))\LL(\T{32}-\|s))+(\|y.\|h\GG%
200 
\|s){}$;\5
201 
${}\|y.\|l\K\\{yhl}+\\{ylh};{}$\6
202 
\4${}\}{}$\2\6
203 
\&{return} \|y;\6
204 
\4${}\}{}$\2\par
205 
\fi
206 
 
207 
\N{1}{8}Multiplication. We need to multiply two unsigned 64-bit integers,
208 
obtaining
209 
an unsigned 128-bit product. It is easy to do this on a 32-bit machine
210 
by using Algorithm 4.3.1M of {\sl Seminumerical Algorithms}, with $b=2^{16}$.
211 
 
212 
The following subroutine returns the lower half of the product, and
213 
puts the upper half into a global octabyte called \PB{\\{aux}}.
214 
 
215 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
216 
\&{octa} \\{omult}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
217 
\hbox{}\6{}\&{octa} ${}\\{omult}(\|y,\39\|z){}$\1\1\6
218 
\&{octa} \|y${},{}$ \|z;\2\2\6
219 
${}\{{}$\1\6
220 
\&{register} \&{int} \|i${},{}$ \|j${},{}$ \|k;\6
221 
\&{tetra} \|u[\T{4}]${},{}$ \|v[\T{4}]${},{}$ \|w[\T{8}];\6
222 
\&{register} \&{tetra} \|t;\6
223 
\&{octa} \\{acc};\7
224 
\X10:Unpack the multiplier and multiplicand to \PB{\|u} and \PB{\|v}\X;\6
225 
\&{for} ${}(\|j\K\T{0};{}$ ${}\|j<\T{4};{}$ ${}\|j\PP){}$\1\5
226 
${}\|w[\|j]\K\T{0};{}$\2\6
227 
\&{for} ${}(\|j\K\T{0};{}$ ${}\|j<\T{4};{}$ ${}\|j\PP){}$\1\6
228 
\&{if} ${}(\R\|v[\|j]){}$\1\5
229 
${}\|w[\|j+\T{4}]\K\T{0};{}$\2\6
230 
\&{else}\5
231 
${}\{{}$\1\6
232 
\&{for} ${}(\|i\K\|k\K\T{0};{}$ ${}\|i<\T{4};{}$ ${}\|i\PP){}$\5
233 
${}\{{}$\1\6
234 
${}\|t\K\|u[\|i]*\|v[\|j]+\|w[\|i+\|j]+\|k;{}$\6
235 
${}\|w[\|i+\|j]\K\|t\AND\T{\^ffff},\39\|k\K\|t\GG\T{16};{}$\6
236 
\4${}\}{}$\2\6
237 
${}\|w[\|j+\T{4}]\K\|k;{}$\6
238 
\4${}\}{}$\2\2\6
239 
\X11:Pack \PB{\|w} into the outputs \PB{\\{aux}} and \PB{\\{acc}}\X;\6
240 
\&{return} \\{acc};\6
241 
\4${}\}{}$\2\par
242 
\fi
243 
 
244 
\M{9}\B\X4:Global variables\X${}\mathrel+\E{}$\6
245 
\&{extern} \&{octa} \\{aux};\C{ secondary output of subroutines with multiple
246 
outputs }\6
247 
\&{extern} \&{bool} \\{overflow};\par
248 
\fi
249 
 
250 
\M{10}\B\X10:Unpack the multiplier and multiplicand to \PB{\|u} and \PB{\|v}%
251 
\X${}\E{}$\6
252 
$\|u[\T{3}]\K\|y.\|h\GG\T{16},\39\|u[\T{2}]\K\|y.\|h\AND\T{\^ffff},\39\|u[%
253 
\T{1}]\K\|y.\|l\GG\T{16},\39\|u[\T{0}]\K\|y.\|l\AND\T{\^ffff};{}$\6
254 
${}\|v[\T{3}]\K\|z.\|h\GG\T{16},\39\|v[\T{2}]\K\|z.\|h\AND\T{\^ffff},\39\|v[%
255 
\T{1}]\K\|z.\|l\GG\T{16},\39\|v[\T{0}]\K\|z.\|l\AND\T{\^ffff}{}$;\par
256 
\U8.\fi
257 
 
258 
\M{11}\B\X11:Pack \PB{\|w} into the outputs \PB{\\{aux}} and \PB{\\{acc}}\X${}%
259 
\E{}$\6
260 
$\\{aux}.\|h\K(\|w[\T{7}]\LL\T{16})+\|w[\T{6}],\39\\{aux}.\|l\K(\|w[\T{5}]\LL%
261 
\T{16})+\|w[\T{4}];{}$\6
262 
${}\\{acc}.\|h\K(\|w[\T{3}]\LL\T{16})+\|w[\T{2}],\39\\{acc}.\|l\K(\|w[\T{1}]\LL%
263 
\T{16})+\|w[\T{0}]{}$;\par
264 
\U8.\fi
265 
 
266 
\M{12}Signed multiplication has the same lower half product as unsigned
267 
multiplication. The signed upper half product is obtained with at most two
268 
further subtractions, after which the result has overflowed if and only if
269 
the upper half is unequal to 64 copies of the sign bit in the lower half.
270 
 
271 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
272 
\&{octa} \\{signed\_omult}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
273 
\hbox{}\6{}\&{octa} ${}\\{signed\_omult}(\|y,\39\|z){}$\1\1\6
274 
\&{octa} \|y${},{}$ \|z;\2\2\6
275 
${}\{{}$\1\6
276 
\&{octa} \\{acc};\7
277 
${}\\{acc}\K\\{omult}(\|y,\39\|z);{}$\6
278 
\&{if} ${}(\|y.\|h\AND\\{sign\_bit}){}$\1\5
279 
${}\\{aux}\K\\{ominus}(\\{aux},\39\|z);{}$\2\6
280 
\&{if} ${}(\|z.\|h\AND\\{sign\_bit}){}$\1\5
281 
${}\\{aux}\K\\{ominus}(\\{aux},\39\|y);{}$\2\6
282 
${}\\{overflow}\K(\\{aux}.\|h\I\\{aux}.\|l\V(\\{aux}.\|h\XOR(\\{aux}.\|h\GG%
283 
\T{1})\XOR(\\{acc}.\|h\AND\\{sign\_bit})));{}$\6
284 
\&{return} \\{acc};\6
285 
\4${}\}{}$\2\par
286 
\fi
287 
 
288 
\N{1}{13}Division. Long division of an unsigned 128-bit integer by an unsigned
289 
64-bit integer is, of course, one of the most challenging routines
290 
needed for \MMIX\ arithmetic. The following program, based on
291 
Algorithm 4.3.1D of {\sl Seminumerical Algorithms}, computes
292 
octabytes $q$ and $r$ such that $(2^{64}x+y)=qz+r$ and $0\le r<z$,
293 
given octabytes $x$, $y$, and~$z$, assuming that $x<z$.
294 
(If $x\ge z$, it simply sets $q=x$ and $r=y$.)
295 
The quotient~$q$ is returned by the subroutine;
296 
the remainder~$r$ is stored in \PB{\\{aux}}.
297 
 
298 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
299 
\&{octa} \\{odiv}\,\,${}\.{ARGS}((\&{octa},\39\&{octa},\39\&{octa})){}$;\5
300 
\hbox{}\6{}\&{octa} ${}\\{odiv}(\|x,\39\|y,\39\|z){}$\1\1\6
301 
\&{octa} \|x${},{}$ \|y${},{}$ \|z;\2\2\6
302 
${}\{{}$\1\6
303 
\&{register} \&{int} \|i${},{}$ \|j${},{}$ \|k${},{}$ \|n${},{}$ \|d;\6
304 
\&{tetra} \|u[\T{8}]${},{}$ \|v[\T{4}]${},{}$ \|q[\T{4}]${},{}$ \\{mask}${},{}$
305 
\\{qhat}${},{}$ \\{rhat}${},{}$ \\{vh}${},{}$ \\{vmh};\6
306 
\&{register} \&{tetra} \|t;\6
307 
\&{octa} \\{acc};\7
308 
\X14:Check that \PB{$\|x<\|z$}; otherwise give trivial answer\X;\6
309 
\X15:Unpack the dividend and divisor to \PB{\|u} and \PB{\|v}\X;\6
310 
\X16:Determine the number of significant places \PB{\|n} in the divisor \PB{%
311 
\|v}\X;\6
312 
\X17:Normalize the divisor\X;\6
313 
\&{for} ${}(\|j\K\T{3};{}$ ${}\|j\G\T{0};{}$ ${}\|j\MM){}$\1\5
314 
\X20:Determine the quotient digit \PB{\|q[\|j]}\X;\2\6
315 
\X18:Unnormalize the remainder\X;\6
316 
\X19:Pack \PB{\|q} and \PB{\|u} to \PB{\\{acc}} and \PB{\\{aux}}\X;\6
317 
\&{return} \\{acc};\6
318 
\4${}\}{}$\2\par
319 
\fi
320 
 
321 
\M{14}\B\X14:Check that \PB{$\|x<\|z$}; otherwise give trivial answer\X${}\E{}$%
322 
\6
323 
\&{if} ${}(\|x.\|h>\|z.\|h\V(\|x.\|h\E\|z.\|h\W\|x.\|l\G\|z.\|l)){}$\5
324 
${}\{{}$\1\6
325 
${}\\{aux}\K\|y{}$;\5
326 
\&{return} \|x;\6
327 
\4${}\}{}$\2\par
328 
\U13.\fi
329 
 
330 
\M{15}\B\X15:Unpack the dividend and divisor to \PB{\|u} and \PB{\|v}\X${}\E{}$%
331 
\6
332 
$\|u[\T{7}]\K\|x.\|h\GG\T{16},\39\|u[\T{6}]\K\|x.\|h\AND\T{\^ffff},\39\|u[%
333 
\T{5}]\K\|x.\|l\GG\T{16},\39\|u[\T{4}]\K\|x.\|l\AND\T{\^ffff};{}$\6
334 
${}\|u[\T{3}]\K\|y.\|h\GG\T{16},\39\|u[\T{2}]\K\|y.\|h\AND\T{\^ffff},\39\|u[%
335 
\T{1}]\K\|y.\|l\GG\T{16},\39\|u[\T{0}]\K\|y.\|l\AND\T{\^ffff};{}$\6
336 
${}\|v[\T{3}]\K\|z.\|h\GG\T{16},\39\|v[\T{2}]\K\|z.\|h\AND\T{\^ffff},\39\|v[%
337 
\T{1}]\K\|z.\|l\GG\T{16},\39\|v[\T{0}]\K\|z.\|l\AND\T{\^ffff}{}$;\par
338 
\U13.\fi
339 
 
340 
\M{16}\B\X16:Determine the number of significant places \PB{\|n} in the divisor
341 
\PB{\|v}\X${}\E{}$\6
342 
\&{for} ${}(\|n\K\T{4};{}$ ${}\|v[\|n-\T{1}]\E\T{0};{}$ ${}\|n\MM){}$\1\5
343 
;\2\par
344 
\U13.\fi
345 
 
346 
\M{17}We shift \PB{\|u} and \PB{\|v} left by \PB{\|d} places, where \PB{\|d} is
347 
chosen to
348 
make $2^{15}\le v_{n-1}<2^{16}$.
349 
 
350 
\Y\B\4\X17:Normalize the divisor\X${}\E{}$\6
351 
$\\{vh}\K\|v[\|n-\T{1}];{}$\6
352 
\&{for} ${}(\|d\K\T{0};{}$ ${}\\{vh}<\T{\^8000};{}$ ${}\|d\PP,\39\\{vh}\MRL{{%
353 
\LL}{\K}}\T{1}){}$\1\5
354 
;\2\6
355 
\&{for} ${}(\|j\K\|k\K\T{0};{}$ ${}\|j<\|n+\T{4};{}$ ${}\|j\PP){}$\5
356 
${}\{{}$\1\6
357 
${}\|t\K(\|u[\|j]\LL\|d)+\|k;{}$\6
358 
${}\|u[\|j]\K\|t\AND\T{\^ffff},\39\|k\K\|t\GG\T{16};{}$\6
359 
\4${}\}{}$\2\6
360 
\&{for} ${}(\|j\K\|k\K\T{0};{}$ ${}\|j<\|n;{}$ ${}\|j\PP){}$\5
361 
${}\{{}$\1\6
362 
${}\|t\K(\|v[\|j]\LL\|d)+\|k;{}$\6
363 
${}\|v[\|j]\K\|t\AND\T{\^ffff},\39\|k\K\|t\GG\T{16};{}$\6
364 
\4${}\}{}$\2\6
365 
${}\\{vh}\K\|v[\|n-\T{1}];{}$\6
366 
${}\\{vmh}\K(\|n>\T{1}\?\|v[\|n-\T{2}]:\T{0}){}$;\par
367 
\U13.\fi
368 
 
369 
\M{18}\B\X18:Unnormalize the remainder\X${}\E{}$\6
370 
$\\{mask}\K(\T{1}\LL\|d)-\T{1};{}$\6
371 
\&{for} ${}(\|j\K\T{3};{}$ ${}\|j\G\|n;{}$ ${}\|j\MM){}$\1\5
372 
${}\|u[\|j]\K\T{0};{}$\2\6
373 
\&{for} ${}(\|k\K\T{0};{}$ ${}\|j\G\T{0};{}$ ${}\|j\MM){}$\5
374 
${}\{{}$\1\6
375 
${}\|t\K(\|k\LL\T{16})+\|u[\|j];{}$\6
376 
${}\|u[\|j]\K\|t\GG\|d,\39\|k\K\|t\AND\\{mask};{}$\6
377 
\4${}\}{}$\2\par
378 
\U13.\fi
379 
 
380 
\M{19}\B\X19:Pack \PB{\|q} and \PB{\|u} to \PB{\\{acc}} and \PB{\\{aux}}\X${}%
381 
\E{}$\6
382 
$\\{acc}.\|h\K(\|q[\T{3}]\LL\T{16})+\|q[\T{2}],\39\\{acc}.\|l\K(\|q[\T{1}]\LL%
383 
\T{16})+\|q[\T{0}];{}$\6
384 
${}\\{aux}.\|h\K(\|u[\T{3}]\LL\T{16})+\|u[\T{2}],\39\\{aux}.\|l\K(\|u[\T{1}]\LL%
385 
\T{16})+\|u[\T{0}]{}$;\par
386 
\U13.\fi
387 
 
388 
\M{20}\B\X20:Determine the quotient digit \PB{\|q[\|j]}\X${}\E{}$\6
389 
${}\{{}$\1\6
390 
\X21:Find the trial quotient, $\hat q$\X;\6
391 
\X22:Subtract $b^j\hat q v$ from \PB{\|u}\X;\6
392 
\X23:If the result was negative, decrease $\hat q$ by 1\X;\6
393 
${}\|q[\|j]\K\\{qhat};{}$\6
394 
\4${}\}{}$\2\par
395 
\U13.\fi
396 
 
397 
\M{21}\B\X21:Find the trial quotient, $\hat q$\X${}\E{}$\6
398 
$\|t\K(\|u[\|j+\|n]\LL\T{16})+\|u[\|j+\|n-\T{1}];{}$\6
399 
${}\\{qhat}\K\|t/\\{vh},\39\\{rhat}\K\|t-\\{vh}*\\{qhat};{}$\6
400 
\&{if} ${}(\|n>\T{1}){}$\1\6
401 
\&{while} ${}(\\{qhat}\E\T{\^10000}\V\\{qhat}*\\{vmh}>(\\{rhat}\LL\T{16})+\|u[%
402 
\|j+\|n-\T{2}]){}$\5
403 
${}\{{}$\1\6
404 
${}\\{qhat}\MM,\39\\{rhat}\MRL{+{\K}}\\{vh};{}$\6
405 
\&{if} ${}(\\{rhat}\G\T{\^10000}){}$\1\5
406 
\&{break};\2\6
407 
\4${}\}{}$\2\2\par
408 
\U20.\fi
409 
 
410 
\M{22}After this step, \PB{$\|u[\|j+\|n]$} will either equal \PB{\|k} or \PB{$%
411 
\|k-\T{1}$}. The
412 
true value of~\PB{\|u} would be obtained by subtracting~\PB{\|k} from \PB{$\|u[%
413 
\|j+\|n]$};
414 
but we don't have to fuss over \PB{$\|u[\|j+\|n]$}, because it won't be
415 
examined later.
416 
 
417 
\Y\B\4\X22:Subtract $b^j\hat q v$ from \PB{\|u}\X${}\E{}$\6
418 
\&{for} ${}(\|i\K\|k\K\T{0};{}$ ${}\|i<\|n;{}$ ${}\|i\PP){}$\5
419 
${}\{{}$\1\6
420 
${}\|t\K\|u[\|i+\|j]+\T{\^ffff0000}-\|k-\\{qhat}*\|v[\|i];{}$\6
421 
${}\|u[\|i+\|j]\K\|t\AND\T{\^ffff},\39\|k\K\T{\^ffff}-(\|t\GG\T{16});{}$\6
422 
\4${}\}{}$\2\par
423 
\U20.\fi
424 
 
425 
\M{23}The correction here occurs only rarely, but it can be necessary---for
426 
example, when dividing the number \Hex{7fff800100000000} by \Hex{800080020005}.
427 
 
428 
\Y\B\4\X23:If the result was negative, decrease $\hat q$ by 1\X${}\E{}$\6
429 
\&{if} ${}(\|u[\|j+\|n]\I\|k){}$\5
430 
${}\{{}$\1\6
431 
${}\\{qhat}\MM;{}$\6
432 
\&{for} ${}(\|i\K\|k\K\T{0};{}$ ${}\|i<\|n;{}$ ${}\|i\PP){}$\5
433 
${}\{{}$\1\6
434 
${}\|t\K\|u[\|i+\|j]+\|v[\|i]+\|k;{}$\6
435 
${}\|u[\|i+\|j]\K\|t\AND\T{\^ffff},\39\|k\K\|t\GG\T{16};{}$\6
436 
\4${}\}{}$\2\6
437 
\4${}\}{}$\2\par
438 
\U20.\fi
439 
 
440 
\M{24}Signed division can be reduced to unsigned division in a tedious
441 
but straightforward manner. We assume that the divisor isn't zero.
442 
 
443 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
444 
\&{octa} \\{signed\_odiv}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
445 
\hbox{}\6{}\&{octa} ${}\\{signed\_odiv}(\|y,\39\|z){}$\1\1\6
446 
\&{octa} \|y${},{}$ \|z;\2\2\6
447 
${}\{{}$\1\6
448 
\&{octa} \\{yy}${},{}$ \\{zz}${},{}$ \|q;\6
449 
\&{register} \&{int} \\{sy}${},{}$ \\{sz};\7
450 
\&{if} ${}(\|y.\|h\AND\\{sign\_bit}){}$\1\5
451 
${}\\{sy}\K\T{2},\39\\{yy}\K\\{ominus}(\\{zero\_octa},\39\|y);{}$\2\6
452 
\&{else}\1\5
453 
${}\\{sy}\K\T{0},\39\\{yy}\K\|y;{}$\2\6
454 
\&{if} ${}(\|z.\|h\AND\\{sign\_bit}){}$\1\5
455 
${}\\{sz}\K\T{1},\39\\{zz}\K\\{ominus}(\\{zero\_octa},\39\|z);{}$\2\6
456 
\&{else}\1\5
457 
${}\\{sz}\K\T{0},\39\\{zz}\K\|z;{}$\2\6
458 
${}\|q\K\\{odiv}(\\{zero\_octa},\39\\{yy},\39\\{zz});{}$\6
459 
${}\\{overflow}\K\\{false};{}$\6
460 
\&{switch} ${}(\\{sy}+\\{sz}){}$\5
461 
${}\{{}$\1\6
462 
\4\&{case} \T{2}${}+\T{1}{}$:\5
463 
${}\\{aux}\K\\{ominus}(\\{zero\_octa},\39\\{aux});{}$\6
464 
\&{if} ${}(\|q.\|h\E\\{sign\_bit}){}$\1\5
465 
${}\\{overflow}\K\\{true};{}$\2\6
466 
\4\&{case} \T{0}${}+\T{0}{}$:\5
467 
\&{return} \|q;\6
468 
\4\&{case} \T{2}${}+\T{0}{}$:\5
469 
\&{if} ${}(\\{aux}.\|h\V\\{aux}.\|l){}$\1\5
470 
${}\\{aux}\K\\{ominus}(\\{zz},\39\\{aux});{}$\2\6
471 
\&{goto} \\{negate\_q};\6
472 
\4\&{case} \T{0}${}+\T{1}{}$:\5
473 
\&{if} ${}(\\{aux}.\|h\V\\{aux}.\|l){}$\1\5
474 
${}\\{aux}\K\\{ominus}(\\{aux},\39\\{zz});{}$\2\6
475 
\4\\{negate\_q}:\5
476 
\&{if} ${}(\\{aux}.\|h\V\\{aux}.\|l){}$\1\5
477 
\&{return} \\{ominus}${}(\\{neg\_one},\39\|q);{}$\2\6
478 
\&{else}\1\5
479 
\&{return} \\{ominus}${}(\\{zero\_octa},\39\|q);{}$\2\6
480 
\4${}\}{}$\2\6
481 
\4${}\}{}$\2\par
482 
\fi
483 
 
484 
\N{1}{25}Bit fiddling. The bitwise operators of \MMIX\ are fairly easy to
485 
implement directly, but three of them occur often enough to deserve
486 
packaging as subroutines.
487 
 
488 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
489 
\&{octa} \\{oand}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
490 
\hbox{}\6{}\&{octa} ${}\\{oand}(\|y,\39\|z{}$)\C{ compute $y\land z$ }\1\1\6
491 
\&{octa} \|y${},{}$ \|z;\2\2\6
492 
${}\{{}$\5
493 
\1\&{octa} \|x;\7
494 
${}\|x.\|h\K\|y.\|h\AND\|z.\|h{}$;\5
495 
${}\|x.\|l\K\|y.\|l\AND\|z.\|l;{}$\6
496 
\&{return} \|x;\6
497 
\4${}\}{}$\2\7
498 
\&{octa} \\{oandn}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
499 
\hbox{}\6{}\&{octa} ${}\\{oandn}(\|y,\39\|z{}$)\C{ compute $y\land\bar z$ }\1\1%
500 
\6
501 
\&{octa} \|y${},{}$ \|z;\2\2\6
502 
${}\{{}$\5
503 
\1\&{octa} \|x;\7
504 
${}\|x.\|h\K\|y.\|h\AND\CM\|z.\|h{}$;\5
505 
${}\|x.\|l\K\|y.\|l\AND\CM\|z.\|l;{}$\6
506 
\&{return} \|x;\6
507 
\4${}\}{}$\2\7
508 
\&{octa} \\{oxor}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
509 
\hbox{}\6{}\&{octa} ${}\\{oxor}(\|y,\39\|z{}$)\C{ compute $y\oplus z$ }\1\1\6
510 
\&{octa} \|y${},{}$ \|z;\2\2\6
511 
${}\{{}$\5
512 
\1\&{octa} \|x;\7
513 
${}\|x.\|h\K\|y.\|h\XOR\|z.\|h{}$;\5
514 
${}\|x.\|l\K\|y.\|l\XOR\|z.\|l;{}$\6
515 
\&{return} \|x;\6
516 
\4${}\}{}$\2\par
517 
\fi
518 
 
519 
\M{26}Here's a fun way to count the number of bits in a tetrabyte.
520 
[This classical trick is called the ``Gillies--Miller method
521 
for sideways addition'' in {\sl The Preparation of Programs
522 
for an Electronic Digital Computer\/} by Wilkes, Wheeler, and
523 
Gill, second edition (Reading, Mass.:\ Addison--Wesley, 1957),
524 
191--193. Some of the tricks used here were suggested by
525 
Balbir Singh, Peter Rossmanith, and Stefan Schwoon.]
526 
 
527 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
528 
\&{int} \\{count\_bits}\,\,\.{ARGS}((\&{tetra}));\5
529 
\hbox{}\6{}\&{int} \\{count\_bits}(\|x)\1\1\6
530 
\&{tetra} \|x;\2\2\6
531 
${}\{{}$\1\6
532 
\&{register} \&{int} \\{xx}${}\K\|x;{}$\7
533 
${}\\{xx}\K\\{xx}-((\\{xx}\GG\T{1})\AND\T{\^55555555});{}$\6
534 
${}\\{xx}\K(\\{xx}\AND\T{\^33333333})+((\\{xx}\GG\T{2})\AND\T{\^33333333});{}$\6
535 
${}\\{xx}\K(\\{xx}+(\\{xx}\GG\T{4}))\AND\T{\^0f0f0f0f};{}$\6
536 
${}\\{xx}\K\\{xx}+(\\{xx}\GG\T{8});{}$\6
537 
\&{return} ${}(\\{xx}+(\\{xx}\GG\T{16}))\AND\T{\^ff};{}$\6
538 
\4${}\}{}$\2\par
539 
\fi
540 
 
541 
\M{27}To compute the nonnegative byte differences of two given tetrabytes,
542 
we can carry out the following 20-step branchless computation:
543 
 
544 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
545 
\&{tetra} \\{byte\_diff}\,\,${}\.{ARGS}((\&{tetra},\39\&{tetra})){}$;\5
546 
\hbox{}\6{}\&{tetra} ${}\\{byte\_diff}(\|y,\39\|z){}$\1\1\6
547 
\&{tetra} \|y${},{}$ \|z;\2\2\6
548 
${}\{{}$\1\6
549 
\&{register} \&{tetra} \|d${}\K(\|y\AND\T{\^00ff00ff})+\T{\^01000100}-(\|z\AND%
550 
\T{\^00ff00ff});{}$\6
551 
\&{register} \&{tetra} \|m${}\K\|d\AND\T{\^01000100};{}$\6
552 
\&{register} \&{tetra} \|x${}\K\|d\AND(\|m-(\|m\GG\T{8}));{}$\7
553 
${}\|d\K((\|y\GG\T{8})\AND\T{\^00ff00ff})+\T{\^01000100}-((\|z\GG\T{8})\AND\T{%
554 
\^00ff00ff});{}$\6
555 
${}\|m\K\|d\AND\T{\^01000100};{}$\6
556 
\&{return} \|x${}+((\|d\AND(\|m-(\|m\GG\T{8})))\LL\T{8});{}$\6
557 
\4${}\}{}$\2\par
558 
\fi
559 
 
560 
\M{28}To compute the nonnegative wyde differences of two tetrabytes,
561 
another trick leads to a 15-step branchless computation.
562 
(Research problem: Can \PB{\\{count\_bits}}, \PB{\\{byte\_diff}}, or \PB{%
563 
\\{wyde\_diff}} be done
564 
with fewer operations?)
565 
 
566 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
567 
\&{tetra} \\{wyde\_diff}\,\,${}\.{ARGS}((\&{tetra},\39\&{tetra})){}$;\5
568 
\hbox{}\6{}\&{tetra} ${}\\{wyde\_diff}(\|y,\39\|z){}$\1\1\6
569 
\&{tetra} \|y${},{}$ \|z;\2\2\6
570 
${}\{{}$\1\6
571 
\&{register} \&{tetra} \|a${}\K((\|y\GG\T{16})-(\|z\GG\T{16}))\AND\T{%
572 
\^10000};{}$\6
573 
\&{register} \&{tetra} \|b${}\K((\|y\AND\T{\^ffff})-(\|z\AND\T{\^ffff}))\AND\T{%
574 
\^10000};{}$\7
575 
\&{return} \|y${}-(\|z\XOR((\|y\XOR\|z)\AND(\|b-\|a-(\|b\GG\T{16}))));{}$\6
576 
\4${}\}{}$\2\par
577 
\fi
578 
 
579 
\M{29}The last bitwise subroutine we need is the most interesting:
580 
It implements \MMIX's \.{MOR} and \.{MXOR} operations.
581 
 
582 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
583 
\&{octa} \\{bool\_mult}\,\,${}\.{ARGS}((\&{octa},\39\&{octa},\39\&{bool})){}$;\5
584 
\hbox{}\6{}\&{octa} ${}\\{bool\_mult}(\|y,\39\|z,\39\\{xor}){}$\1\1\6
585 
\&{octa} \|y${},{}$ \|z;\C{ the operands }\6
586 
\&{bool} \\{xor};\C{ do we do xor instead of or? }\2\2\6
587 
${}\{{}$\1\6
588 
\&{octa} \|o${},{}$ \|x;\6
589 
\&{register} \&{tetra} \|a${},{}$ \|b${},{}$ \|c;\6
590 
\&{register} \&{int} \|k;\7
591 
\&{for} ${}(\|k\K\T{0},\39\|o\K\|y,\39\|x\K\\{zero\_octa};{}$ ${}\|o.\|h\V\|o.%
592 
\|l;{}$ ${}\|k\PP,\39\|o\K\\{shift\_right}(\|o,\39\T{8},\39\T{1})){}$\1\6
593 
\&{if} ${}(\|o.\|l\AND\T{\^ff}){}$\5
594 
${}\{{}$\1\6
595 
${}\|a\K((\|z.\|h\GG\|k)\AND\T{\^01010101})*\T{\^ff};{}$\6
596 
${}\|b\K((\|z.\|l\GG\|k)\AND\T{\^01010101})*\T{\^ff};{}$\6
597 
${}\|c\K(\|o.\|l\AND\T{\^ff})*\T{\^01010101};{}$\6
598 
\&{if} (\\{xor})\1\5
599 
${}\|x.\|h\MRL{{\XOR}{\K}}\|a\AND\|c,\39\|x.\|l\MRL{{\XOR}{\K}}\|b\AND\|c;{}$\2%
600 
\6
601 
\&{else}\1\5
602 
${}\|x.\|h\MRL{{\OR}{\K}}\|a\AND\|c,\39\|x.\|l\MRL{{\OR}{\K}}\|b\AND\|c;{}$\2\6
603 
\4${}\}{}$\2\2\6
604 
\&{return} \|x;\6
605 
\4${}\}{}$\2\par
606 
\fi
607 
 
608 
\N{1}{30}Floating point packing and unpacking. Standard IEEE floating binary
609 
numbers pack a sign, exponent, and fraction into a tetrabyte
610 
or octabyte. In this section we consider basic subroutines that
611 
convert between IEEE format and the separate unpacked components.
612 
 
613 
\Y\B\4\D$\.{ROUND\_OFF}$ \5
614 
\T{1}\par
615 
\B\4\D$\.{ROUND\_UP}$ \5
616 
\T{2}\par
617 
\B\4\D$\.{ROUND\_DOWN}$ \5
618 
\T{3}\par
619 
\B\4\D$\.{ROUND\_NEAR}$ \5
620 
\T{4}\par
621 
\Y\B\4\X4:Global variables\X${}\mathrel+\E{}$\6
622 
\&{int} \\{cur\_round};\C{ the current rounding mode }\par
623 
\fi
624 
 
625 
\M{31}The \PB{\\{fpack}} routine takes an octabyte $f$, a raw exponent~$e$,
626 
and a sign~\PB{\|s}, and packs them
627 
into the floating binary number that corresponds to
628 
$\pm2^{e-1076}f$, using a given rounding mode.
629 
The value of $f$ should satisfy $2^{54}\le f\le 2^{55}$.
630 
 
631 
Thus, for example, the floating binary number $+1.0=\Hex{3ff0000000000000}$
632 
is obtained when $f=2^{54}$, $e=\Hex{3fe}$, and \PB{$\|s\K\.{'+'}$}.
633 
The raw exponent~$e$ is usually one less than
634 
the final exponent value; the leading bit of~$f$ is essentially added
635 
to the exponent. (This trick works nicely for subnormal numbers, when
636 
$e<0$, or in cases where the value of $f$ is rounded upwards to $2^{55}$.)
637 
 
638 
Exceptional events are noted by oring appropriate bits into
639 
the global variable \PB{\\{exceptions}}. Special considerations apply to
640 
underflow, which is not fully specified by Section 7.4 of the IEEE standard:
641 
Implementations of the standard are free to choose between two definitions
642 
of ``tininess'' and two definitions of ``accuracy loss.''
643 
\MMIX\ determines tininess {\it after\/} rounding, hence a result with
644 
$e<0$ is not necessarily tiny; \MMIX\ treats accuracy loss as equivalent
645 
to inexactness. Thus, a result underflows if and only if
646 
it is tiny and either (i)~it is inexact or (ii)~the underflow trap is enabled.
647 
The \PB{\\{fpack}} routine sets \PB{\.{U\_BIT}} in \PB{\\{exceptions}} if and
648 
only if the result is
649 
tiny, \PB{\.{X\_BIT}} if and only if the result is inexact.
650 
 
651 
\Y\B\4\D$\.{X\_BIT}$ \5
652 
$(\T{1}\LL\T{8}{}$)\C{ floating inexact }\par
653 
\B\4\D$\.{Z\_BIT}$ \5
654 
$(\T{1}\LL\T{9}{}$)\C{ floating division by zero }\par
655 
\B\4\D$\.{U\_BIT}$ \5
656 
$(\T{1}\LL\T{10}{}$)\C{ floating underflow }\par
657 
\B\4\D$\.{O\_BIT}$ \5
658 
$(\T{1}\LL\T{11}{}$)\C{ floating overflow }\par
659 
\B\4\D$\.{I\_BIT}$ \5
660 
$(\T{1}\LL\T{12}{}$)\C{ floating invalid operation }\par
661 
\B\4\D$\.{W\_BIT}$ \5
662 
$(\T{1}\LL\T{13}{}$)\C{ float-to-fix overflow }\par
663 
\B\4\D$\.{V\_BIT}$ \5
664 
$(\T{1}\LL\T{14}{}$)\C{ integer overflow }\par
665 
\B\4\D$\.{D\_BIT}$ \5
666 
$(\T{1}\LL\T{15}{}$)\C{ integer divide check }\par
667 
\B\4\D$\.{E\_BIT}$ \5
668 
$(\T{1}\LL\T{18}{}$)\C{ external (dynamic) trap bit }\par
669 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
670 
\&{octa} \\{fpack}\,\,${}\.{ARGS}((\&{octa},\39\&{int},\39\&{char},\39%
671 
\&{int})){}$;\5
672 
\hbox{}\6{}\&{octa} ${}\\{fpack}(\|f,\39\|e,\39\|s,\39\|r){}$\1\1\6
673 
\&{octa} \|f;\C{ the normalized fraction part }\6
674 
\&{int} \|e;\C{ the raw exponent }\6
675 
\&{char} \|s;\C{ the sign }\6
676 
\&{int} \|r;\C{ the rounding mode }\2\2\6
677 
${}\{{}$\1\6
678 
\&{octa} \|o;\7
679 
\&{if} ${}(\|e>\T{\^7fd}){}$\1\5
680 
${}\|e\K\T{\^7ff},\39\|o\K\\{zero\_octa};{}$\2\6
681 
\&{else}\5
682 
${}\{{}$\1\6
683 
\&{if} ${}(\|e<\T{0}){}$\5
684 
${}\{{}$\1\6
685 
\&{if} ${}(\|e<{-}\T{54}){}$\1\5
686 
${}\|o.\|h\K\T{0},\39\|o.\|l\K\T{1};{}$\2\6
687 
\&{else}\5
688 
${}\{{}$\5
689 
\1\&{octa} \\{oo};\7
690 
${}\|o\K\\{shift\_right}(\|f,\39{-}\|e,\39\T{1});{}$\6
691 
${}\\{oo}\K\\{shift\_left}(\|o,\39{-}\|e);{}$\6
692 
\&{if} ${}(\\{oo}.\|l\I\|f.\|l\V\\{oo}.\|h\I\|f.\|h){}$\1\5
693 
${}\|o.\|l\MRL{{\OR}{\K}}\T{1}{}$;\C{ sticky bit }\2\6
694 
\4${}\}{}$\2\6
695 
${}\|e\K\T{0};{}$\6
696 
\4${}\}{}$\5
697 
\2\&{else}\1\5
698 
${}\|o\K\|f;{}$\2\6
699 
\4${}\}{}$\2\6
700 
\X33:Round and return the result\X;\6
701 
\4${}\}{}$\2\par
702 
\fi
703 
 
704 
\M{32}\B\X4:Global variables\X${}\mathrel+\E{}$\6
705 
\&{int} \\{exceptions};\C{ bits possibly destined for rA }\par
706 
\fi
707 
 
708 
\M{33}Everything falls together so nicely here, it's almost too good to be
709 
true!
710 
 
711 
\Y\B\4\X33:Round and return the result\X${}\E{}$\6
712 
\&{if} ${}(\|o.\|l\AND\T{3}){}$\1\5
713 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{X\_BIT};{}$\2\6
714 
\&{switch} (\|r)\5
715 
${}\{{}$\1\6
716 
\4\&{case} \.{ROUND\_DOWN}:\5
717 
\&{if} ${}(\|s\E\.{'-'}){}$\1\5
718 
${}\|o\K\\{incr}(\|o,\39\T{3}){}$;\5
719 
\2\&{break};\6
720 
\4\&{case} \.{ROUND\_UP}:\5
721 
\&{if} ${}(\|s\I\.{'-'}){}$\1\5
722 
${}\|o\K\\{incr}(\|o,\39\T{3});{}$\2\6
723 
\4\&{case} \.{ROUND\_OFF}:\5
724 
\&{break};\6
725 
\4\&{case} \.{ROUND\_NEAR}:\5
726 
${}\|o\K\\{incr}(\|o,\39\|o.\|l\AND\T{4}\?\T{2}:\T{1}){}$;\5
727 
\&{break};\6
728 
\4${}\}{}$\2\6
729 
${}\|o\K\\{shift\_right}(\|o,\39\T{2},\39\T{1});{}$\6
730 
${}\|o.\|h\MRL{+{\K}}\|e\LL\T{20};{}$\6
731 
\&{if} ${}(\|o.\|h\G\T{\^7ff00000}){}$\1\5
732 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{O\_BIT}+\.{X\_BIT}{}$;\C{ overflow }\2\6
733 
\&{else} \&{if} ${}(\|o.\|h<\T{\^100000}){}$\1\5
734 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{U\_BIT}{}$;\C{ tininess }\2\6
735 
\&{if} ${}(\|s\E\.{'-'}){}$\1\5
736 
${}\|o.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
737 
\&{return} \|o;\par
738 
\U31.\fi
739 
 
740 
\M{34}Similarly, \PB{\\{sfpack}} packs a short float, from inputs
741 
having the same conventions as \PB{\\{fpack}}.
742 
 
743 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
744 
\&{tetra} \\{sfpack}\,\,${}\.{ARGS}((\&{octa},\39\&{int},\39\&{char},\39%
745 
\&{int})){}$;\5
746 
\hbox{}\6{}\&{tetra} ${}\\{sfpack}(\|f,\39\|e,\39\|s,\39\|r){}$\1\1\6
747 
\&{octa} \|f;\C{ the fraction part }\6
748 
\&{int} \|e;\C{ the raw exponent }\6
749 
\&{char} \|s;\C{ the sign }\6
750 
\&{int} \|r;\C{ the rounding mode }\2\2\6
751 
${}\{{}$\1\6
752 
\&{register} \&{tetra} \|o;\7
753 
\&{if} ${}(\|e>\T{\^47d}){}$\1\5
754 
${}\|e\K\T{\^47f},\39\|o\K\T{0};{}$\2\6
755 
\&{else}\5
756 
${}\{{}$\1\6
757 
${}\|o\K\\{shift\_left}(\|f,\39\T{3}).\|h;{}$\6
758 
\&{if} ${}(\|f.\|l\AND\T{\^1fffffff}){}$\1\5
759 
${}\|o\MRL{{\OR}{\K}}\T{1};{}$\2\6
760 
\&{if} ${}(\|e<\T{\^380}){}$\5
761 
${}\{{}$\1\6
762 
\&{if} ${}(\|e<\T{\^380}-\T{25}){}$\1\5
763 
${}\|o\K\T{1};{}$\2\6
764 
\&{else}\5
765 
${}\{{}$\5
766 
\1\&{register} \&{tetra} \\{o0}${},{}$ \\{oo};\7
767 
${}\\{o0}\K\|o;{}$\6
768 
${}\|o\K\|o\GG(\T{\^380}-\|e);{}$\6
769 
${}\\{oo}\K\|o\LL(\T{\^380}-\|e);{}$\6
770 
\&{if} ${}(\\{oo}\I\\{o0}){}$\1\5
771 
${}\|o\MRL{{\OR}{\K}}\T{1}{}$;\C{ sticky bit }\2\6
772 
\4${}\}{}$\2\6
773 
${}\|e\K\T{\^380};{}$\6
774 
\4${}\}{}$\2\6
775 
\4${}\}{}$\2\6
776 
\X35:Round and return the short result\X;\6
777 
\4${}\}{}$\2\par
778 
\fi
779 
 
780 
\M{35}\B\X35:Round and return the short result\X${}\E{}$\6
781 
\&{if} ${}(\|o\AND\T{3}){}$\1\5
782 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{X\_BIT};{}$\2\6
783 
\&{switch} (\|r)\5
784 
${}\{{}$\1\6
785 
\4\&{case} \.{ROUND\_DOWN}:\5
786 
\&{if} ${}(\|s\E\.{'-'}){}$\1\5
787 
${}\|o\MRL{+{\K}}\T{3}{}$;\5
788 
\2\&{break};\6
789 
\4\&{case} \.{ROUND\_UP}:\5
790 
\&{if} ${}(\|s\I\.{'-'}){}$\1\5
791 
${}\|o\MRL{+{\K}}\T{3};{}$\2\6
792 
\4\&{case} \.{ROUND\_OFF}:\5
793 
\&{break};\6
794 
\4\&{case} \.{ROUND\_NEAR}:\5
795 
${}\|o\MRL{+{\K}}(\|o\AND\T{4}\?\T{2}:\T{1}){}$;\5
796 
\&{break};\6
797 
\4${}\}{}$\2\6
798 
${}\|o\K\|o\GG\T{2};{}$\6
799 
${}\|o\MRL{+{\K}}(\|e-\T{\^380})\LL\T{23};{}$\6
800 
\&{if} ${}(\|o\G\T{\^7f800000}){}$\1\5
801 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{O\_BIT}+\.{X\_BIT}{}$;\C{ overflow }\2\6
802 
\&{else} \&{if} ${}(\|o<\T{\^100000}){}$\1\5
803 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{U\_BIT}{}$;\C{ tininess }\2\6
804 
\&{if} ${}(\|s\E\.{'-'}){}$\1\5
805 
${}\|o\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
806 
\&{return} \|o;\par
807 
\U34.\fi
808 
 
809 
\M{36}The \PB{\\{funpack}} routine is, roughly speaking, the opposite of \PB{%
810 
\\{fpack}}.
811 
It takes a given floating point number~$x$ and separates out its
812 
fraction part~$f$, exponent~$e$, and sign~$s$. It clears \PB{\\{exceptions}}
813 
to zero. It returns the type of value found: \PB{\\{zro}}, \PB{\\{num}}, \PB{%
814 
\\{inf}},
815 
or \PB{\\{nan}}. When it returns \PB{\\{num}},
816 
it will have set $f$, $e$, and~$s$
817 
to the values from which \PB{\\{fpack}} would produce the original number~$x$
818 
without exceptions.
819 
 
820 
\Y\B\4\D$\\{zero\_exponent}$ \5
821 
$({-}\T{1000}{}$)\C{ zero is assumed to have this exponent }\par
822 
\Y\B\4\X36:Other type definitions\X${}\E{}$\6
823 
\&{typedef} \&{enum} ${}\{{}$\1\6
824 
${}\\{zro},\39\\{num},\39\\{inf},\39\\{nan}{}$\2\6
825 
${}\}{}$\5
826 
\&{ftype};\par
827 
\A59.
828 
\U1.\fi
829 
 
830 
\M{37}\B\X5:Subroutines\X${}\mathrel+\E{}$\6
831 
\&{ftype} \\{funpack}\,\,${}\.{ARGS}((\&{octa},\39{}$\&{octa} ${}{*},\39{}$%
832 
\&{int} ${}{*},\39{}$\&{char} ${}{*})){}$;\5
833 
\hbox{}\6{}\&{ftype} ${}\\{funpack}(\|x,\39\|f,\39\|e,\39\|s){}$\1\1\6
834 
\&{octa} \|x;\C{ the given floating point value }\6
835 
\&{octa} ${}{*}\|f{}$;\C{ address where the fraction part should be stored }\6
836 
\&{int} ${}{*}\|e{}$;\C{ address where the exponent part should be stored }\6
837 
\&{char} ${}{*}\|s{}$;\C{ address where the sign should be stored }\2\2\6
838 
${}\{{}$\1\6
839 
\&{register} \&{int} \\{ee};\7
840 
${}\\{exceptions}\K\T{0};{}$\6
841 
${}{*}\|s\K(\|x.\|h\AND\\{sign\_bit}\?\.{'-'}:\.{'+'});{}$\6
842 
${}{*}\|f\K\\{shift\_left}(\|x,\39\T{2});{}$\6
843 
${}\|f\MG\|h\MRL{\AND{\K}}\T{\^3fffff};{}$\6
844 
${}\\{ee}\K(\|x.\|h\GG\T{20})\AND\T{\^7ff};{}$\6
845 
\&{if} (\\{ee})\5
846 
${}\{{}$\1\6
847 
${}{*}\|e\K\\{ee}-\T{1};{}$\6
848 
${}\|f\MG\|h\MRL{{\OR}{\K}}\T{\^400000};{}$\6
849 
\&{return} ${}(\\{ee}<\T{\^7ff}\?\\{num}:\|f\MG\|h\E\T{\^400000}\W\R\|f\MG\|l\?%
850 
\\{inf}:\\{nan});{}$\6
851 
\4${}\}{}$\2\6
852 
\&{if} ${}(\R\|x.\|l\W\R\|f\MG\|h){}$\5
853 
${}\{{}$\1\6
854 
${}{*}\|e\K\\{zero\_exponent}{}$;\5
855 
\&{return} \\{zro};\6
856 
\4${}\}{}$\2\6
857 
\&{do}\5
858 
${}\{{}$\5
859 
\1${}\\{ee}\MM{}$;\5
860 
${}{*}\|f\K\\{shift\_left}({*}\|f,\39\T{1}){}$;\5
861 
${}\}{}$\2\5
862 
\&{while} ${}(\R(\|f\MG\|h\AND\T{\^400000}));{}$\6
863 
${}{*}\|e\K\\{ee}{}$;\5
864 
\&{return} \\{num};\6
865 
\4${}\}{}$\2\par
866 
\fi
867 
 
868 
\M{38}\B\X5:Subroutines\X${}\mathrel+\E{}$\6
869 
\&{ftype} \\{sfunpack}\,\,${}\.{ARGS}((\&{tetra},\39{}$\&{octa} ${}{*},\39{}$%
870 
\&{int} ${}{*},\39{}$\&{char} ${}{*})){}$;\5
871 
\hbox{}\6{}\&{ftype} ${}\\{sfunpack}(\|x,\39\|f,\39\|e,\39\|s){}$\1\1\6
872 
\&{tetra} \|x;\C{ the given floating point value }\6
873 
\&{octa} ${}{*}\|f{}$;\C{ address where the fraction part should be stored }\6
874 
\&{int} ${}{*}\|e{}$;\C{ address where the exponent part should be stored }\6
875 
\&{char} ${}{*}\|s{}$;\C{ address where the sign should be stored }\2\2\6
876 
${}\{{}$\1\6
877 
\&{register} \&{int} \\{ee};\7
878 
${}\\{exceptions}\K\T{0};{}$\6
879 
${}{*}\|s\K(\|x\AND\\{sign\_bit}\?\.{'-'}:\.{'+'});{}$\6
880 
${}\|f\MG\|h\K(\|x\GG\T{1})\AND\T{\^3fffff},\39\|f\MG\|l\K\|x\LL\T{31};{}$\6
881 
${}\\{ee}\K(\|x\GG\T{23})\AND\T{\^ff};{}$\6
882 
\&{if} (\\{ee})\5
883 
${}\{{}$\1\6
884 
${}{*}\|e\K\\{ee}+\T{\^380}-\T{1};{}$\6
885 
${}\|f\MG\|h\MRL{{\OR}{\K}}\T{\^400000};{}$\6
886 
\&{return} ${}(\\{ee}<\T{\^ff}\?\\{num}:(\|x\AND\T{\^7fffffff})\E\T{\^7f800000}%
887 
\?\\{inf}:\\{nan});{}$\6
888 
\4${}\}{}$\2\6
889 
\&{if} ${}(\R(\|x\AND\T{\^7fffffff})){}$\5
890 
${}\{{}$\1\6
891 
${}{*}\|e\K\\{zero\_exponent}{}$;\5
892 
\&{return} \\{zro};\6
893 
\4${}\}{}$\2\6
894 
\&{do}\5
895 
${}\{{}$\5
896 
\1${}\\{ee}\MM{}$;\5
897 
${}{*}\|f\K\\{shift\_left}({*}\|f,\39\T{1}){}$;\5
898 
${}\}{}$\2\5
899 
\&{while} ${}(\R(\|f\MG\|h\AND\T{\^400000}));{}$\6
900 
${}{*}\|e\K\\{ee}+\T{\^380}{}$;\5
901 
\&{return} \\{num};\6
902 
\4${}\}{}$\2\par
903 
\fi
904 
 
905 
\M{39}Since \MMIX\ downplays 32-bit operations, it uses \PB{\\{sfpack}} and %
906 
\PB{\\{sfunpack}}
907 
only when loading and storing short floats, or when converting
908 
from fixed point to floating point.
909 
 
910 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
911 
\&{octa} \\{load\_sf}\,\,\.{ARGS}((\&{tetra}));\5
912 
\hbox{}\6{}\&{octa} \\{load\_sf}(\|z)\1\1\6
913 
\&{tetra} \|z;\C{ 32 bits to be loaded into a 64-bit register }\2\2\6
914 
${}\{{}$\1\6
915 
\&{octa} \|f${},{}$ \|x;\5
916 
\&{int} \|e;\5
917 
\&{char} \|s;\5
918 
\&{ftype} \|t;\7
919 
${}\|t\K\\{sfunpack}(\|z,\39{\AND}\|f,\39{\AND}\|e,\39{\AND}\|s);{}$\6
920 
\&{switch} (\|t)\5
921 
${}\{{}$\1\6
922 
\4\&{case} \\{zro}:\5
923 
${}\|x\K\\{zero\_octa}{}$;\5
924 
\&{break};\6
925 
\4\&{case} \\{num}:\5
926 
\&{return} \\{fpack}${}(\|f,\39\|e,\39\|s,\39\.{ROUND\_OFF});{}$\6
927 
\4\&{case} \\{inf}:\5
928 
${}\|x\K\\{inf\_octa}{}$;\5
929 
\&{break};\6
930 
\4\&{case} \\{nan}:\5
931 
${}\|x\K\\{shift\_right}(\|f,\39\T{2},\39\T{1}){}$;\5
932 
${}\|x.\|h\MRL{{\OR}{\K}}\T{\^7ff00000}{}$;\5
933 
\&{break};\6
934 
\4${}\}{}$\2\6
935 
\&{if} ${}(\|s\E\.{'-'}){}$\1\5
936 
${}\|x.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
937 
\&{return} \|x;\6
938 
\4${}\}{}$\2\par
939 
\fi
940 
 
941 
\M{40}\B\X5:Subroutines\X${}\mathrel+\E{}$\6
942 
\&{tetra} \\{store\_sf}\,\,\.{ARGS}((\&{octa}));\5
943 
\hbox{}\6{}\&{tetra} \\{store\_sf}(\|x)\1\1\6
944 
\&{octa} \|x;\C{ 64 bits to be loaded into a 32-bit word }\2\2\6
945 
${}\{{}$\1\6
946 
\&{octa} \|f;\5
947 
\&{tetra} \|z;\5
948 
\&{int} \|e;\5
949 
\&{char} \|s;\5
950 
\&{ftype} \|t;\7
951 
${}\|t\K\\{funpack}(\|x,\39{\AND}\|f,\39{\AND}\|e,\39{\AND}\|s);{}$\6
952 
\&{switch} (\|t)\5
953 
${}\{{}$\1\6
954 
\4\&{case} \\{zro}:\5
955 
${}\|z\K\T{0}{}$;\5
956 
\&{break};\6
957 
\4\&{case} \\{num}:\5
958 
\&{return} \\{sfpack}${}(\|f,\39\|e,\39\|s,\39\\{cur\_round});{}$\6
959 
\4\&{case} \\{inf}:\5
960 
${}\|z\K\T{\^7f800000}{}$;\5
961 
\&{break};\6
962 
\4\&{case} \\{nan}:\5
963 
\&{if} ${}(\R(\|f.\|h\AND\T{\^200000})){}$\5
964 
${}\{{}$\1\6
965 
${}\|f.\|h\MRL{{\OR}{\K}}\T{\^200000}{}$;\5
966 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\C{ NaN was signaling }\6
967 
\4${}\}{}$\2\6
968 
${}\|z\K\T{\^7f800000}\OR(\|f.\|h\LL\T{1})\OR(\|f.\|l\GG\T{31}){}$;\5
969 
\&{break};\6
970 
\4${}\}{}$\2\6
971 
\&{if} ${}(\|s\E\.{'-'}){}$\1\5
972 
${}\|z\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
973 
\&{return} \|z;\6
974 
\4${}\}{}$\2\par
975 
\fi
976 
 
977 
\N{1}{41}Floating multiplication and division.
978 
The hardest fixed point operations were multiplication and division;
979 
but these two operations are the {\it easiest\/} to implement in floating point
980 
arithmetic, once their fixed point counterparts are available.
981 
 
982 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
983 
\&{octa} \\{fmult}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
984 
\hbox{}\6{}\&{octa} ${}\\{fmult}(\|y,\39\|z){}$\1\1\6
985 
\&{octa} \|y${},{}$ \|z;\2\2\6
986 
${}\{{}$\1\6
987 
\&{ftype} \\{yt}${},{}$ \\{zt};\6
988 
\&{int} \\{ye}${},{}$ \\{ze};\6
989 
\&{char} \\{ys}${},{}$ \\{zs};\6
990 
\&{octa} \|x${},{}$ \\{xf}${},{}$ \\{yf}${},{}$ \\{zf};\6
991 
\&{register} \&{int} \\{xe};\6
992 
\&{register} \&{char} \\{xs};\7
993 
${}\\{yt}\K\\{funpack}(\|y,\39{\AND}\\{yf},\39{\AND}\\{ye},\39{\AND}\\{ys});{}$%
994 
\6
995 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
996 
\6
997 
${}\\{xs}\K\\{ys}+\\{zs}-\.{'+'}{}$;\C{ will be \PB{\.{'-'}} when the result is
998 
negative }\6
999 
\&{switch} ${}(\T{4}*\\{yt}+\\{zt}){}$\5
1000 
${}\{{}$\1\6
1001 
\hbox{\4}\X42:The usual NaN cases\X;\6
1002 
\4\&{case} \T{4}${}*\\{zro}+\\{zro}{}$:\5
1003 
\&{case} \T{4}${}*\\{zro}+\\{num}{}$:\5
1004 
\&{case} \T{4}${}*\\{num}+\\{zro}{}$:\5
1005 
${}\|x\K\\{zero\_octa}{}$;\5
1006 
\&{break};\6
1007 
\4\&{case} \T{4}${}*\\{num}+\\{inf}{}$:\5
1008 
\&{case} \T{4}${}*\\{inf}+\\{num}{}$:\5
1009 
\&{case} \T{4}${}*\\{inf}+\\{inf}{}$:\5
1010 
${}\|x\K\\{inf\_octa}{}$;\5
1011 
\&{break};\6
1012 
\4\&{case} \T{4}${}*\\{zro}+\\{inf}{}$:\5
1013 
\&{case} \T{4}${}*\\{inf}+\\{zro}{}$:\5
1014 
${}\|x\K\\{standard\_NaN};{}$\6
1015 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\5
1016 
\&{break};\6
1017 
\4\&{case} \T{4}${}*\\{num}+\\{num}{}$:\5
1018 
\X43:Multiply nonzero numbers and \PB{\&{return}}\X;\6
1019 
\4${}\}{}$\2\6
1020 
\&{if} ${}(\\{xs}\E\.{'-'}){}$\1\5
1021 
${}\|x.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
1022 
\&{return} \|x;\6
1023 
\4${}\}{}$\2\par
1024 
\fi
1025 
 
1026 
\M{42}\B\X42:The usual NaN cases\X${}\E{}$\6
1027 
\4\&{case} \T{4}${}*\\{nan}+\\{nan}{}$:\5
1028 
\&{if} ${}(\R(\|y.\|h\AND\T{\^80000})){}$\1\5
1029 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\C{ \PB{\|y} is signaling }\2\6
1030 
\4\&{case} \T{4}${}*\\{zro}+\\{nan}{}$:\5
1031 
\&{case} \T{4}${}*\\{num}+\\{nan}{}$:\5
1032 
\&{case} \T{4}${}*\\{inf}+\\{nan}{}$:\6
1033 
\&{if} ${}(\R(\|z.\|h\AND\T{\^80000})){}$\1\5
1034 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT},\39\|z.\|h\MRL{{\OR}{\K}}\T{%
1035 
\^80000};{}$\2\6
1036 
\&{return} \|z;\6
1037 
\4\&{case} \T{4}${}*\\{nan}+\\{zro}{}$:\5
1038 
\&{case} \T{4}${}*\\{nan}+\\{num}{}$:\5
1039 
\&{case} \T{4}${}*\\{nan}+\\{inf}{}$:\6
1040 
\&{if} ${}(\R(\|y.\|h\AND\T{\^80000})){}$\1\5
1041 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT},\39\|y.\|h\MRL{{\OR}{\K}}\T{%
1042 
\^80000};{}$\2\6
1043 
\&{return} \|y;\par
1044 
\Us41, 44, 46\ETs93.\fi
1045 
 
1046 
\M{43}\B\X43:Multiply nonzero numbers and \PB{\&{return}}\X${}\E{}$\6
1047 
$\\{xe}\K\\{ye}+\\{ze}-\T{\^3fd}{}$;\C{ the raw exponent }\6
1048 
${}\|x\K\\{omult}(\\{yf},\39\\{shift\_left}(\\{zf},\39\T{9}));{}$\6
1049 
\&{if} ${}(\\{aux}.\|h\G\T{\^400000}){}$\1\5
1050 
${}\\{xf}\K\\{aux};{}$\2\6
1051 
\&{else}\1\5
1052 
${}\\{xf}\K\\{shift\_left}(\\{aux},\39\T{1}),\39\\{xe}\MM;{}$\2\6
1053 
\&{if} ${}(\|x.\|h\V\|x.\|l){}$\1\5
1054 
${}\\{xf}.\|l\MRL{{\OR}{\K}}\T{1}{}$;\C{ adjust the sticky bit }\2\6
1055 
\&{return} \\{fpack}${}(\\{xf},\39\\{xe},\39\\{xs},\39\\{cur\_round}){}$;\par
1056 
\U41.\fi
1057 
 
1058 
\M{44}\B\X5:Subroutines\X${}\mathrel+\E{}$\6
1059 
\&{octa} \\{fdivide}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
1060 
\hbox{}\6{}\&{octa} ${}\\{fdivide}(\|y,\39\|z){}$\1\1\6
1061 
\&{octa} \|y${},{}$ \|z;\2\2\6
1062 
${}\{{}$\1\6
1063 
\&{ftype} \\{yt}${},{}$ \\{zt};\6
1064 
\&{int} \\{ye}${},{}$ \\{ze};\6
1065 
\&{char} \\{ys}${},{}$ \\{zs};\6
1066 
\&{octa} \|x${},{}$ \\{xf}${},{}$ \\{yf}${},{}$ \\{zf};\6
1067 
\&{register} \&{int} \\{xe};\6
1068 
\&{register} \&{char} \\{xs};\7
1069 
${}\\{yt}\K\\{funpack}(\|y,\39{\AND}\\{yf},\39{\AND}\\{ye},\39{\AND}\\{ys});{}$%
1070 
\6
1071 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
1072 
\6
1073 
${}\\{xs}\K\\{ys}+\\{zs}-\.{'+'}{}$;\C{ will be \PB{\.{'-'}} when the result is
1074 
negative }\6
1075 
\&{switch} ${}(\T{4}*\\{yt}+\\{zt}){}$\5
1076 
${}\{{}$\1\6
1077 
\hbox{\4}\X42:The usual NaN cases\X;\6
1078 
\4\&{case} \T{4}${}*\\{zro}+\\{inf}{}$:\5
1079 
\&{case} \T{4}${}*\\{zro}+\\{num}{}$:\5
1080 
\&{case} \T{4}${}*\\{num}+\\{inf}{}$:\5
1081 
${}\|x\K\\{zero\_octa}{}$;\5
1082 
\&{break};\6
1083 
\4\&{case} \T{4}${}*\\{num}+\\{zro}{}$:\5
1084 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{Z\_BIT};{}$\6
1085 
\4\&{case} \T{4}${}*\\{inf}+\\{num}{}$:\5
1086 
\&{case} \T{4}${}*\\{inf}+\\{zro}{}$:\5
1087 
${}\|x\K\\{inf\_octa}{}$;\5
1088 
\&{break};\6
1089 
\4\&{case} \T{4}${}*\\{zro}+\\{zro}{}$:\5
1090 
\&{case} \T{4}${}*\\{inf}+\\{inf}{}$:\5
1091 
${}\|x\K\\{standard\_NaN};{}$\6
1092 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\5
1093 
\&{break};\6
1094 
\4\&{case} \T{4}${}*\\{num}+\\{num}{}$:\5
1095 
\X45:Divide nonzero numbers and \PB{\&{return}}\X;\6
1096 
\4${}\}{}$\2\6
1097 
\&{if} ${}(\\{xs}\E\.{'-'}){}$\1\5
1098 
${}\|x.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
1099 
\&{return} \|x;\6
1100 
\4${}\}{}$\2\par
1101 
\fi
1102 
 
1103 
\M{45}\B\X45:Divide nonzero numbers and \PB{\&{return}}\X${}\E{}$\6
1104 
$\\{xe}\K\\{ye}-\\{ze}+\T{\^3fd}{}$;\C{ the raw exponent }\6
1105 
${}\\{xf}\K\\{odiv}(\\{yf},\39\\{zero\_octa},\39\\{shift\_left}(\\{zf},\39%
1106 
\T{9}));{}$\6
1107 
\&{if} ${}(\\{xf}.\|h\G\T{\^800000}){}$\5
1108 
${}\{{}$\1\6
1109 
${}\\{aux}.\|l\MRL{{\OR}{\K}}\\{xf}.\|l\AND\T{1};{}$\6
1110 
${}\\{xf}\K\\{shift\_right}(\\{xf},\39\T{1},\39\T{1});{}$\6
1111 
${}\\{xe}\PP;{}$\6
1112 
\4${}\}{}$\2\6
1113 
\&{if} ${}(\\{aux}.\|h\V\\{aux}.\|l){}$\1\5
1114 
${}\\{xf}.\|l\MRL{{\OR}{\K}}\T{1}{}$;\C{ adjust the sticky bit }\2\6
1115 
\&{return} \\{fpack}${}(\\{xf},\39\\{xe},\39\\{xs},\39\\{cur\_round}){}$;\par
1116 
\U44.\fi
1117 
 
1118 
\N{1}{46}Floating addition and subtraction. Now for the bread-and-butter
1119 
operation, the sum of two floating point numbers.
1120 
It is not terribly difficult, but many cases need to be handled carefully.
1121 
 
1122 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
1123 
\&{octa} \\{fplus}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
1124 
\hbox{}\6{}\&{octa} ${}\\{fplus}(\|y,\39\|z){}$\1\1\6
1125 
\&{octa} \|y${},{}$ \|z;\2\2\6
1126 
${}\{{}$\1\6
1127 
\&{ftype} \\{yt}${},{}$ \\{zt};\6
1128 
\&{int} \\{ye}${},{}$ \\{ze};\6
1129 
\&{char} \\{ys}${},{}$ \\{zs};\6
1130 
\&{octa} \|x${},{}$ \\{xf}${},{}$ \\{yf}${},{}$ \\{zf};\6
1131 
\&{register} \&{int} \\{xe}${},{}$ \|d;\6
1132 
\&{register} \&{char} \\{xs};\7
1133 
${}\\{yt}\K\\{funpack}(\|y,\39{\AND}\\{yf},\39{\AND}\\{ye},\39{\AND}\\{ys});{}$%
1134 
\6
1135 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
1136 
\6
1137 
\&{switch} ${}(\T{4}*\\{yt}+\\{zt}){}$\5
1138 
${}\{{}$\1\6
1139 
\hbox{\4}\X42:The usual NaN cases\X;\6
1140 
\4\&{case} \T{4}${}*\\{zro}+\\{num}{}$:\5
1141 
\&{return} \\{fpack}${}(\\{zf},\39\\{ze},\39\\{zs},\39\.{ROUND\_OFF}){}$;\5
1142 
\&{break};\C{ may underflow }\6
1143 
\4\&{case} \T{4}${}*\\{num}+\\{zro}{}$:\5
1144 
\&{return} \\{fpack}${}(\\{yf},\39\\{ye},\39\\{ys},\39\.{ROUND\_OFF}){}$;\5
1145 
\&{break};\C{ may underflow }\6
1146 
\4\&{case} \T{4}${}*\\{inf}+\\{inf}{}$:\5
1147 
\&{if} ${}(\\{ys}\I\\{zs}){}$\5
1148 
${}\{{}$\1\6
1149 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\5
1150 
${}\|x\K\\{standard\_NaN}{}$;\5
1151 
${}\\{xs}\K\\{zs}{}$;\5
1152 
\&{break};\6
1153 
\4${}\}{}$\2\6
1154 
\4\&{case} \T{4}${}*\\{num}+\\{inf}{}$:\5
1155 
\&{case} \T{4}${}*\\{zro}+\\{inf}{}$:\5
1156 
${}\|x\K\\{inf\_octa}{}$;\5
1157 
${}\\{xs}\K\\{zs}{}$;\5
1158 
\&{break};\6
1159 
\4\&{case} \T{4}${}*\\{inf}+\\{num}{}$:\5
1160 
\&{case} \T{4}${}*\\{inf}+\\{zro}{}$:\5
1161 
${}\|x\K\\{inf\_octa}{}$;\5
1162 
${}\\{xs}\K\\{ys}{}$;\5
1163 
\&{break};\6
1164 
\4\&{case} \T{4}${}*\\{num}+\\{num}{}$:\5
1165 
\&{if} ${}(\|y.\|h\I(\|z.\|h\XOR\T{\^80000000})\V\|y.\|l\I\|z.\|l){}$\1\5
1166 
\X47:Add nonzero numbers and \PB{\&{return}}\X;\2\6
1167 
\4\&{case} \T{4}${}*\\{zro}+\\{zro}{}$:\5
1168 
${}\|x\K\\{zero\_octa};{}$\6
1169 
${}\\{xs}\K(\\{ys}\E\\{zs}\?\\{ys}:\\{cur\_round}\E\.{ROUND\_DOWN}\?\.{'-'}:%
1170 
\.{'+'}){}$;\5
1171 
\&{break};\6
1172 
\4${}\}{}$\2\6
1173 
\&{if} ${}(\\{xs}\E\.{'-'}){}$\1\5
1174 
${}\|x.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
1175 
\&{return} \|x;\6
1176 
\4${}\}{}$\2\par
1177 
\fi
1178 
 
1179 
\M{47}\B\X47:Add nonzero numbers and \PB{\&{return}}\X${}\E{}$\6
1180 
${}\{{}$\5
1181 
\1\&{octa} \|o${},{}$ \\{oo};\7
1182 
\&{if} ${}(\\{ye}<\\{ze}\V(\\{ye}\E\\{ze}\W(\\{yf}.\|h<\\{zf}.\|h\V(\\{yf}.\|h%
1183 
\E\\{zf}.\|h\W\\{yf}.\|l<\\{zf}.\|l)))){}$\1\5
1184 
\X48:Exchange \PB{\|y} with \PB{\|z}\X;\2\6
1185 
${}\|d\K\\{ye}-\\{ze};{}$\6
1186 
${}\\{xs}\K\\{ys},\39\\{xe}\K\\{ye};{}$\6
1187 
\&{if} (\|d)\1\5
1188 
\X49:Adjust for difference in exponents\X;\2\6
1189 
\&{if} ${}(\\{ys}\E\\{zs}){}$\5
1190 
${}\{{}$\1\6
1191 
${}\\{xf}\K\\{oplus}(\\{yf},\39\\{zf});{}$\6
1192 
\&{if} ${}(\\{xf}.\|h\G\T{\^800000}){}$\1\5
1193 
${}\\{xe}\PP,\39\|d\K\\{xf}.\|l\AND\T{1},\39\\{xf}\K\\{shift\_right}(\\{xf},\39%
1194 
\T{1},\39\T{1}),\39\\{xf}.\|l\MRL{{\OR}{\K}}\|d;{}$\2\6
1195 
\4${}\}{}$\5
1196 
\2\&{else}\5
1197 
${}\{{}$\1\6
1198 
${}\\{xf}\K\\{ominus}(\\{yf},\39\\{zf});{}$\6
1199 
\&{if} ${}(\\{xf}.\|h\G\T{\^800000}){}$\1\5
1200 
${}\\{xe}\PP,\39\|d\K\\{xf}.\|l\AND\T{1},\39\\{xf}\K\\{shift\_right}(\\{xf},\39%
1201 
\T{1},\39\T{1}),\39\\{xf}.\|l\MRL{{\OR}{\K}}\|d;{}$\2\6
1202 
\&{else}\5
1203 
\1\&{while} ${}(\\{xf}.\|h<\T{\^400000}){}$\1\5
1204 
${}\\{xe}\MM,\39\\{xf}\K\\{shift\_left}(\\{xf},\39\T{1});{}$\2\2\6
1205 
\4${}\}{}$\2\6
1206 
\&{return} \\{fpack}${}(\\{xf},\39\\{xe},\39\\{xs},\39\\{cur\_round});{}$\6
1207 
\4${}\}{}$\2\par
1208 
\U46.\fi
1209 
 
1210 
\M{48}\B\X48:Exchange \PB{\|y} with \PB{\|z}\X${}\E{}$\6
1211 
${}\{{}$\1\6
1212 
${}\|o\K\\{yf},\39\\{yf}\K\\{zf},\39\\{zf}\K\|o;{}$\6
1213 
${}\|d\K\\{ye},\39\\{ye}\K\\{ze},\39\\{ze}\K\|d;{}$\6
1214 
${}\|d\K\\{ys},\39\\{ys}\K\\{zs},\39\\{zs}\K\|d;{}$\6
1215 
\4${}\}{}$\2\par
1216 
\Us47\ET51.\fi
1217 
 
1218 
\M{49}Proper rounding requires two bits to the right of the fraction delivered
1219 
to~\PB{\\{fpack}}. The first is the true next bit of the result;
1220 
the other is a ``sticky'' bit, which is nonzero if any further bits of the
1221 
true result are nonzero. Sticky rounding to an integer takes
1222 
$x$ into the number $\lfloor x/2\rfloor+\lceil x/2\rceil$.
1223 
 
1224 
Some subtleties need to be observed here, in order to
1225 
prevent the sticky bit from being shifted left. If we did not
1226 
shift \PB{\\{yf}} left~1 before shifting \PB{\\{zf}} to the right, an incorrect
1227 
answer would be obtained in certain cases---for example, if
1228 
$\PB{\\{yf}}=2^{54}$, $\PB{\\{zf}}=2^{54}+2^{53}-1$, $d=52$.
1229 
 
1230 
\Y\B\4\X49:Adjust for difference in exponents\X${}\E{}$\6
1231 
${}\{{}$\1\6
1232 
\&{if} ${}(\|d\Z\T{2}){}$\1\5
1233 
${}\\{zf}\K\\{shift\_right}(\\{zf},\39\|d,\39\T{1}){}$;\C{ exact result }\2\6
1234 
\&{else} \&{if} ${}(\|d>\T{53}){}$\1\5
1235 
${}\\{zf}.\|h\K\T{0},\39\\{zf}.\|l\K\T{1}{}$;\C{ tricky but OK }\2\6
1236 
\&{else}\5
1237 
${}\{{}$\1\6
1238 
\&{if} ${}(\\{ys}\I\\{zs}){}$\1\5
1239 
${}\|d\MM,\39\\{xe}\MM,\39\\{yf}\K\\{shift\_left}(\\{yf},\39\T{1});{}$\2\6
1240 
${}\|o\K\\{zf};{}$\6
1241 
${}\\{zf}\K\\{shift\_right}(\|o,\39\|d,\39\T{1});{}$\6
1242 
${}\\{oo}\K\\{shift\_left}(\\{zf},\39\|d);{}$\6
1243 
\&{if} ${}(\\{oo}.\|l\I\|o.\|l\V\\{oo}.\|h\I\|o.\|h){}$\1\5
1244 
${}\\{zf}.\|l\MRL{{\OR}{\K}}\T{1};{}$\2\6
1245 
\4${}\}{}$\2\6
1246 
\4${}\}{}$\2\par
1247 
\U47.\fi
1248 
 
1249 
\M{50}The comparison of floating point numbers with respect to $\epsilon$
1250 
shares some of the characteristics of floating point addition/subtraction.
1251 
In some ways it is simpler, and in other ways it is more difficult;
1252 
we might as well deal with it now. % anyways
1253 
 
1254 
Subroutine \PB{$\\{fepscomp}(\|y,\|z,\|e,\|s)$} returns 2 if \PB{\|y}, \PB{%
1255 
\|z}, or \PB{\|e} is a NaN
1256 
or \PB{\|e} is negative. It returns 1 if \PB{$\|s\K\T{0}$} and $y\approx z\
1257 
(e)$ or if
1258 
\PB{$\|s\I\T{0}$} and $y\sim z\ (e)$,
1259 
as defined in Section~4.2.2 of {\sl Seminumerical Algorithms\/};
1260 
otherwise it returns~0.
1261 
 
1262 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
1263 
\&{int} \\{fepscomp}\,\,${}\.{ARGS}((\&{octa},\39\&{octa},\39\&{octa},\39%
1264 
\&{int})){}$;\5
1265 
\hbox{}\6{}\&{int} ${}\\{fepscomp}(\|y,\39\|z,\39\|e,\39\|s){}$\1\1\6
1266 
\&{octa} \|y${},{}$ \|z${},{}$ \|e;\C{ the operands }\6
1267 
\&{int} \|s;\C{ test similarity? }\2\2\6
1268 
${}\{{}$\1\6
1269 
\&{octa} \\{yf}${},{}$ \\{zf}${},{}$ \\{ef}${},{}$ \|o${},{}$ \\{oo};\6
1270 
\&{int} \\{ye}${},{}$ \\{ze}${},{}$ \\{ee};\6
1271 
\&{char} \\{ys}${},{}$ \\{zs}${},{}$ \\{es};\6
1272 
\&{register} \&{int} \\{yt}${},{}$ \\{zt}${},{}$ \\{et}${},{}$ \|d;\7
1273 
${}\\{et}\K\\{funpack}(\|e,\39{\AND}\\{ef},\39{\AND}\\{ee},\39{\AND}\\{es});{}$%
1274 
\6
1275 
\&{if} ${}(\\{es}\E\.{'-'}){}$\1\5
1276 
\&{return} \T{2};\2\6
1277 
\&{switch} (\\{et})\5
1278 
${}\{{}$\1\6
1279 
\4\&{case} \\{nan}:\5
1280 
\&{return} \T{2};\6
1281 
\4\&{case} \\{inf}:\5
1282 
${}\\{ee}\K\T{10000};{}$\6
1283 
\4\&{case} \\{num}:\5
1284 
\&{case} \\{zro}:\5
1285 
\&{break};\6
1286 
\4${}\}{}$\2\6
1287 
${}\\{yt}\K\\{funpack}(\|y,\39{\AND}\\{yf},\39{\AND}\\{ye},\39{\AND}\\{ys});{}$%
1288 
\6
1289 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
1290 
\6
1291 
\&{switch} ${}(\T{4}*\\{yt}+\\{zt}){}$\5
1292 
${}\{{}$\1\6
1293 
\4\&{case} \T{4}${}*\\{nan}+\\{nan}{}$:\5
1294 
\&{case} \T{4}${}*\\{nan}+\\{inf}{}$:\5
1295 
\&{case} \T{4}${}*\\{nan}+\\{num}{}$:\5
1296 
\&{case} \T{4}${}*\\{nan}+\\{zro}{}$:\5
1297 
\&{case} \T{4}${}*\\{inf}+\\{nan}{}$:\5
1298 
\&{case} \T{4}${}*\\{num}+\\{nan}{}$:\5
1299 
\&{case} \T{4}${}*\\{zro}+\\{nan}{}$:\5
1300 
\&{return} \T{2};\6
1301 
\4\&{case} \T{4}${}*\\{inf}+\\{inf}{}$:\5
1302 
\&{return} ${}(\\{ys}\E\\{zs}\V\\{ee}\G\T{1023});{}$\6
1303 
\4\&{case} \T{4}${}*\\{inf}+\\{num}{}$:\5
1304 
\&{case} \T{4}${}*\\{inf}+\\{zro}{}$:\5
1305 
\&{case} \T{4}${}*\\{num}+\\{inf}{}$:\5
1306 
\&{case} \T{4}${}*\\{zro}+\\{inf}{}$:\5
1307 
\&{return} ${}(\|s\W\\{ee}\G\T{1022});{}$\6
1308 
\4\&{case} \T{4}${}*\\{zro}+\\{zro}{}$:\5
1309 
\&{return} \T{1};\6
1310 
\4\&{case} \T{4}${}*\\{zro}+\\{num}{}$:\5
1311 
\&{case} \T{4}${}*\\{num}+\\{zro}{}$:\5
1312 
\&{if} ${}(\R\|s){}$\1\5
1313 
\&{return} \T{0};\2\6
1314 
\4\&{case} \T{4}${}*\\{num}+\\{num}{}$:\5
1315 
\&{break};\6
1316 
\4${}\}{}$\2\6
1317 
\X51:Compare two numbers with respect to epsilon and \PB{\&{return}}\X;\6
1318 
\4${}\}{}$\2\par
1319 
\fi
1320 
 
1321 
\M{51}The relation $y\approx z\ (\epsilon)$ reduces to
1322 
$y\sim z\ (\epsilon/2^d)$, if $d$~is the difference between the
1323 
larger and smaller exponents of $y$ and~$z$.
1324 
 
1325 
\Y\B\4\X51:Compare two numbers with respect to epsilon and \PB{\&{return}}\X${}%
1326 
\E{}$\6
1327 
\X52:Unsubnormalize \PB{\|y} and \PB{\|z}, if they are subnormal\X;\6
1328 
\&{if} ${}(\\{ye}<\\{ze}\V(\\{ye}\E\\{ze}\W(\\{yf}.\|h<\\{zf}.\|h\V(\\{yf}.\|h%
1329 
\E\\{zf}.\|h\W\\{yf}.\|l<\\{zf}.\|l)))){}$\1\5
1330 
\X48:Exchange \PB{\|y} with \PB{\|z}\X;\2\6
1331 
\&{if} ${}(\\{ze}\E\\{zero\_exponent}){}$\1\5
1332 
${}\\{ze}\K\\{ye};{}$\2\6
1333 
${}\|d\K\\{ye}-\\{ze};{}$\6
1334 
\&{if} ${}(\R\|s){}$\1\5
1335 
${}\\{ee}\MRL{-{\K}}\|d;{}$\2\6
1336 
\&{if} ${}(\\{ee}\G\T{1023}){}$\1\5
1337 
\&{return} \T{1};\C{ if $\epsilon\ge2$, $z\in N_\epsilon(y)$ }\2\6
1338 
\X53:Compute the difference of fraction parts, \PB{\|o}\X;\6
1339 
\&{if} ${}(\R\|o.\|h\W\R\|o.\|l){}$\1\5
1340 
\&{return} \T{1};\2\6
1341 
\&{if} ${}(\\{ee}<\T{968}){}$\1\5
1342 
\&{return} \T{0};\C{ if $y\ne z$ and $\epsilon<2^{-54}$, $y\not\sim z$ }\2\6
1343 
\&{if} ${}(\\{ee}\G\T{1021}){}$\1\5
1344 
${}\\{ef}\K\\{shift\_left}(\\{ef},\39\\{ee}-\T{1021});{}$\2\6
1345 
\&{else}\1\5
1346 
${}\\{ef}\K\\{shift\_right}(\\{ef},\39\T{1021}-\\{ee},\39\T{1});{}$\2\6
1347 
\&{return} \|o${}.\|h<\\{ef}.\|h\V(\|o.\|h\E\\{ef}.\|h\W\|o.\|l\Z\\{ef}.%
1348 
\|l){}$;\par
1349 
\U50.\fi
1350 
 
1351 
\M{52}\B\X52:Unsubnormalize \PB{\|y} and \PB{\|z}, if they are subnormal\X${}%
1352 
\E{}$\6
1353 
\&{if} ${}(\\{ye}<\T{0}\W\\{yt}\I\\{zro}){}$\1\5
1354 
${}\\{yf}\K\\{shift\_left}(\|y,\39\T{2}),\39\\{ye}\K\T{0};{}$\2\6
1355 
\&{if} ${}(\\{ze}<\T{0}\W\\{zt}\I\\{zro}){}$\1\5
1356 
${}\\{zf}\K\\{shift\_left}(\|z,\39\T{2}),\39\\{ze}\K\T{0}{}$;\2\par
1357 
\U51.\fi
1358 
 
1359 
\M{53}At this point $y\sim z$ if and only if
1360 
$$\PB{\\{yf}}+(-1)^{[ys=zs]}\PB{\\{zf}}/2^d\le 2^{ee-1021}\PB{\\{ef}}=2^{55}%
1361 
\epsilon.$$
1362 
We need to evaluate this relation without overstepping the bounds of
1363 
our simulated 64-bit registers.
1364 
 
1365 
When $d>2$, the difference of fraction parts might not fit exactly
1366 
in an octabyte;
1367 
in that case the numbers are not similar unless $\epsilon>3/8$,
1368 
and we replace the difference by the ceiling of the
1369 
true result. When $\epsilon<1/8$, our program essentially replaces
1370 
$2^{55}\epsilon$ by $\lfloor2^{55}\epsilon\rfloor$. These
1371 
truncations are not needed simultaneously. Therefore the logic
1372 
is justified by the facts that, if $n$ is an integer, we have
1373 
$x\le n$ if and only if $\lceil x\rceil\le n$;
1374 
$n\le x$ if and only if $n\le\lfloor x\rfloor$. (Notice that the
1375 
concept of ``sticky bit'' is {\it not\/} appropriate here.)
1376 
 
1377 
\Y\B\4\X53:Compute the difference of fraction parts, \PB{\|o}\X${}\E{}$\6
1378 
\&{if} ${}(\|d>\T{54}){}$\1\5
1379 
${}\|o\K\\{zero\_octa},\39\\{oo}\K\\{zf};{}$\2\6
1380 
\&{else}\1\5
1381 
${}\|o\K\\{shift\_right}(\\{zf},\39\|d,\39\T{1}),\39\\{oo}\K\\{shift\_left}(%
1382 
\|o,\39\|d);{}$\2\6
1383 
\&{if} ${}(\\{oo}.\|h\I\\{zf}.\|h\V\\{oo}.\|l\I\\{zf}.\|l){}$\5
1384 
${}\{{}$\C{ truncated result, hence $d>2$ }\1\6
1385 
\&{if} ${}(\\{ee}<\T{1020}){}$\1\5
1386 
\&{return} \T{0};\C{ difference is too large for similarity }\2\6
1387 
${}\|o\K\\{incr}(\|o,\39\\{ys}\E\\{zs}\?\T{0}:\T{1}){}$;\C{ adjust for ceiling
1388 
}\6
1389 
\4${}\}{}$\2\6
1390 
${}\|o\K(\\{ys}\E\\{zs}\?\\{ominus}(\\{yf},\39\|o):\\{oplus}(\\{yf},\39%
1391 
\|o)){}$;\par
1392 
\U51.\fi
1393 
 
1394 
\N{1}{54}Floating point output conversion.
1395 
The \PB{\\{print\_float}} routine converts an octabyte to a floating decimal
1396 
representation that will be input as precisely the same value.
1397 
 
1398 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
1399 
\&{static} \&{void} \\{bignum\_times\_ten}\,\,${}\.{ARGS}((\\{bignum}*));{}$\6
1400 
\&{static} \&{void} \\{bignum\_dec}\,\,${}\.{ARGS}((\\{bignum}*,\\{bignum}*,%
1401 
\&{tetra}));{}$\6
1402 
\&{static} \&{int} \\{bignum\_compare}\,\,${}\.{ARGS}((\\{bignum}*,%
1403 
\\{bignum}*));{}$\6
1404 
\&{void} \\{print\_float}\,\,\.{ARGS}((\&{octa}));\5
1405 
\hbox{}\6{}\&{void} \\{print\_float}(\|x)\1\1\6
1406 
\&{octa} \|x;\2\2\6
1407 
${}\{{}$\1\6
1408 
\X56:Local variables for \PB{\\{print\_float}}\X;\6
1409 
\&{if} ${}(\|x.\|h\AND\\{sign\_bit}){}$\1\5
1410 
\\{printf}(\.{"-"});\2\6
1411 
\X55:Extract the exponent \PB{\|e} and determine the fraction interval $[f\dts
1412 
g]$ or $(f\dts g)$\X;\6
1413 
\X63:Store $f$ and $g$ as multiprecise integers\X;\6
1414 
\X64:Compute the significant digits \PB{\|s} and decimal exponent \PB{\|e}\X;\6
1415 
\X67:Print the significant digits with proper context\X;\6
1416 
\4${}\}{}$\2\par
1417 
\fi
1418 
 
1419 
\M{55}One way to visualize the problem being solved here is to consider
1420 
the vastly simpler case in which there are only 2-bit exponents
1421 
and 2-bit fractions. Then the sixteen possible 4-bit combinations
1422 
have the following interpretations:
1423 
$$\def\\{\;\dts\;}
1424 
\vbox{\halign{#\qquad&$#$\hfil\cr
1425 
0000&[0\\0.125]\cr
1426 
0001&(0.125\\0.375)\cr
1427 
0010&[0.375\\0.625]\cr
1428 
0011&(0.625\\0.875)\cr
1429 
0100&[0.875\\1.125]\cr
1430 
0101&(1.125\\1.375)\cr
1431 
0110&[1.375\\1.625]\cr
1432 
0111&(1.625\\1.875)\cr
1433 
1000&[1.875\\2.25]\cr
1434 
1001&(2.25\\2.75)\cr
1435 
1010&[2.75\\3.25]\cr
1436 
1011&(3.25\\3.75)\cr
1437 
1100&[3.75\\\infty]\cr
1438 
1101&\rm NaN(0\\0.375)\cr
1439 
1110&\rm NaN[0.375\\0.625]\cr
1440 
1111&\rm NaN(0.625\\1)\cr}}$$
1441 
Notice that the interval is closed, $[f\dts g]$, when the fraction part
1442 
is even; it is open, $(f\dts g)$, when the fraction part is odd.
1443 
The printed outputs for these sixteen values, if we actually were
1444 
dealing with such short exponents and fractions, would be
1445 
\.{0.}, \.{.2}, \.{.5}, \.{.7}, \.{1.}, \.{1.2}, \.{1.5}, \.{1.7},
1446 
\.{2.}, \.{2.5}, \.{3.}, \.{3.5}, \.{Inf}, \.{NaN.2}, \.{NaN}, \.{NaN.8},
1447 
respectively.
1448 
 
1449 
\Y\B\4\X55:Extract the exponent \PB{\|e} and determine the fraction interval
1450 
$[f\dts g]$ or $(f\dts g)$\X${}\E{}$\6
1451 
$\|f\K\\{shift\_left}(\|x,\39\T{1});{}$\6
1452 
${}\|e\K\|f.\|h\GG\T{21};{}$\6
1453 
${}\|f.\|h\MRL{\AND{\K}}\T{\^1fffff};{}$\6
1454 
\&{if} ${}(\R\|f.\|h\W\R\|f.\|l){}$\1\5
1455 
\X57:Handle the special case when the fraction part is zero\X\2\6
1456 
\&{else}\5
1457 
${}\{{}$\1\6
1458 
${}\|g\K\\{incr}(\|f,\39\T{1});{}$\6
1459 
${}\|f\K\\{incr}(\|f,\39{-}\T{1});{}$\6
1460 
\&{if} ${}(\R\|e){}$\1\5
1461 
${}\|e\K\T{1}{}$;\C{ subnormal }\2\6
1462 
\&{else} \&{if} ${}(\|e\E\T{\^7ff}){}$\5
1463 
${}\{{}$\1\6
1464 
\\{printf}(\.{"NaN"});\6
1465 
\&{if} ${}(\|g.\|h\E\T{\^100000}\W\|g.\|l\E\T{1}){}$\1\5
1466 
\&{return};\C{ the ``standard'' NaN }\2\6
1467 
${}\|e\K\T{\^3ff}{}$;\C{ extreme NaNs come out OK even without adjusting \PB{%
1468 
\|f} or \PB{\|g} }\6
1469 
\4${}\}{}$\5
1470 
\2\&{else}\1\5
1471 
${}\|f.\|h\MRL{{\OR}{\K}}\T{\^200000},\39\|g.\|h\MRL{{\OR}{\K}}\T{\^200000};{}$%
1472 
\2\6
1473 
\4${}\}{}$\2\par
1474 
\U54.\fi
1475 
 
1476 
\M{56}\B\X56:Local variables for \PB{\\{print\_float}}\X${}\E{}$\6
1477 
\&{octa} \|f${},{}$ \|g;\C{ lower and upper bounds on the fraction part }\6
1478 
\&{register} \&{int} \|e;\C{ exponent part }\6
1479 
\&{register} \&{int} \|j${},{}$ \|k;\C{ all purpose indices }\par
1480 
\A66.
1481 
\U54.\fi
1482 
 
1483 
\M{57}The transition points between exponents correspond to powers of~2. At
1484 
such points the interval extends only half as far to the left of that
1485 
power of~2 as it does to the right. For example, in the 4-bit minifloat numbers
1486 
considered above, case 1000 corresponds to the interval $[1.875\;\dts\;2.25]$.
1487 
 
1488 
\Y\B\4\X57:Handle the special case when the fraction part is zero\X${}\E{}$\6
1489 
${}\{{}$\1\6
1490 
\&{if} ${}(\R\|e){}$\5
1491 
${}\{{}$\1\6
1492 
\\{printf}(\.{"0."});\5
1493 
\&{return};\6
1494 
\4${}\}{}$\2\6
1495 
\&{if} ${}(\|e\E\T{\^7ff}){}$\5
1496 
${}\{{}$\1\6
1497 
\\{printf}(\.{"Inf"});\5
1498 
\&{return};\6
1499 
\4${}\}{}$\2\6
1500 
${}\|e\MM;{}$\6
1501 
${}\|f.\|h\K\T{\^3fffff},\39\|f.\|l\K\T{\^ffffffff};{}$\6
1502 
${}\|g.\|h\K\T{\^400000},\39\|g.\|l\K\T{2};{}$\6
1503 
\4${}\}{}$\2\par
1504 
\U55.\fi
1505 
 
1506 
\M{58}We want to find the ``simplest'' value in the interval corresponding
1507 
to the given number, in the sense that it has fewest significant
1508 
digits when expressed in decimal notation. Thus, for example,
1509 
if the floating point number can be described by a relatively
1510 
short string such as `\.{.1}' or `\.{37e100}', we want to discover that
1511 
representation.
1512 
 
1513 
The basic idea is to generate the decimal representations of the
1514 
two endpoints of the interval, outputting the leading digits where
1515 
both endpoints agree, then making a final decision at the first place where
1516 
they disagree.
1517 
 
1518 
The ``simplest'' value is not always unique. For example, in the
1519 
case of 4-bit minifloat numbers we could represent the bit pattern 0001 as
1520 
either \.{.2} or \.{.3}, and we could represent 1001 in five equally short
1521 
ways: \.{2.3} or \.{2.4} or \.{2.5} or \.{2.6} or \.{2.7}. The
1522 
algorithm below tries to choose the middle possibility in such cases.
1523 
 
1524 
[A solution to the analogous problem for fixed-point representations,
1525 
without the additional complication of round-to-even, was used by
1526 
the author in the program for \TeX; see {\sl Beauty is Our Business\/}
1527 
(Springer, 1990), 233--242.]
1528 
 
1529 
Suppose we are given two fractions $f$ and $g$, where $0\le f<g<1$, and
1530 
we want to compute the shortest decimal in the closed interval $[f\dts g]$.
1531 
If $f=0$, we are done. Otherwise let $10f=d+f'$ and $10g=e+g'$, where
1532 
$0\le f'<1$ and $0\le g'<1$. If $d<e$, we can terminate by outputting
1533 
any of the digits $d+1$, \dots,~$e$; otherwise we output the
1534 
common digit $d=e$, and repeat the process on the fractions $0\le f'<g'<1$.
1535 
A similar procedure works with respect to the open interval $(f\dts g)$.
1536 
 
1537 
\fi
1538 
 
1539 
\M{59}The program below carries out the stated algorithm by using
1540 
multiprecision
1541 
arithmetic on 77-place integers with 28 bits each. This choice
1542 
facilitates multiplication by~10, and allows us to deal with the whole range of
1543 
floating binary numbers using fixed point arithmetic. We keep track of
1544 
the leading and trailing digit positions so that trivial operations on
1545 
zeros are avoided.
1546 
 
1547 
If \PB{\|f} points to a \&{bignum}, its radix-$2^{28}$ digits are
1548 
\PB{$\|f\MG\\{dat}[\T{0}]$} through \PB{$\|f\MG\\{dat}[\T{76}]$}, from most
1549 
significant to least significant.
1550 
We assume that all digit positions are zero unless they lie in the
1551 
subarray between indices \PB{$\|f\MG\|a$} and \PB{$\|f\MG\|b$}, inclusive.
1552 
Furthermore, both \PB{$\|f\MG\\{dat}[\|f\MG\|a]$} and \PB{$\|f\MG\\{dat}[\|f\MG%
1553 
\|b]$} are nonzero,
1554 
unless \PB{$\|f\MG\|a\K\|f\MG\|b\K\\{bignum\_prec}-\T{1}$}.
1555 
 
1556 
The \&{bignum} data type can be used with any radix less than
1557 
$2^{32}$; we will use it later with radix~$10^9$. The \PB{\\{dat}} array
1558 
is made large enough to accommodate both applications.
1559 
 
1560 
\Y\B\4\D$\\{bignum\_prec}$ \5
1561 
\T{157}\C{ would be 77 if we cared only about \PB{\\{print\_float}} }\par
1562 
\Y\B\4\X36:Other type definitions\X${}\mathrel+\E{}$\6
1563 
\&{typedef} \&{struct} ${}\{{}$\1\6
1564 
\&{int} \|a;\C{ index of the most significant digit }\6
1565 
\&{int} \|b;\C{ index of the least significant digit; must be $\ge a$ }\6
1566 
\&{tetra} \\{dat}[\\{bignum\_prec}];\C{ the digits; undefined except between %
1567 
\PB{\|a} and \PB{\|b} }\2\6
1568 
${}\}{}$ \&{bignum};\par
1569 
\fi
1570 
 
1571 
\M{60}Here, for example, is how we go from $f$ to $10f$, assuming that
1572 
overflow will not occur and that the radix is $2^{28}$:
1573 
 
1574 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
1575 
\&{static} \&{void} \\{bignum\_times\_ten}(\|f)\1\1\6
1576 
\&{bignum} ${}{*}\|f;\2\2{}$\6
1577 
${}\{{}$\1\6
1578 
\&{register} \&{tetra} ${}{*}\|p,{}$ ${}{*}\|q;{}$\6
1579 
\&{register} \&{tetra} \|x${},{}$ \\{carry};\7
1580 
\&{for} ${}(\|p\K{\AND}\|f\MG\\{dat}[\|f\MG\|b],\39\|q\K{\AND}\|f\MG\\{dat}[\|f%
1581 
\MG\|a],\39\\{carry}\K\T{0};{}$ ${}\|p\G\|q;{}$ ${}\|p\MM){}$\5
1582 
${}\{{}$\1\6
1583 
${}\|x\K{*}\|p*\T{10}+\\{carry};{}$\6
1584 
${}{*}\|p\K\|x\AND\T{\^fffffff};{}$\6
1585 
${}\\{carry}\K\|x\GG\T{28};{}$\6
1586 
\4${}\}{}$\2\6
1587 
${}{*}\|p\K\\{carry};{}$\6
1588 
\&{if} (\\{carry})\1\5
1589 
${}\|f\MG\|a\MM;{}$\2\6
1590 
\&{if} ${}(\|f\MG\\{dat}[\|f\MG\|b]\E\T{0}\W\|f\MG\|b>\|f\MG\|a){}$\1\5
1591 
${}\|f\MG\|b\MM;{}$\2\6
1592 
\4${}\}{}$\2\par
1593 
\fi
1594 
 
1595 
\M{61}And here is how we test whether $f<g$, $f=g$, or $f>g$, using any
1596 
radix whatever:
1597 
 
1598 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
1599 
\&{static} \&{int} ${}\\{bignum\_compare}(\|f,\39\|g){}$\1\1\6
1600 
\&{bignum} ${}{*}\|f,{}$ ${}{*}\|g;\2\2{}$\6
1601 
${}\{{}$\1\6
1602 
\&{register} \&{tetra} ${}{*}\|p,{}$ ${}{*}\\{pp},{}$ ${}{*}\|q,{}$ ${}{*}%
1603 
\\{qq};{}$\7
1604 
\&{if} ${}(\|f\MG\|a\I\|g\MG\|a){}$\1\5
1605 
\&{return} \|f${}\MG\|a>\|g\MG\|a\?{-}\T{1}:\T{1};{}$\2\6
1606 
${}\\{pp}\K{\AND}\|f\MG\\{dat}[\|f\MG\|b],\39\\{qq}\K{\AND}\|g\MG\\{dat}[\|g\MG%
1607 
\|b];{}$\6
1608 
\&{for} ${}(\|p\K{\AND}\|f\MG\\{dat}[\|f\MG\|a],\39\|q\K{\AND}\|g\MG\\{dat}[\|g%
1609 
\MG\|a];{}$ ${}\|p\Z\\{pp};{}$ ${}\|p\PP,\39\|q\PP){}$\5
1610 
${}\{{}$\1\6
1611 
\&{if} ${}({*}\|p\I{*}\|q){}$\1\5
1612 
\&{return} ${}{*}\|p<{*}\|q\?{-}\T{1}:\T{1};{}$\2\6
1613 
\&{if} ${}(\|q\E\\{qq}){}$\1\5
1614 
\&{return} \|p${}<\\{pp};{}$\2\6
1615 
\4${}\}{}$\2\6
1616 
\&{return} ${}{-}\T{1};{}$\6
1617 
\4${}\}{}$\2\par
1618 
\fi
1619 
 
1620 
\M{62}The following subroutine subtracts $g$ from~$f$, assuming that
1621 
$f\ge g>0$ and using a given radix.
1622 
 
1623 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
1624 
\&{static} \&{void} ${}\\{bignum\_dec}(\|f,\39\|g,\39\|r){}$\1\1\6
1625 
\&{bignum} ${}{*}\|f,{}$ ${}{*}\|g;{}$\6
1626 
\&{tetra} \|r;\C{ the radix }\2\2\6
1627 
${}\{{}$\1\6
1628 
\&{register} \&{tetra} ${}{*}\|p,{}$ ${}{*}\|q,{}$ ${}{*}\\{qq};{}$\6
1629 
\&{register} \&{int} \|x${},{}$ \\{borrow};\7
1630 
\&{while} ${}(\|g\MG\|b>\|f\MG\|b){}$\1\5
1631 
${}\|f\MG\\{dat}[\PP\|f\MG\|b]\K\T{0};{}$\2\6
1632 
${}\\{qq}\K{\AND}\|g\MG\\{dat}[\|g\MG\|a];{}$\6
1633 
\&{for} ${}(\|p\K{\AND}\|f\MG\\{dat}[\|g\MG\|b],\39\|q\K{\AND}\|g\MG\\{dat}[\|g%
1634 
\MG\|b],\39\\{borrow}\K\T{0};{}$ ${}\|q\G\\{qq};{}$ ${}\|p\MM,\39\|q\MM){}$\5
1635 
${}\{{}$\1\6
1636 
${}\|x\K{*}\|p-{*}\|q-\\{borrow};{}$\6
1637 
\&{if} ${}(\|x\G\T{0}){}$\1\5
1638 
${}\\{borrow}\K\T{0},\39{*}\|p\K\|x;{}$\2\6
1639 
\&{else}\1\5
1640 
${}\\{borrow}\K\T{1},\39{*}\|p\K\|x+\|r;{}$\2\6
1641 
\4${}\}{}$\2\6
1642 
\&{for} ( ; \\{borrow}; ${}\|p\MM){}$\1\6
1643 
\&{if} ${}({*}\|p){}$\1\5
1644 
${}\\{borrow}\K\T{0},\39{*}\|p\K{*}\|p-\T{1};{}$\2\6
1645 
\&{else}\1\5
1646 
${}{*}\|p\K\|r-\T{1};{}$\2\2\6
1647 
\&{while} ${}(\|f\MG\\{dat}[\|f\MG\|a]\E\T{0}){}$\5
1648 
${}\{{}$\1\6
1649 
\&{if} ${}(\|f\MG\|a\E\|f\MG\|b){}$\5
1650 
${}\{{}$\C{ the result is zero }\1\6
1651 
${}\|f\MG\|a\K\|f\MG\|b\K\\{bignum\_prec}-\T{1},\39\|f\MG\\{dat}[\\{bignum%
1652 
\_prec}-\T{1}]\K\T{0};{}$\6
1653 
\&{return};\6
1654 
\4${}\}{}$\2\6
1655 
${}\|f\MG\|a\PP;{}$\6
1656 
\4${}\}{}$\2\6
1657 
\&{while} ${}(\|f\MG\\{dat}[\|f\MG\|b]\E\T{0}){}$\1\5
1658 
${}\|f\MG\|b\MM;{}$\2\6
1659 
\4${}\}{}$\2\par
1660 
\fi
1661 
 
1662 
\M{63}Armed with these subroutines, we are ready to solve the problem.
1663 
The first task is to put the numbers into \&{bignum} form.
1664 
If the exponent is \PB{\|e}, the number destined for digit \PB{\\{dat}[\|k]}
1665 
will
1666 
consist of the rightmost 28 bits of the given fraction after it has
1667 
been shifted right $c-e-28k$ bits, for some constant~$c$.
1668 
We choose $c$ so that,
1669 
when $e$ has its maximum value \Hex{7ff}, the leading digit will
1670 
go into position \PB{\\{dat}[\T{1}]}, and so that when the number to be printed
1671 
is exactly~1 the integer part of~$g$ will also be exactly~1.
1672 
 
1673 
\Y\B\4\D$\\{magic\_offset}$ \5
1674 
\T{2112}\C{ the constant $c$ that makes it work }\par
1675 
\B\4\D$\\{origin}$ \5
1676 
\T{37}\C{ the radix point follows \PB{\\{dat}[\T{37}]} }\par
1677 
\Y\B\4\X63:Store $f$ and $g$ as multiprecise integers\X${}\E{}$\6
1678 
$\|k\K(\\{magic\_offset}-\|e)/\T{28};{}$\6
1679 
${}\ff.\\{dat}[\|k-\T{1}]\K\\{shift\_right}(\|f,\39\\{magic\_offset}+\T{28}-%
1680 
\|e-\T{28}*\|k,\39\T{1}).\|l\AND\T{\^fffffff};{}$\6
1681 
${}\\{gg}.\\{dat}[\|k-\T{1}]\K\\{shift\_right}(\|g,\39\\{magic\_offset}+\T{28}-%
1682 
\|e-\T{28}*\|k,\39\T{1}).\|l\AND\T{\^fffffff};{}$\6
1683 
${}\ff.\\{dat}[\|k]\K\\{shift\_right}(\|f,\39\\{magic\_offset}-\|e-\T{28}*\|k,%
1684 
\39\T{1}).\|l\AND\T{\^fffffff};{}$\6
1685 
${}\\{gg}.\\{dat}[\|k]\K\\{shift\_right}(\|g,\39\\{magic\_offset}-\|e-\T{28}*%
1686 
\|k,\39\T{1}).\|l\AND\T{\^fffffff};{}$\6
1687 
${}\ff.\\{dat}[\|k+\T{1}]\K\\{shift\_left}(\|f,\39\|e+\T{28}*\|k-(\\{magic%
1688 
\_offset}-\T{28})).\|l\AND\T{\^fffffff};{}$\6
1689 
${}\\{gg}.\\{dat}[\|k+\T{1}]\K\\{shift\_left}(\|g,\39\|e+\T{28}*\|k-(\\{magic%
1690 
\_offset}-\T{28})).\|l\AND\T{\^fffffff};{}$\6
1691 
${}\ff.\|a\K(\ff.\\{dat}[\|k-\T{1}]\?\|k-\T{1}:\|k);{}$\6
1692 
${}\ff.\|b\K(\ff.\\{dat}[\|k+\T{1}]\?\|k+\T{1}:\|k);{}$\6
1693 
${}\\{gg}.\|a\K(\\{gg}.\\{dat}[\|k-\T{1}]\?\|k-\T{1}:\|k);{}$\6
1694 
${}\\{gg}.\|b\K(\\{gg}.\\{dat}[\|k+\T{1}]\?\|k+\T{1}:\|k){}$;\par
1695 
\U54.\fi
1696 
 
1697 
\M{64}If $e$ is sufficiently small, the fractions $f$ and $g$ will be less
1698 
than~1,
1699 
and we can use the stated algorithm directly. Of course, if $e$ is
1700 
extremely small, a lot of leading zeros need to be lopped off; in the
1701 
worst case, we may have to multiply $f$ and~$g$ by~10 more than 300 times.
1702 
But hey, we don't need to do that extremely often, and computers are
1703 
pretty fast nowadays.
1704 
 
1705 
In the small-exponent case, the computation always terminates before
1706 
$f$ becomes zero, because the interval endpoints are fractions with
1707 
denominator $2^t$ for some $t>50$.
1708 
 
1709 
The invariant relations \PB{$\ff.\\{dat}[\ff.\|a]\I\T{0}$} and \PB{$\\{gg}.%
1710 
\\{dat}[\\{gg}.\|a]\I\T{0}$} are
1711 
not maintained by the computation here, when \PB{$\ff.\|a\K\\{origin}$} or %
1712 
\PB{$\\{gg}.\|a\K\\{origin}$}.
1713 
But no harm is done, because \PB{\\{bignum\_compare}} is not used.
1714 
 
1715 
\Y\B\4\X64:Compute the significant digits \PB{\|s} and decimal exponent \PB{%
1716 
\|e}\X${}\E{}$\6
1717 
\&{if} ${}(\|e>\T{\^401}){}$\1\5
1718 
\X65:Compute the significant digits in the large-exponent case\X\2\6
1719 
\&{else}\5
1720 
${}\{{}$\C{ if \PB{$\|e\Z\T{\^401}$} we have \PB{$\\{gg}.\|a\G\\{origin}$} and %
1721 
\PB{$\\{gg}.\\{dat}[\\{origin}]\Z\T{8}$} }\1\6
1722 
\&{if} ${}(\ff.\|a>\\{origin}){}$\1\5
1723 
${}\ff.\\{dat}[\\{origin}]\K\T{0};{}$\2\6
1724 
\&{for} ${}(\|e\K\T{1},\39\|p\K\|s;{}$ ${}\\{gg}.\|a>\\{origin}\V\ff.\\{dat}[%
1725 
\\{origin}]\E\\{gg}.\\{dat}[\\{origin}];{}$ \,)\5
1726 
${}\{{}$\1\6
1727 
\&{if} ${}(\\{gg}.\|a>\\{origin}){}$\1\5
1728 
${}\|e\MM;{}$\2\6
1729 
\&{else}\1\5
1730 
${}{*}\|p\PP\K\ff.\\{dat}[\\{origin}]+\.{'0'},\39\ff.\\{dat}[\\{origin}]\K%
1731 
\T{0},\39\\{gg}.\\{dat}[\\{origin}]\K\T{0};{}$\2\6
1732 
${}\\{bignum\_times\_ten}({\AND}\ff);{}$\6
1733 
${}\\{bignum\_times\_ten}({\AND}\\{gg});{}$\6
1734 
\4${}\}{}$\2\6
1735 
${}{*}\|p\PP\K((\ff.\\{dat}[\\{origin}]+\T{1}+\\{gg}.\\{dat}[\\{origin}])\GG%
1736 
\T{1})+\.{'0'}{}$;\C{ the middle digit }\6
1737 
\4${}\}{}$\2\6
1738 
${}{*}\|p\K\.{'\\0'}{}$;\C{ terminate the string \PB{\|s} }\par
1739 
\U54.\fi
1740 
 
1741 
\M{65}When \PB{\|e} is large, we use the stated algorithm by considering $f$
1742 
and
1743 
$g$ to be fractions whose denominator is a power of~10.
1744 
 
1745 
An interesting case arises when the number to be converted is
1746 
\Hex{44ada56a4b0835bf}, since the interval turns out to be
1747 
$$ (69999999999999991611392\ \ \dts\ \ 70000000000000000000000).$$
1748 
If this were a closed interval, we could simply give the answer
1749 
\.{7e22}; but the number \.{7e22} actually corresponds to
1750 
\Hex{44ada56a4b0835c0}
1751 
because of the round-to-even rule. Therefore the correct answer is, say,
1752 
\.{6.9999999999999995e22}. This example shows that we need a slightly
1753 
different strategy in the case of open intervals; we cannot simply
1754 
look at the first position in which the endpoints have different
1755 
decimal digits. Therefore we change the invariant relation to $0\le f<g\le 1$,
1756 
when open intervals are involved,
1757 
and we do not terminate the process when $f=0$ or $g=1$.
1758 
 
1759 
\Y\B\4\X65:Compute the significant digits in the large-exponent case\X${}\E{}$\6
1760 
${}\{{}$\5
1761 
\1\&{register} \&{int} \\{open}${}\K\|x.\|l\AND\T{1};{}$\7
1762 
${}\\{tt}.\\{dat}[\\{origin}]\K\T{10};{}$\6
1763 
${}\\{tt}.\|a\K\\{tt}.\|b\K\\{origin};{}$\6
1764 
\&{for} ${}(\|e\K\T{1};{}$ ${}\\{bignum\_compare}({\AND}\\{gg},\39{\AND}\\{tt})%
1765 
\G\\{open};{}$ ${}\|e\PP){}$\1\5
1766 
${}\\{bignum\_times\_ten}({\AND}\\{tt});{}$\2\6
1767 
${}\|p\K\|s;{}$\6
1768 
\&{while} (\T{1})\5
1769 
${}\{{}$\1\6
1770 
${}\\{bignum\_times\_ten}({\AND}\ff);{}$\6
1771 
${}\\{bignum\_times\_ten}({\AND}\\{gg});{}$\6
1772 
\&{for} ${}(\|j\K\.{'0'};{}$ ${}\\{bignum\_compare}({\AND}\ff,\39{\AND}\\{tt})%
1773 
\G\T{0};{}$ ${}\|j\PP){}$\1\5
1774 
${}\\{bignum\_dec}({\AND}\ff,\39{\AND}\\{tt},\39\T{\^10000000}),\39\\{bignum%
1775 
\_dec}({\AND}\\{gg},\39{\AND}\\{tt},\39\T{\^10000000});{}$\2\6
1776 
\&{if} ${}(\\{bignum\_compare}({\AND}\\{gg},\39{\AND}\\{tt})\G\\{open}){}$\1\5
1777 
\&{break};\2\6
1778 
${}{*}\|p\PP\K\|j;{}$\6
1779 
\&{if} ${}(\ff.\|a\E\\{bignum\_prec}-\T{1}\W\R\\{open}){}$\1\5
1780 
\&{goto} \\{done};\C{ $f=0$ in a closed interval }\2\6
1781 
\4${}\}{}$\2\6
1782 
\&{for} ${}(\|k\K\|j;{}$ ${}\\{bignum\_compare}({\AND}\\{gg},\39{\AND}\\{tt})\G%
1783 
\\{open};{}$ ${}\|k\PP){}$\1\5
1784 
${}\\{bignum\_dec}({\AND}\\{gg},\39{\AND}\\{tt},\39\T{\^10000000});{}$\2\6
1785 
${}{*}\|p\PP\K(\|j+\T{1}+\|k)\GG\T{1}{}$;\C{ the middle digit }\6
1786 
\4\\{done}:\5
1787 
;\6
1788 
\4${}\}{}$\2\par
1789 
\U64.\fi
1790 
 
1791 
\M{66}The length of string~\PB{\|s} will be at most 17. For if $f$ and $g$
1792 
agree to 17 places, we have $g/f<1+10^{-16}$; but the
1793 
ratio $g/f$ is always $\ge(1+2^{-52}+2^{-53})/(1+2^{-52}-2^{-53})
1794 
>1+2\times10^{-16}$.
1795 
 
1796 
\Y\B\4\X56:Local variables for \PB{\\{print\_float}}\X${}\mathrel+\E{}$\6
1797 
\&{bignum} ${}\ff,{}$ \\{gg};\C{ fractions or numerators of fractions }\6
1798 
\&{bignum} \\{tt};\C{ power of ten (used as the denominator) }\6
1799 
\&{char} \|s[\T{18}];\6
1800 
\&{register} \&{char} ${}{*}\|p{}$;\par
1801 
\fi
1802 
 
1803 
\M{67}At this point the significant digits are in string \PB{\|s}, and \PB{$%
1804 
\|s[\T{0}]\I\.{'0'}$}.
1805 
If we put a decimal point at the left of~\PB{\|s}, the result should
1806 
be multiplied by $10^e$.
1807 
 
1808 
We prefer the output `\.{300.}' to the form `\.{3e2}', and we prefer
1809 
`\.{.03}' to `\.{3e-2}'. In general, the output will use an
1810 
explicit exponent only if the alternative would take more than
1811 
18~characters.
1812 
 
1813 
\Y\B\4\X67:Print the significant digits with proper context\X${}\E{}$\6
1814 
\&{if} ${}(\|e>\T{17}\V\|e<{}$(\&{int}) \\{strlen}(\|s)${}-\T{17}){}$\1\5
1815 
${}\\{printf}(\.{"\%c\%s\%se\%d"},\39\|s[\T{0}],\39(\|s[\T{1}]\?\.{"."}:%
1816 
\.{""}),\39\|s+\T{1},\39\|e-\T{1});{}$\2\6
1817 
\&{else} \&{if} ${}(\|e<\T{0}){}$\1\5
1818 
${}\\{printf}(\.{".\%0*d\%s"},\39{-}\|e,\39\T{0},\39\|s);{}$\2\6
1819 
\&{else} \&{if} ${}(\\{strlen}(\|s)\G\|e){}$\1\5
1820 
${}\\{printf}(\.{"\%.*s.\%s"},\39\|e,\39\|s,\39\|s+\|e);{}$\2\6
1821 
\&{else}\1\5
1822 
${}\\{printf}(\.{"\%s\%0*d."},\39\|s,\39\|e-{}$(\&{int}) \\{strlen}(\|s)${},\39%
1823 
\T{0}){}$;\2\par
1824 
\U54.\fi
1825 
 
1826 
\N{1}{68}Floating point input conversion. Going the other way, we want to
1827 
be able to convert a given decimal number into its floating binary
1828 
equivalent. The following syntax is supported:
1829 
$$\vbox{\halign{$#$\hfil\cr
1830 
\<digit>\is\.0\mid\.1\mid\.2\mid\.3\mid\.4\mid
1831 
\.5\mid\.6\mid\.7\mid\.8\mid\.9\cr
1832 
\<digit string>\is\<digit>\mid\<digit string>\<digit>\cr
1833 
\<decimal string>\is\<digit string>\..\mid\..\<digit string>\mid
1834 
\<digit string>\..\<digit string>\cr
1835 
\<optional sign>\is\<empty>\mid\.+\mid\.-\cr
1836 
\<exponent>\is\.e\<optional sign>\<digit string>\cr
1837 
\<optional exponent>\is\<empty>\mid\<exponent>\cr
1838 
\<floating magnitude>\is\<digit string>\<exponent>\mid
1839 
\<decimal string>\<optional exponent>\mid\cr
1840 
\hskip12em          \.{Inf}\mid\.{NaN}\mid\.{NaN.}\<digit string>\cr
1841 
\<floating constant>\is\<optional sign>\<floating magnitude>\cr
1842 
\<decimal constant>\is\<optional sign>\<digit string>\cr
1843 
}}$$
1844 
For example, `\.{-3.}' is the floating constant \Hex{c008000000000000}%
1845 
\thinspace;
1846 
`\.{1e3}' and `\.{1000}' are both equivalent to \Hex{408f400000000000}%
1847 
\thinspace;
1848 
`\.{NaN}' and `\.{+NaN.5}' are both equivalent to \Hex{7ff8000000000000}.
1849 
 
1850 
The \PB{\\{scan\_const}} routine looks at a given string and finds the
1851 
longest initial substring that matches the syntax of either \<decimal
1852 
constant> or \<floating constant>. It puts the corresponding value
1853 
into the global octabyte variable~\PB{\\{val}}; it also puts the position of
1854 
the first
1855 
unscanned character in the global pointer variable \PB{\\{next\_char}}.
1856 
It returns 1 if a floating constant was found, 0~if a decimal constant
1857 
was found, $-1$ if nothing was found. A decimal constant that doesn't
1858 
fit in an octabyte is computed modulo~$2^{64}$.
1859 
 
1860 
The value of \PB{\\{exceptions}} set by \PB{\\{scan\_const}} is not necessarily
1861 
correct.
1862 
 
1863 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
1864 
\&{static} \&{void} \\{bignum\_double}\,\,\.{ARGS}((\&{bignum} ${}{*}));{}$\6
1865 
\&{int} \\{scan\_const}\,\,\.{ARGS}((\&{char} ${}{*})){}$;\5
1866 
\hbox{}\6{}\&{int} \\{scan\_const}(\|s)\1\1\6
1867 
\&{char} ${}{*}\|s;\2\2{}$\6
1868 
${}\{{}$\1\6
1869 
\X70:Local variables for \PB{\\{scan\_const}}\X;\6
1870 
${}\\{val}.\|h\K\\{val}.\|l\K\T{0};{}$\6
1871 
${}\|p\K\|s;{}$\6
1872 
\&{if} ${}({*}\|p\E\.{'+'}\V{*}\|p\E\.{'-'}){}$\1\5
1873 
${}\\{sign}\K{*}\|p\PP{}$;\5
1874 
\2\&{else}\1\5
1875 
${}\\{sign}\K\.{'+'};{}$\2\6
1876 
\&{if} ${}(\\{strncmp}(\|p,\39\.{"NaN"},\39\T{3})\E\T{0}){}$\1\5
1877 
${}\\{NaN}\K\\{true},\39\|p\MRL{+{\K}}\T{3};{}$\2\6
1878 
\&{else}\1\5
1879 
${}\\{NaN}\K\\{false};{}$\2\6
1880 
\&{if} ${}((\\{isdigit}({*}\|p)\W\R\\{NaN})\V({*}\|p\E\.{'.'}\W\\{isdigit}({*}(%
1881 
\|p+\T{1})))){}$\1\5
1882 
\X73:Scan a number and \PB{\&{return}}\X;\2\6
1883 
\&{if} (\\{NaN})\1\5
1884 
\X71:Return the standard NaN\X;\2\6
1885 
\&{if} ${}(\\{strncmp}(\|p,\39\.{"Inf"},\39\T{3})\E\T{0}){}$\1\5
1886 
\X72:Return infinity\X;\2\6
1887 
\4\\{no\_const\_found}:\5
1888 
${}\\{next\_char}\K\|s{}$;\5
1889 
\&{return} ${}{-}\T{1};{}$\6
1890 
\4${}\}{}$\2\par
1891 
\fi
1892 
 
1893 
\M{69}\B\X4:Global variables\X${}\mathrel+\E{}$\6
1894 
\&{octa} \\{val};\C{ value returned by \PB{\\{scan\_const}} }\6
1895 
\&{char} ${}{*}\\{next\_char}{}$;\C{ pointer returned by \PB{\\{scan\_const}} }%
1896 
\par
1897 
\fi
1898 
 
1899 
\M{70}\B\X70:Local variables for \PB{\\{scan\_const}}\X${}\E{}$\6
1900 
\&{register} \&{char} ${}{*}\|p,{}$ ${}{*}\|q{}$;\C{ for string manipulations }%
1901 
\6
1902 
\&{register} \&{bool} \\{NaN};\C{ are we processing a NaN? }\6
1903 
\&{int} \\{sign};\C{ \PB{\.{'+'}} or \PB{\.{'-'}} }\par
1904 
\As76\ET81.
1905 
\U68.\fi
1906 
 
1907 
\M{71}\B\X71:Return the standard NaN\X${}\E{}$\6
1908 
${}\{{}$\1\6
1909 
${}\\{next\_char}\K\|p;{}$\6
1910 
${}\\{val}.\|h\K\T{\^600000},\39\\{exp}\K\T{\^3fe};{}$\6
1911 
\&{goto} \\{packit};\6
1912 
\4${}\}{}$\2\par
1913 
\U68.\fi
1914 
 
1915 
\M{72}\B\X72:Return infinity\X${}\E{}$\6
1916 
${}\{{}$\1\6
1917 
${}\\{next\_char}\K\|p+\T{3};{}$\6
1918 
\&{goto} \\{make\_it\_infinite};\6
1919 
\4${}\}{}$\2\par
1920 
\U68.\fi
1921 
 
1922 
\M{73}We saw above that a string of at most 17 digits is enough to characterize
1923 
a floating point number, for purposes of output. But a much longer buffer
1924 
for digits is needed when we're doing input. For example, consider the
1925 
borderline quantity $(1+2^{-53})/2^{1022}$; its decimal expansion, when
1926 
written out exactly, is a number with more than 750 significant digits:
1927 
\.{2.2250738585...8125e-308}.
1928 
If {\it any one\/} of those digits is increased, or if
1929 
additional nonzero digits are added as in
1930 
\.{2.2250738585...81250000001e-308},
1931 
the rounded value is supposed to change from \Hex{0010000000000000}
1932 
to \Hex{0010000000000001}.
1933 
 
1934 
We assume here that the user prefers a perfectly correct answer to
1935 
a speedy almost-correct one, so we implement the most general case.
1936 
 
1937 
\Y\B\4\X73:Scan a number and \PB{\&{return}}\X${}\E{}$\6
1938 
${}\{{}$\1\6
1939 
\&{for} ${}(\|q\K\\{buf0},\39\\{dec\_pt}\K{}$(\&{char} ${}{*}){}$ \T{0}; ${}%
1940 
\\{isdigit}({*}\|p);{}$ ${}\|p\PP){}$\5
1941 
${}\{{}$\1\6
1942 
${}\\{val}\K\\{oplus}(\\{val},\39\\{shift\_left}(\\{val},\39\T{2})){}$;\C{
1943 
multiply by 5 }\6
1944 
${}\\{val}\K\\{incr}(\\{shift\_left}(\\{val},\39\T{1}),\39{*}\|p-\.{'0'});{}$\6
1945 
\&{if} ${}(\|q>\\{buf0}\V{*}\|p\I\.{'0'}){}$\1\6
1946 
\&{if} ${}(\|q<\\{buf\_max}){}$\1\5
1947 
${}{*}\|q\PP\K{*}\|p;{}$\2\6
1948 
\&{else} \&{if} ${}({*}(\|q-\T{1})\E\.{'0'}){}$\1\5
1949 
${}{*}(\|q-\T{1})\K{*}\|p;{}$\2\2\6
1950 
\4${}\}{}$\2\6
1951 
\&{if} (\\{NaN})\1\5
1952 
${}{*}\|q\PP\K\.{'1'};{}$\2\6
1953 
\&{if} ${}({*}\|p\E\.{'.'}){}$\1\5
1954 
\X74:Scan a fraction part\X;\2\6
1955 
${}\\{next\_char}\K\|p;{}$\6
1956 
\&{if} ${}({*}\|p\E\.{'e'}\W\R\\{NaN}){}$\1\5
1957 
\X77:Scan an exponent\X\2\6
1958 
\&{else}\1\5
1959 
${}\\{exp}\K\T{0};{}$\2\6
1960 
\&{if} (\\{dec\_pt})\1\5
1961 
\X78:Return a floating point constant\X;\2\6
1962 
\&{if} ${}(\\{sign}\E\.{'-'}){}$\1\5
1963 
${}\\{val}\K\\{ominus}(\\{zero\_octa},\39\\{val});{}$\2\6
1964 
\&{return} \T{0};\6
1965 
\4${}\}{}$\2\par
1966 
\U68.\fi
1967 
 
1968 
\M{74}\B\X74:Scan a fraction part\X${}\E{}$\6
1969 
${}\{{}$\1\6
1970 
${}\\{dec\_pt}\K\|q;{}$\6
1971 
${}\|p\PP;{}$\6
1972 
\&{for} ${}(\\{zeros}\K\T{0};{}$ ${}\\{isdigit}({*}\|p);{}$ ${}\|p\PP){}$\1\6
1973 
\&{if} ${}({*}\|p\E\.{'0'}\W\|q\E\\{buf0}){}$\1\5
1974 
${}\\{zeros}\PP;{}$\2\6
1975 
\&{else} \&{if} ${}(\|q<\\{buf\_max}){}$\1\5
1976 
${}{*}\|q\PP\K{*}\|p;{}$\2\6
1977 
\&{else} \&{if} ${}({*}(\|q-\T{1})\E\.{'0'}){}$\1\5
1978 
${}{*}(\|q-\T{1})\K{*}\|p;{}$\2\2\6
1979 
\4${}\}{}$\2\par
1980 
\U73.\fi
1981 
 
1982 
\M{75}The buffer needs room for eight digits of padding at the left, followed
1983 
by up to $1022+53-307$ significant digits, followed by a ``sticky'' digit
1984 
at position \PB{$\\{buf\_max}-\T{1}$}, and eight more digits of padding.
1985 
 
1986 
\Y\B\4\D$\\{buf0}$ \5
1987 
$(\\{buf}+\T{8}{}$)\par
1988 
\B\4\D$\\{buf\_max}$ \5
1989 
$(\\{buf}+\T{777}{}$)\par
1990 
\Y\B\4\X4:Global variables\X${}\mathrel+\E{}$\6
1991 
\&{static} \&{char} \\{buf}[\T{785}]${}\K\.{"00000000"}{}$;\C{ where we put
1992 
significant input digits }\par
1993 
\fi
1994 
 
1995 
\M{76}\B\X70:Local variables for \PB{\\{scan\_const}}\X${}\mathrel+\E{}$\6
1996 
\&{register} \&{char} ${}{*}\\{dec\_pt}{}$;\C{ position of decimal point in %
1997 
\PB{\\{buf}} }\6
1998 
\&{register} \&{int} \\{exp};\C{ scanned exponent; later used for raw binary
1999 
exponent }\6
2000 
\&{register} \&{int} \\{zeros};\C{ leading zeros removed after decimal point }%
2001 
\par
2002 
\fi
2003 
 
2004 
\M{77}Here we don't advance \PB{\\{next\_char}} and force a decimal point until
2005 
we
2006 
know that a syntactically correct exponent exists.
2007 
 
2008 
The code here will convert extra-large inputs like
2009 
`\.{9e+9999999999999999}' into $\infty$ and extra-small inputs into zero.
2010 
Strange inputs like `\.{-00.0e9999999}' must also be accommodated.
2011 
 
2012 
\Y\B\4\X77:Scan an exponent\X${}\E{}$\6
2013 
${}\{{}$\5
2014 
\1\&{register} \&{char} \\{exp\_sign};\7
2015 
${}\|p\PP;{}$\6
2016 
\&{if} ${}({*}\|p\E\.{'+'}\V{*}\|p\E\.{'-'}){}$\1\5
2017 
${}\\{exp\_sign}\K{*}\|p\PP{}$;\5
2018 
\2\&{else}\1\5
2019 
${}\\{exp\_sign}\K\.{'+'};{}$\2\6
2020 
\&{if} ${}(\\{isdigit}({*}\|p)){}$\5
2021 
${}\{{}$\1\6
2022 
\&{for} ${}(\\{exp}\K{*}\|p\PP-\.{'0'};{}$ ${}\\{isdigit}({*}\|p);{}$ ${}\|p%
2023 
\PP){}$\1\6
2024 
\&{if} ${}(\\{exp}<\T{1000}){}$\1\5
2025 
${}\\{exp}\K\T{10}*\\{exp}+{*}\|p-\.{'0'};{}$\2\2\6
2026 
\&{if} ${}(\R\\{dec\_pt}){}$\1\5
2027 
${}\\{dec\_pt}\K\|q,\39\\{zeros}\K\T{0};{}$\2\6
2028 
\&{if} ${}(\\{exp\_sign}\E\.{'-'}){}$\1\5
2029 
${}\\{exp}\K{-}\\{exp};{}$\2\6
2030 
${}\\{next\_char}\K\|p;{}$\6
2031 
\4${}\}{}$\2\6
2032 
\4${}\}{}$\2\par
2033 
\U73.\fi
2034 
 
2035 
\M{78}\B\X78:Return a floating point constant\X${}\E{}$\6
2036 
${}\{{}$\1\6
2037 
\X79:Move the digits from \PB{\\{buf}} to \PB{$\ff$}\X;\6
2038 
\X83:Determine the binary fraction and binary exponent\X;\6
2039 
\4\\{packit}:\5
2040 
\X84:Pack and round the answer\X;\6
2041 
\&{return} \T{1};\6
2042 
\4${}\}{}$\2\par
2043 
\U73.\fi
2044 
 
2045 
\M{79}Now we get ready to compute the binary fraction bits, by putting the
2046 
scanned input digits into a multiprecision fixed-point
2047 
accumulator \PB{$\ff$} that spans the full necessary range.
2048 
After this step, the number that we want to convert to floating binary
2049 
will appear in \PB{$\ff.\\{dat}[\ff.\|a]$}, \PB{$\ff.\\{dat}[\ff.\|a+\T{1}]$}, %
2050 
\dots,
2051 
\PB{$\ff.\\{dat}[\ff.\|b]$}.
2052 
The radix-$10^9$ digit in ${\it ff}[36-k]$ is understood to be multiplied
2053 
by $10^{9k}$, for $36\ge k\ge-120$.
2054 
 
2055 
\Y\B\4\X79:Move the digits from \PB{\\{buf}} to \PB{$\ff$}\X${}\E{}$\6
2056 
$\|x\K\\{buf}+\T{341}+\\{zeros}-\\{dec\_pt}-\\{exp};{}$\6
2057 
\&{if} ${}(\|q\E\\{buf0}\V\|x\G\T{1413}){}$\5
2058 
${}\{{}$\1\6
2059 
\4\\{make\_it\_zero}:\5
2060 
${}\\{exp}\K{-}\T{99999}{}$;\5
2061 
\&{goto} \\{packit};\6
2062 
\4${}\}{}$\2\6
2063 
\&{if} ${}(\|x<\T{0}){}$\5
2064 
${}\{{}$\1\6
2065 
\4\\{make\_it\_infinite}:\5
2066 
${}\\{exp}\K\T{99999}{}$;\5
2067 
\&{goto} \\{packit};\6
2068 
\4${}\}{}$\2\6
2069 
${}\ff.\|a\K\|x/\T{9};{}$\6
2070 
\&{for} ${}(\|p\K\|q;{}$ ${}\|p<\|q+\T{8};{}$ ${}\|p\PP){}$\1\5
2071 
${}{*}\|p\K\.{'0'}{}$;\C{ pad with trailing zeros }\2\6
2072 
${}\|q\K\|q-\T{1}-(\|q+\T{341}+\\{zeros}-\\{dec\_pt}-\\{exp})\MOD\T{9}{}$;\C{
2073 
compute stopping place in \PB{\\{buf}} }\6
2074 
\&{for} ${}(\|p\K\\{buf0}-\|x\MOD\T{9},\39\|k\K\ff.\|a;{}$ ${}\|p\Z\|q\W\|k\Z%
2075 
\T{156};{}$ ${}\|p\MRL{+{\K}}\T{9},\39\|k\PP){}$\1\5
2076 
\X80:Put the 9-digit number \PB{${*}\|p$}\thinspace\dots\thinspace\PB{${*}(\|p+%
2077 
\T{8})$} into \PB{$\ff.\\{dat}[\|k]$}\X;\2\6
2078 
${}\ff.\|b\K\|k-\T{1};{}$\6
2079 
\&{for} ${}(\|x\K\T{0};{}$ ${}\|p\Z\|q;{}$ ${}\|p\MRL{+{\K}}\T{9}){}$\1\6
2080 
\&{if} ${}(\\{strncmp}(\|p,\39\.{"000000000"},\39\T{9})\I\T{0}){}$\1\5
2081 
${}\|x\K\T{1};{}$\2\2\6
2082 
${}\ff.\\{dat}[\T{156}]\MRL{+{\K}}\|x{}$;\C{ nonzero digits that fall off the
2083 
right are sticky }\6
2084 
\&{while} ${}(\ff.\\{dat}[\ff.\|b]\E\T{0}){}$\1\5
2085 
${}\ff.\|b\MM{}$;\2\par
2086 
\U78.\fi
2087 
 
2088 
\M{80}\B\X80:Put the 9-digit number \PB{${*}\|p$}\thinspace\dots\thinspace%
2089 
\PB{${*}(\|p+\T{8})$} into \PB{$\ff.\\{dat}[\|k]$}\X${}\E{}$\6
2090 
${}\{{}$\1\6
2091 
\&{for} ${}(\|x\K{*}\|p-\.{'0'},\39\\{pp}\K\|p+\T{1};{}$ ${}\\{pp}<\|p+%
2092 
\T{9};{}$ ${}\\{pp}\PP){}$\1\5
2093 
${}\|x\K\T{10}*\|x+{*}\\{pp}-\.{'0'};{}$\2\6
2094 
${}\ff.\\{dat}[\|k]\K\|x;{}$\6
2095 
\4${}\}{}$\2\par
2096 
\U79.\fi
2097 
 
2098 
\M{81}\B\X70:Local variables for \PB{\\{scan\_const}}\X${}\mathrel+\E{}$\6
2099 
\&{register} \&{int} \|k${},{}$ \|x;\6
2100 
\&{register} \&{char} ${}{*}\\{pp};{}$\6
2101 
\&{bignum} ${}\ff,{}$ \\{tt};\par
2102 
\fi
2103 
 
2104 
\M{82}Here's a subroutine that is dual to \PB{\\{bignum\_times\_ten}}. It
2105 
changes $f$
2106 
to~$2f$, assuming that overflow will not occur and that the radix is $10^9$.
2107 
 
2108 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
2109 
\&{static} \&{void} \\{bignum\_double}(\|f)\1\1\6
2110 
\&{bignum} ${}{*}\|f;\2\2{}$\6
2111 
${}\{{}$\1\6
2112 
\&{register} \&{tetra} ${}{*}\|p,{}$ ${}{*}\|q;{}$\6
2113 
\&{register} \&{int} \|x${},{}$ \\{carry};\7
2114 
\&{for} ${}(\|p\K{\AND}\|f\MG\\{dat}[\|f\MG\|b],\39\|q\K{\AND}\|f\MG\\{dat}[\|f%
2115 
\MG\|a],\39\\{carry}\K\T{0};{}$ ${}\|p\G\|q;{}$ ${}\|p\MM){}$\5
2116 
${}\{{}$\1\6
2117 
${}\|x\K{*}\|p+{*}\|p+\\{carry};{}$\6
2118 
\&{if} ${}(\|x\G\T{1000000000}){}$\1\5
2119 
${}\\{carry}\K\T{1},\39{*}\|p\K\|x-\T{1000000000};{}$\2\6
2120 
\&{else}\1\5
2121 
${}\\{carry}\K\T{0},\39{*}\|p\K\|x;{}$\2\6
2122 
\4${}\}{}$\2\6
2123 
${}{*}\|p\K\\{carry};{}$\6
2124 
\&{if} (\\{carry})\1\5
2125 
${}\|f\MG\|a\MM;{}$\2\6
2126 
\&{if} ${}(\|f\MG\\{dat}[\|f\MG\|b]\E\T{0}\W\|f\MG\|b>\|f\MG\|a){}$\1\5
2127 
${}\|f\MG\|b\MM;{}$\2\6
2128 
\4${}\}{}$\2\par
2129 
\fi
2130 
 
2131 
\M{83}\B\X83:Determine the binary fraction and binary exponent\X${}\E{}$\6
2132 
$\\{val}\K\\{zero\_octa};{}$\6
2133 
\&{if} ${}(\ff.\|a>\T{36}){}$\5
2134 
${}\{{}$\1\6
2135 
\&{for} ${}(\\{exp}\K\T{\^3fe};{}$ ${}\ff.\|a>\T{36};{}$ ${}\\{exp}\MM){}$\1\5
2136 
${}\\{bignum\_double}({\AND}\ff);{}$\2\6
2137 
\&{for} ${}(\|k\K\T{54};{}$ \|k; ${}\|k\MM){}$\5
2138 
${}\{{}$\1\6
2139 
\&{if} ${}(\ff.\\{dat}[\T{36}]){}$\5
2140 
${}\{{}$\1\6
2141 
\&{if} ${}(\|k\G\T{32}){}$\1\5
2142 
${}\\{val}.\|h\MRL{{\OR}{\K}}\T{1}\LL(\|k-\T{32}){}$;\5
2143 
\2\&{else}\1\5
2144 
${}\\{val}.\|l\MRL{{\OR}{\K}}\T{1}\LL\|k;{}$\2\6
2145 
${}\ff.\\{dat}[\T{36}]\K\T{0};{}$\6
2146 
\&{if} ${}(\ff.\|b\E\T{36}){}$\1\5
2147 
\&{break};\C{ break if \PB{$\ff$} now zero }\2\6
2148 
\4${}\}{}$\2\6
2149 
${}\\{bignum\_double}({\AND}\ff);{}$\6
2150 
\4${}\}{}$\2\6
2151 
\4${}\}{}$\5
2152 
\2\&{else}\5
2153 
${}\{{}$\1\6
2154 
${}\\{tt}.\|a\K\\{tt}.\|b\K\T{36},\39\\{tt}.\\{dat}[\T{36}]\K\T{2};{}$\6
2155 
\&{for} ${}(\\{exp}\K\T{\^3fe};{}$ ${}\\{bignum\_compare}({\AND}\ff,\39{\AND}%
2156 
\\{tt})\G\T{0};{}$ ${}\\{exp}\PP){}$\1\5
2157 
${}\\{bignum\_double}({\AND}\\{tt});{}$\2\6
2158 
\&{for} ${}(\|k\K\T{54};{}$ \|k; ${}\|k\MM){}$\5
2159 
${}\{{}$\1\6
2160 
${}\\{bignum\_double}({\AND}\ff);{}$\6
2161 
\&{if} ${}(\\{bignum\_compare}({\AND}\ff,\39{\AND}\\{tt})\G\T{0}){}$\5
2162 
${}\{{}$\1\6
2163 
\&{if} ${}(\|k\G\T{32}){}$\1\5
2164 
${}\\{val}.\|h\MRL{{\OR}{\K}}\T{1}\LL(\|k-\T{32}){}$;\5
2165 
\2\&{else}\1\5
2166 
${}\\{val}.\|l\MRL{{\OR}{\K}}\T{1}\LL\|k;{}$\2\6
2167 
${}\\{bignum\_dec}({\AND}\ff,\39{\AND}\\{tt},\39\T{1000000000});{}$\6
2168 
\&{if} ${}(\ff.\|a\E\\{bignum\_prec}-\T{1}){}$\1\5
2169 
\&{break};\C{ break if \PB{$\ff$} now zero }\2\6
2170 
\4${}\}{}$\2\6
2171 
\4${}\}{}$\2\6
2172 
\4${}\}{}$\2\6
2173 
\&{if} ${}(\|k\E\T{0}){}$\1\5
2174 
${}\\{val}.\|l\MRL{{\OR}{\K}}\T{1}{}$;\C{ add sticky bit if \PB{$\ff$} nonzero
2175 
}\2\par
2176 
\U78.\fi
2177 
 
2178 
\M{84}We need to be careful that the input `\.{NaN.999999999999999999999}'
2179 
doesn't
2180 
get rounded up; it is supposed to yield \Hex{7fffffffffffffff}.
2181 
 
2182 
Although the input `\.{NaN.0}' is illegal, strictly speaking, we silently
2183 
convert it to \Hex{7ff0000000000001}---a number that would be
2184 
output as `\.{NaN.0000000000000002}'.
2185 
 
2186 
\Y\B\4\X84:Pack and round the answer\X${}\E{}$\6
2187 
$\\{val}\K\\{fpack}(\\{val},\39\\{exp},\39\\{sign},\39\.{ROUND\_NEAR});{}$\6
2188 
\&{if} (\\{NaN})\5
2189 
${}\{{}$\1\6
2190 
\&{if} ${}((\\{val}.\|h\AND\T{\^7fffffff})\E\T{\^40000000}){}$\1\5
2191 
${}\\{val}.\|h\MRL{{\OR}{\K}}\T{\^7fffffff},\39\\{val}.\|l\K\T{\^ffffffff};{}$%
2192 
\2\6
2193 
\&{else} \&{if} ${}((\\{val}.\|h\AND\T{\^7fffffff})\E\T{\^3ff00000}\W\R\\{val}.%
2194 
\|l){}$\1\5
2195 
${}\\{val}.\|h\MRL{{\OR}{\K}}\T{\^40000000},\39\\{val}.\|l\K\T{1};{}$\2\6
2196 
\&{else}\1\5
2197 
${}\\{val}.\|h\MRL{{\OR}{\K}}\T{\^40000000};{}$\2\6
2198 
\4${}\}{}$\2\par
2199 
\U78.\fi
2200 
 
2201 
\N{1}{85}Floating point remainders. In this section we implement the remainder
2202 
of the floating point operations---one of which happens to be the
2203 
operation of taking the remainder.
2204 
 
2205 
The easiest task remaining is to compare two floating point quantities.
2206 
Routine \PB{\\{fcomp}} returns $-1$~if~$y<z$, 0~if~$y=z$, $+1$~if~$y>z$, and
2207 
$+2$~if $y$ and~$z$ are unordered.
2208 
 
2209 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
2210 
\&{int} \\{fcomp}\,\,${}\.{ARGS}((\&{octa},\39\&{octa})){}$;\5
2211 
\hbox{}\6{}\&{int} ${}\\{fcomp}(\|y,\39\|z){}$\1\1\6
2212 
\&{octa} \|y${},{}$ \|z;\2\2\6
2213 
${}\{{}$\1\6
2214 
\&{ftype} \\{yt}${},{}$ \\{zt};\6
2215 
\&{int} \\{ye}${},{}$ \\{ze};\6
2216 
\&{char} \\{ys}${},{}$ \\{zs};\6
2217 
\&{octa} \\{yf}${},{}$ \\{zf};\6
2218 
\&{register} \&{int} \|x;\7
2219 
${}\\{yt}\K\\{funpack}(\|y,\39{\AND}\\{yf},\39{\AND}\\{ye},\39{\AND}\\{ys});{}$%
2220 
\6
2221 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
2222 
\6
2223 
\&{switch} ${}(\T{4}*\\{yt}+\\{zt}){}$\5
2224 
${}\{{}$\1\6
2225 
\4\&{case} \T{4}${}*\\{nan}+\\{nan}{}$:\5
2226 
\&{case} \T{4}${}*\\{zro}+\\{nan}{}$:\5
2227 
\&{case} \T{4}${}*\\{num}+\\{nan}{}$:\5
2228 
\&{case} \T{4}${}*\\{inf}+\\{nan}{}$:\5
2229 
\&{case} \T{4}${}*\\{nan}+\\{zro}{}$:\5
2230 
\&{case} \T{4}${}*\\{nan}+\\{num}{}$:\5
2231 
\&{case} \T{4}${}*\\{nan}+\\{inf}{}$:\5
2232 
\&{return} \T{2};\6
2233 
\4\&{case} \T{4}${}*\\{zro}+\\{zro}{}$:\5
2234 
\&{return} \T{0};\6
2235 
\4\&{case} \T{4}${}*\\{zro}+\\{num}{}$:\5
2236 
\&{case} \T{4}${}*\\{num}+\\{zro}{}$:\5
2237 
\&{case} \T{4}${}*\\{zro}+\\{inf}{}$:\5
2238 
\&{case} \T{4}${}*\\{inf}+\\{zro}{}$:\5
2239 
\&{case} \T{4}${}*\\{num}+\\{num}{}$:\5
2240 
\&{case} \T{4}${}*\\{num}+\\{inf}{}$:\5
2241 
\&{case} \T{4}${}*\\{inf}+\\{num}{}$:\5
2242 
\&{case} \T{4}${}*\\{inf}+\\{inf}{}$:\6
2243 
\&{if} ${}(\\{ys}\I\\{zs}){}$\1\5
2244 
${}\|x\K\T{1};{}$\2\6
2245 
\&{else} \&{if} ${}(\|y.\|h>\|z.\|h){}$\1\5
2246 
${}\|x\K\T{1};{}$\2\6
2247 
\&{else} \&{if} ${}(\|y.\|h<\|z.\|h){}$\1\5
2248 
${}\|x\K{-}\T{1};{}$\2\6
2249 
\&{else} \&{if} ${}(\|y.\|l>\|z.\|l){}$\1\5
2250 
${}\|x\K\T{1};{}$\2\6
2251 
\&{else} \&{if} ${}(\|y.\|l<\|z.\|l){}$\1\5
2252 
${}\|x\K{-}\T{1};{}$\2\6
2253 
\&{else}\1\5
2254 
\&{return} \T{0};\2\6
2255 
\&{break};\6
2256 
\4${}\}{}$\2\6
2257 
\&{return} ${}(\\{ys}\E\.{'-'}\?{-}\|x:\|x);{}$\6
2258 
\4${}\}{}$\2\par
2259 
\fi
2260 
 
2261 
\M{86}Several \MMIX\ operations act on a single floating point number and
2262 
accept an arbitrary rounding mode. For example, consider the
2263 
operation of rounding to the nearest floating point integer:
2264 
 
2265 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
2266 
\&{octa} \\{fintegerize}\,\,${}\.{ARGS}((\&{octa},\39\&{int})){}$;\5
2267 
\hbox{}\6{}\&{octa} ${}\\{fintegerize}(\|z,\39\|r){}$\1\1\6
2268 
\&{octa} \|z;\C{ the operand }\6
2269 
\&{int} \|r;\C{ the rounding mode }\2\2\6
2270 
${}\{{}$\1\6
2271 
\&{ftype} \\{zt};\6
2272 
\&{int} \\{ze};\6
2273 
\&{char} \\{zs};\6
2274 
\&{octa} \\{xf}${},{}$ \\{zf};\7
2275 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
2276 
\6
2277 
\&{if} ${}(\R\|r){}$\1\5
2278 
${}\|r\K\\{cur\_round};{}$\2\6
2279 
\&{switch} (\\{zt})\5
2280 
${}\{{}$\1\6
2281 
\4\&{case} \\{nan}:\5
2282 
\&{if} ${}(\R(\|z.\|h\AND\T{\^80000})){}$\5
2283 
${}\{{}$\5
2284 
\1${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\5
2285 
${}\|z.\|h\MRL{{\OR}{\K}}\T{\^80000}{}$;\5
2286 
${}\}{}$\2\6
2287 
\4\&{case} \\{inf}:\5
2288 
\&{case} \\{zro}:\5
2289 
\&{return} \|z;\6
2290 
\4\&{case} \\{num}:\5
2291 
\X87:Integerize and \PB{\&{return}}\X;\6
2292 
\4${}\}{}$\2\6
2293 
\4${}\}{}$\2\par
2294 
\fi
2295 
 
2296 
\M{87}\B\X87:Integerize and \PB{\&{return}}\X${}\E{}$\6
2297 
\&{if} ${}(\\{ze}\G\T{1074}){}$\1\5
2298 
\&{return} \\{fpack}${}(\\{zf},\39\\{ze},\39\\{zs},\39\.{ROUND\_OFF}){}$;\C{
2299 
already an integer }\2\6
2300 
\&{if} ${}(\\{ze}\Z\T{1020}){}$\1\5
2301 
${}\\{xf}.\|h\K\T{0},\39\\{xf}.\|l\K\T{1};{}$\2\6
2302 
\&{else}\5
2303 
${}\{{}$\5
2304 
\1\&{octa} \\{oo};\7
2305 
${}\\{xf}\K\\{shift\_right}(\\{zf},\39\T{1074}-\\{ze},\39\T{1});{}$\6
2306 
${}\\{oo}\K\\{shift\_left}(\\{xf},\39\T{1074}-\\{ze});{}$\6
2307 
\&{if} ${}(\\{oo}.\|l\I\\{zf}.\|l\V\\{oo}.\|h\I\\{zf}.\|h){}$\1\5
2308 
${}\\{xf}.\|l\MRL{{\OR}{\K}}\T{1}{}$;\C{ sticky bit }\2\6
2309 
\4${}\}{}$\2\6
2310 
\&{switch} (\|r)\5
2311 
${}\{{}$\1\6
2312 
\4\&{case} \.{ROUND\_DOWN}:\5
2313 
\&{if} ${}(\\{zs}\E\.{'-'}){}$\1\5
2314 
${}\\{xf}\K\\{incr}(\\{xf},\39\T{3}){}$;\5
2315 
\2\&{break};\6
2316 
\4\&{case} \.{ROUND\_UP}:\5
2317 
\&{if} ${}(\\{zs}\I\.{'-'}){}$\1\5
2318 
${}\\{xf}\K\\{incr}(\\{xf},\39\T{3});{}$\2\6
2319 
\4\&{case} \.{ROUND\_OFF}:\5
2320 
\&{break};\6
2321 
\4\&{case} \.{ROUND\_NEAR}:\5
2322 
${}\\{xf}\K\\{incr}(\\{xf},\39\\{xf}.\|l\AND\T{4}\?\T{2}:\T{1}){}$;\5
2323 
\&{break};\6
2324 
\4${}\}{}$\2\6
2325 
${}\\{xf}.\|l\MRL{\AND{\K}}\T{\^fffffffc};{}$\6
2326 
\&{if} ${}(\\{ze}\G\T{1022}){}$\1\5
2327 
\&{return} \\{fpack}${}(\\{shift\_left}(\\{xf},\39\T{1074}-\\{ze}),\39\\{ze},%
2328 
\39\\{zs},\39\.{ROUND\_OFF});{}$\2\6
2329 
\&{if} ${}(\\{xf}.\|l){}$\1\5
2330 
${}\\{xf}.\|h\K\T{\^3ff00000},\39\\{xf}.\|l\K\T{0};{}$\2\6
2331 
\&{if} ${}(\\{zs}\E\.{'-'}){}$\1\5
2332 
${}\\{xf}.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
2333 
\&{return} \\{xf};\par
2334 
\U86.\fi
2335 
 
2336 
\M{88}To convert floating point to fixed point, we use \PB{\\{fixit}}.
2337 
 
2338 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
2339 
\&{octa} \\{fixit}\,\,${}\.{ARGS}((\&{octa},\39\&{int})){}$;\5
2340 
\hbox{}\6{}\&{octa} ${}\\{fixit}(\|z,\39\|r){}$\1\1\6
2341 
\&{octa} \|z;\C{ the operand }\6
2342 
\&{int} \|r;\C{ the rounding mode }\2\2\6
2343 
${}\{{}$\1\6
2344 
\&{ftype} \\{zt};\6
2345 
\&{int} \\{ze};\6
2346 
\&{char} \\{zs};\6
2347 
\&{octa} \\{zf}${},{}$ \|o;\7
2348 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
2349 
\6
2350 
\&{if} ${}(\R\|r){}$\1\5
2351 
${}\|r\K\\{cur\_round};{}$\2\6
2352 
\&{switch} (\\{zt})\5
2353 
${}\{{}$\1\6
2354 
\4\&{case} \\{nan}:\5
2355 
\&{case} \\{inf}:\5
2356 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\5
2357 
\&{return} \|z;\6
2358 
\4\&{case} \\{zro}:\5
2359 
\&{return} \\{zero\_octa};\6
2360 
\4\&{case} \\{num}:\5
2361 
\&{if} ${}(\\{funpack}(\\{fintegerize}(\|z,\39\|r),\39{\AND}\\{zf},\39{\AND}%
2362 
\\{ze},\39{\AND}\\{zs})\E\\{zro}){}$\1\5
2363 
\&{return} \\{zero\_octa};\2\6
2364 
\&{if} ${}(\\{ze}\Z\T{1076}){}$\1\5
2365 
${}\|o\K\\{shift\_right}(\\{zf},\39\T{1076}-\\{ze},\39\T{1});{}$\2\6
2366 
\&{else}\5
2367 
${}\{{}$\1\6
2368 
\&{if} ${}(\\{ze}>\T{1085}\V(\\{ze}\E\T{1085}\W(\\{zf}.\|h>\T{\^400000}\V%
2369 
\3{-1}(\\{zf}.\|h\E\T{\^400000}\W(\\{zf}.\|l\V\\{zs}\I\.{'-'}))))){}$\1\5
2370 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{W\_BIT};{}$\2\6
2371 
\&{if} ${}(\\{ze}\G\T{1140}){}$\1\5
2372 
\&{return} \\{zero\_octa};\2\6
2373 
${}\|o\K\\{shift\_left}(\\{zf},\39\\{ze}-\T{1076});{}$\6
2374 
\4${}\}{}$\2\6
2375 
\&{return} ${}(\\{zs}\E\.{'-'}\?\\{ominus}(\\{zero\_octa},\39\|o):\|o);{}$\6
2376 
\4${}\}{}$\2\6
2377 
\4${}\}{}$\2\par
2378 
\fi
2379 
 
2380 
\M{89}Going the other way, we can specify not only a rounding mode but whether
2381 
the given fixed point octabyte is signed or unsigned, and whether the
2382 
result should be rounded to short precision.
2383 
 
2384 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
2385 
\&{octa} \\{floatit}\,\,${}\.{ARGS}((\&{octa},\39\&{int},\39\&{int},\39%
2386 
\&{int})){}$;\5
2387 
\hbox{}\6{}\&{octa} ${}\\{floatit}(\|z,\39\|r,\39\|u,\39\|p){}$\1\1\6
2388 
\&{octa} \|z;\C{ octabyte to float }\6
2389 
\&{int} \|r;\C{ rounding mode }\6
2390 
\&{int} \|u;\C{ unsigned? }\6
2391 
\&{int} \|p;\C{ short precision? }\2\2\6
2392 
${}\{{}$\1\6
2393 
\&{int} \|e;\5
2394 
\&{char} \|s;\6
2395 
\&{register} \&{int} \|t;\7
2396 
${}\\{exceptions}\K\T{0};{}$\6
2397 
\&{if} ${}(\R\|z.\|h\W\R\|z.\|l){}$\1\5
2398 
\&{return} \\{zero\_octa};\2\6
2399 
\&{if} ${}(\R\|r){}$\1\5
2400 
${}\|r\K\\{cur\_round};{}$\2\6
2401 
\&{if} ${}(\R\|u\W(\|z.\|h\AND\\{sign\_bit})){}$\1\5
2402 
${}\|s\K\.{'-'},\39\|z\K\\{ominus}(\\{zero\_octa},\39\|z){}$;\5
2403 
\2\&{else}\1\5
2404 
${}\|s\K\.{'+'};{}$\2\6
2405 
${}\|e\K\T{1076};{}$\6
2406 
\&{while} ${}(\|z.\|h<\T{\^400000}){}$\1\5
2407 
${}\|e\MM,\39\|z\K\\{shift\_left}(\|z,\39\T{1});{}$\2\6
2408 
\&{while} ${}(\|z.\|h\G\T{\^800000}){}$\5
2409 
${}\{{}$\1\6
2410 
${}\|e\PP;{}$\6
2411 
${}\|t\K\|z.\|l\AND\T{1};{}$\6
2412 
${}\|z\K\\{shift\_right}(\|z,\39\T{1},\39\T{1});{}$\6
2413 
${}\|z.\|l\MRL{{\OR}{\K}}\|t;{}$\6
2414 
\4${}\}{}$\2\6
2415 
\&{if} (\|p)\1\5
2416 
\X90:Convert to short float\X;\2\6
2417 
\&{return} \\{fpack}${}(\|z,\39\|e,\39\|s,\39\|r);{}$\6
2418 
\4${}\}{}$\2\par
2419 
\fi
2420 
 
2421 
\M{90}\B\X90:Convert to short float\X${}\E{}$\6
2422 
${}\{{}$\1\6
2423 
\&{register} \&{int} \\{ex};\5
2424 
\&{register} \&{tetra} \|t;\7
2425 
${}\|t\K\\{sfpack}(\|z,\39\|e,\39\|s,\39\|r);{}$\6
2426 
${}\\{ex}\K\\{exceptions};{}$\6
2427 
${}\\{sfunpack}(\|t,\39{\AND}\|z,\39{\AND}\|e,\39{\AND}\|s);{}$\6
2428 
${}\\{exceptions}\K\\{ex};{}$\6
2429 
\4${}\}{}$\2\par
2430 
\U89.\fi
2431 
 
2432 
\M{91}The square root operation is more interesting.
2433 
 
2434 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
2435 
\&{octa} \\{froot}\,\,${}\.{ARGS}((\&{octa},\39\&{int})){}$;\5
2436 
\hbox{}\6{}\&{octa} ${}\\{froot}(\|z,\39\|r){}$\1\1\6
2437 
\&{octa} \|z;\C{ the operand }\6
2438 
\&{int} \|r;\C{ the rounding mode }\2\2\6
2439 
${}\{{}$\1\6
2440 
\&{ftype} \\{zt};\6
2441 
\&{int} \\{ze};\6
2442 
\&{char} \\{zs};\6
2443 
\&{octa} \|x${},{}$ \\{xf}${},{}$ \\{rf}${},{}$ \\{zf};\6
2444 
\&{register} \&{int} \\{xe}${},{}$ \|k;\7
2445 
\&{if} ${}(\R\|r){}$\1\5
2446 
${}\|r\K\\{cur\_round};{}$\2\6
2447 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
2448 
\6
2449 
\&{if} ${}(\\{zs}\E\.{'-'}\W\\{zt}\I\\{zro}){}$\1\5
2450 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT},\39\|x\K\\{standard\_NaN};{}$\2\6
2451 
\&{else}\5
2452 
\1\&{switch} (\\{zt})\5
2453 
${}\{{}$\1\6
2454 
\4\&{case} \\{nan}:\5
2455 
\&{if} ${}(\R(\|z.\|h\AND\T{\^80000})){}$\1\5
2456 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT},\39\|z.\|h\MRL{{\OR}{\K}}\T{%
2457 
\^80000};{}$\2\6
2458 
\&{return} \|z;\6
2459 
\4\&{case} \\{inf}:\5
2460 
\&{case} \\{zro}:\5
2461 
${}\|x\K\|z{}$;\5
2462 
\&{break};\6
2463 
\4\&{case} \\{num}:\5
2464 
\X92:Take the square root and \PB{\&{return}}\X;\6
2465 
\4${}\}{}$\2\2\6
2466 
\&{if} ${}(\\{zs}\E\.{'-'}){}$\1\5
2467 
${}\|x.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
2468 
\&{return} \|x;\6
2469 
\4${}\}{}$\2\par
2470 
\fi
2471 
 
2472 
\M{92}The square root can be found by an adaptation of the old pencil-and-paper
2473 
method. If $n=\lfloor\sqrt s\rfloor$, where $s$ is an integer,
2474 
we have $s=n^2+r$ where $0\le r\le2n$;
2475 
this invariant can be maintained if we replace $s$ by $4s+(0,1,2,3)$
2476 
and $n$ by $2n+(0,1)$. The following code implements this idea with
2477 
$2n$ in~\PB{\\{xf}} and $r$ in~\PB{\\{rf}}. (It could easily be made to run
2478 
about
2479 
twice as fast.)
2480 
 
2481 
\Y\B\4\X92:Take the square root and \PB{\&{return}}\X${}\E{}$\6
2482 
$\\{xf}.\|h\K\T{0},\39\\{xf}.\|l\K\T{2};{}$\6
2483 
${}\\{xe}\K(\\{ze}+\T{\^3fe})\GG\T{1};{}$\6
2484 
\&{if} ${}(\\{ze}\AND\T{1}){}$\1\5
2485 
${}\\{zf}\K\\{shift\_left}(\\{zf},\39\T{1});{}$\2\6
2486 
${}\\{rf}.\|h\K\T{0},\39\\{rf}.\|l\K(\\{zf}.\|h\GG\T{22})-\T{1};{}$\6
2487 
\&{for} ${}(\|k\K\T{53};{}$ \|k; ${}\|k\MM){}$\5
2488 
${}\{{}$\1\6
2489 
${}\\{rf}\K\\{shift\_left}(\\{rf},\39\T{2}){}$;\5
2490 
${}\\{xf}\K\\{shift\_left}(\\{xf},\39\T{1});{}$\6
2491 
\&{if} ${}(\|k\G\T{43}){}$\1\5
2492 
${}\\{rf}\K\\{incr}(\\{rf},\39(\\{zf}.\|h\GG(\T{2}*(\|k-\T{43})))\AND\T{3});{}$%
2493 
\2\6
2494 
\&{else} \&{if} ${}(\|k\G\T{27}){}$\1\5
2495 
${}\\{rf}\K\\{incr}(\\{rf},\39(\\{zf}.\|l\GG(\T{2}*(\|k-\T{27})))\AND\T{3});{}$%
2496 
\2\6
2497 
\&{if} ${}((\\{rf}.\|l>\\{xf}.\|l\W\\{rf}.\|h\G\\{xf}.\|h)\V\\{rf}.\|h>\\{xf}.%
2498 
\|h){}$\5
2499 
${}\{{}$\1\6
2500 
${}\\{xf}.\|l\PP{}$;\5
2501 
${}\\{rf}\K\\{ominus}(\\{rf},\39\\{xf}){}$;\5
2502 
${}\\{xf}.\|l\PP;{}$\6
2503 
\4${}\}{}$\2\6
2504 
\4${}\}{}$\2\6
2505 
\&{if} ${}(\\{rf}.\|h\V\\{rf}.\|l){}$\1\5
2506 
${}\\{xf}.\|l\PP{}$;\C{ sticky bit }\2\6
2507 
\&{return} \\{fpack}${}(\\{xf},\39\\{xe},\39\.{'+'},\39\|r){}$;\par
2508 
\U91.\fi
2509 
 
2510 
\M{93}And finally, the genuine floating point remainder. Subroutine \PB{%
2511 
\\{fremstep}}
2512 
either calculates $y\,{\rm rem}\,z$ or reduces $y$ to a smaller number
2513 
having the same remainder with respect to~$z$. In the latter case
2514 
the \PB{\.{E\_BIT}} is set in \PB{\\{exceptions}}. A third parameter, \PB{%
2515 
\\{delta}},
2516 
gives a decrease in exponent that is acceptable for incomplete results;
2517 
if \PB{\\{delta}} is sufficiently large, say 2500, the correct result will
2518 
always be obtained in one step of \PB{\\{fremstep}}.
2519 
 
2520 
\Y\B\4\X5:Subroutines\X${}\mathrel+\E{}$\6
2521 
\&{octa} \\{fremstep}\,\,${}\.{ARGS}((\&{octa},\39\&{octa},\39\&{int})){}$;\5
2522 
\hbox{}\6{}\&{octa} ${}\\{fremstep}(\|y,\39\|z,\39\\{delta}){}$\1\1\6
2523 
\&{octa} \|y${},{}$ \|z;\6
2524 
\&{int} \\{delta};\2\2\6
2525 
${}\{{}$\1\6
2526 
\&{ftype} \\{yt}${},{}$ \\{zt};\6
2527 
\&{int} \\{ye}${},{}$ \\{ze};\6
2528 
\&{char} \\{xs}${},{}$ \\{ys}${},{}$ \\{zs};\6
2529 
\&{octa} \|x${},{}$ \\{xf}${},{}$ \\{yf}${},{}$ \\{zf};\6
2530 
\&{register} \&{int} \\{xe}${},{}$ \\{thresh}${},{}$ \\{odd};\7
2531 
${}\\{yt}\K\\{funpack}(\|y,\39{\AND}\\{yf},\39{\AND}\\{ye},\39{\AND}\\{ys});{}$%
2532 
\6
2533 
${}\\{zt}\K\\{funpack}(\|z,\39{\AND}\\{zf},\39{\AND}\\{ze},\39{\AND}\\{zs});{}$%
2534 
\6
2535 
\&{switch} ${}(\T{4}*\\{yt}+\\{zt}){}$\5
2536 
${}\{{}$\1\6
2537 
\hbox{\4}\X42:The usual NaN cases\X;\6
2538 
\4\&{case} \T{4}${}*\\{zro}+\\{zro}{}$:\5
2539 
\&{case} \T{4}${}*\\{num}+\\{zro}{}$:\5
2540 
\&{case} \T{4}${}*\\{inf}+\\{zro}{}$:\5
2541 
\&{case} \T{4}${}*\\{inf}+\\{num}{}$:\5
2542 
\&{case} \T{4}${}*\\{inf}+\\{inf}{}$:\5
2543 
${}\|x\K\\{standard\_NaN};{}$\6
2544 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{I\_BIT}{}$;\5
2545 
\&{break};\6
2546 
\4\&{case} \T{4}${}*\\{zro}+\\{num}{}$:\5
2547 
\&{case} \T{4}${}*\\{zro}+\\{inf}{}$:\5
2548 
\&{case} \T{4}${}*\\{num}+\\{inf}{}$:\5
2549 
\&{return} \|y;\6
2550 
\4\&{case} \T{4}${}*\\{num}+\\{num}{}$:\5
2551 
\X94:Remainderize nonzero numbers and \PB{\&{return}}\X;\6
2552 
\4\\{zero\_out}:\5
2553 
${}\|x\K\\{zero\_octa};{}$\6
2554 
\4${}\}{}$\2\6
2555 
\&{if} ${}(\\{ys}\E\.{'-'}){}$\1\5
2556 
${}\|x.\|h\MRL{{\OR}{\K}}\\{sign\_bit};{}$\2\6
2557 
\&{return} \|x;\6
2558 
\4${}\}{}$\2\par
2559 
\fi
2560 
 
2561 
\M{94}If there's a huge difference in exponents and the remainder is nonzero,
2562 
this computation will take a long time. One could compute
2563 
$(2^ny)\,{\rm rem}\,z$ much more quickly for large~$n$ by using $O(\log n)$
2564 
multiplications modulo~$z$, but the floating remainder operation isn't
2565 
important enough to justify such expensive hardware.
2566 
 
2567 
Results of floating remainder are always exact, so the rounding mode
2568 
is immaterial.
2569 
 
2570 
\Y\B\4\X94:Remainderize nonzero numbers and \PB{\&{return}}\X${}\E{}$\6
2571 
$\\{odd}\K\T{0}{}$;\C{ will be 1 if we've subtracted an odd multiple of~$z$
2572 
from $y$ }\6
2573 
${}\\{thresh}\K\\{ye}-\\{delta};{}$\6
2574 
\&{if} ${}(\\{thresh}<\\{ze}){}$\1\5
2575 
${}\\{thresh}\K\\{ze};{}$\2\6
2576 
\&{while} ${}(\\{ye}\G\\{thresh}){}$\1\5
2577 
\X95:Reduce \PB{$(\\{ye},\\{yf})$} by a multiple of \PB{\\{zf}}; \PB{\&{goto} %
2578 
\\{zero\_out}} if the remainder is zero, \PB{\&{goto} \\{try\_complement}} if
2579 
appropriate\X;\2\6
2580 
\&{if} ${}(\\{ye}\G\\{ze}){}$\5
2581 
${}\{{}$\1\6
2582 
${}\\{exceptions}\MRL{{\OR}{\K}}\.{E\_BIT}{}$;\5
2583 
\&{return} \\{fpack}${}(\\{yf},\39\\{ye},\39\\{ys},\39\.{ROUND\_OFF});{}$\6
2584 
\4${}\}{}$\2\6
2585 
\&{if} ${}(\\{ye}<\\{ze}-\T{1}){}$\1\5
2586 
\&{return} \\{fpack}${}(\\{yf},\39\\{ye},\39\\{ys},\39\.{ROUND\_OFF});{}$\2\6
2587 
${}\\{yf}\K\\{shift\_right}(\\{yf},\39\T{1},\39\T{1});{}$\6
2588 
\4\\{try\_complement}:\5
2589 
${}\\{xf}\K\\{ominus}(\\{zf},\39\\{yf}),\39\\{xe}\K\\{ze},\39\\{xs}\K\.{'+'}+%
2590 
\.{'-'}-\\{ys};{}$\6
2591 
\&{if} ${}(\\{xf}.\|h>\\{yf}.\|h\V(\\{xf}.\|h\E\\{yf}.\|h\W(\\{xf}.\|l>\\{yf}.%
2592 
\|l\V(\\{xf}.\|l\E\\{yf}.\|l\W\R\\{odd})))){}$\1\5
2593 
${}\\{xf}\K\\{yf},\39\\{xs}\K\\{ys};{}$\2\6
2594 
\&{while} ${}(\\{xf}.\|h<\T{\^400000}){}$\1\5
2595 
${}\\{xe}\MM,\39\\{xf}\K\\{shift\_left}(\\{xf},\39\T{1});{}$\2\6
2596 
\&{return} \\{fpack}${}(\\{xf},\39\\{xe},\39\\{xs},\39\.{ROUND\_OFF}){}$;\par
2597 
\U93.\fi
2598 
 
2599 
\M{95}Here we are careful not to change the sign of \PB{\|y}, because a
2600 
remainder
2601 
of~0 is supposed to inherit the original sign of~\PB{\|y}.
2602 
 
2603 
\Y\B\4\X95:Reduce \PB{$(\\{ye},\\{yf})$} by a multiple of \PB{\\{zf}}; \PB{%
2604 
\&{goto} \\{zero\_out}} if the remainder is zero, \PB{\&{goto} \\{try%
2605 
\_complement}} if appropriate\X${}\E{}$\6
2606 
${}\{{}$\1\6
2607 
\&{if} ${}(\\{yf}.\|h\E\\{zf}.\|h\W\\{yf}.\|l\E\\{zf}.\|l){}$\1\5
2608 
\&{goto} \\{zero\_out};\2\6
2609 
\&{if} ${}(\\{yf}.\|h<\\{zf}.\|h\V(\\{yf}.\|h\E\\{zf}.\|h\W\\{yf}.\|l<\\{zf}.%
2610 
\|l)){}$\5
2611 
${}\{{}$\1\6
2612 
\&{if} ${}(\\{ye}\E\\{ze}){}$\1\5
2613 
\&{goto} \\{try\_complement};\2\6
2614 
${}\\{ye}\MM,\39\\{yf}\K\\{shift\_left}(\\{yf},\39\T{1});{}$\6
2615 
\4${}\}{}$\2\6
2616 
${}\\{yf}\K\\{ominus}(\\{yf},\39\\{zf});{}$\6
2617 
\&{if} ${}(\\{ye}\E\\{ze}){}$\1\5
2618 
${}\\{odd}\K\T{1};{}$\2\6
2619 
\&{while} ${}(\\{yf}.\|h<\T{\^400000}){}$\1\5
2620 
${}\\{ye}\MM,\39\\{yf}\K\\{shift\_left}(\\{yf},\39\T{1});{}$\2\6
2621 
\4${}\}{}$\2\par
2622 
\U94.\fi
2623 
 
2624 
\N{1}{96}Index.
2625 
 
2626 
\fi
2627 
 
2628 
 
2629 
\inx
2630 
\fin
2631 
\con

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.