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

Subversion Repositories eco32

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

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.