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
|