projects / simh.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
First Commit of my working state
[simh.git] / VAX / vax_fpa.c
1 /* vax_fpa.c - VAX f_, d_, g_floating instructions
2
3 Copyright (c) 1998-2008, Robert M Supnik
4
5 Permission is hereby granted, free of charge, to any person obtaining a
6 copy of this software and associated documentation files (the "Software"),
7 to deal in the Software without restriction, including without limitation
8 the rights to use, copy, modify, merge, publish, distribute, sublicense,
9 and/or sell copies of the Software, and to permit persons to whom the
10 Software is furnished to do so, subject to the following conditions:
11
12 The above copyright notice and this permission notice shall be included in
13 all copies or substantial portions of the Software.
14
15 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
18 ROBERT M SUPNIK BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
19 IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
20 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21
22 Except as contained in this notice, the name of Robert M Supnik shall not be
23 used in advertising or otherwise to promote the sale, use or other dealings
24 in this Software without prior written authorization from Robert M Supnik.
25
26 28-May-08 RMS Inlined physical memory routines
27 16-May-06 RMS Fixed bug in 32b floating multiply routine
28 Fixed bug in 64b extended modulus routine
29 03-May-06 RMS Fixed POLYD, POLYG to clear R4, R5
30 Fixed POLYD, POLYG to set R3 correctly
31 Fixed POLYD, POLYG to not exit prematurely if arg = 0
32 Fixed POLYD, POLYG to do full 64b multiply
33 Fixed POLYF, POLYD, POLYG to remove truncation on add
34 Fixed POLYF, POLYD, POLYG to mask mul reslt to 31b/63b/63b
35 Fixed fp add routine to test for zero via fraction
36 to support "denormal" argument from POLYF, POLYD, POLYG
37 (all reported by Tim Stark)
38 27-Sep-05 RMS Fixed bug in 32b structure definitions (from Jason Stevens)
39 30-Sep-04 RMS Comment and formating changes based on vax_octa.c
40 18-Apr-04 RMS Moved format definitions to vax_defs.h
41 19-Jun-03 RMS Simplified add algorithm
42 16-May-03 RMS Fixed bug in floating to integer convert overflow
43 Fixed multiple bugs in EMODx
44 Integrated 32b only code
45 05-Jul-02 RMS Changed internal routine names for C library conflict
46 17-Apr-02 RMS Fixed bug in EDIV zero quotient
47
48 This module contains the instruction simulators for
49
50 - 64 bit arithmetic (ASHQ, EMUL, EDIV)
51 - single precision floating point
52 - double precision floating point, D and G format
53 */
54
55 #include "vax_defs.h"
56 #include <setjmp.h>
57
58 extern int32 R[16];
59 extern int32 PSL;
60 extern int32 p1;
61 extern jmp_buf save_env;
62
63 #if defined (USE_INT64)
64
65 #define M64 0xFFFFFFFFFFFFFFFF /* 64b */
66 #define FD_FRACW (0xFFFF & ~(FD_EXP | FPSIGN))
67 #define FD_FRACL (FD_FRACW | 0xFFFF0000) /* f/d fraction */
68 #define G_FRACW (0xFFFF & ~(G_EXP | FPSIGN))
69 #define G_FRACL (G_FRACW | 0xFFFF0000) /* g fraction */
70 #define UNSCRAM(h,l) (((((t_uint64) (h)) << 48) & 0xFFFF000000000000) | \
71 ((((t_uint64) (h)) << 16) & 0x0000FFFF00000000) | \
72 ((((t_uint64) (l)) << 16) & 0x00000000FFFF0000) | \
73 ((((t_uint64) (l)) >> 16) & 0x000000000000FFFF))
74 #define CONCAT(h,l) ((((t_uint64) (h)) << 32) | ((uint32) (l)))
75
76 typedef struct {
77 int32 sign;
78 int32 exp;
79 t_uint64 frac;
80 } UFP;
81
82 #define UF_NM 0x8000000000000000 /* normalized */
83 #define UF_FRND 0x0000008000000000 /* F round */
84 #define UF_DRND 0x0000000000000080 /* D round */
85 #define UF_GRND 0x0000000000000400 /* G round */
86 #define UF_V_NM 63
87 #define UF_V_FDHI 40
88 #define UF_V_FDLO (UF_V_FDHI - 32)
89 #define UF_V_GHI 43
90 #define UF_V_GLO (UF_V_GHI - 32)
91 #define UF_GETFDHI(x) (int32) ((((x) >> (16 + UF_V_FDHI)) & FD_FRACW) | \
92 (((x) >> (UF_V_FDHI - 16)) & ~0xFFFF))
93 #define UF_GETFDLO(x) (int32) ((((x) >> (16 + UF_V_FDLO)) & 0xFFFF) | \
94 (((x) << (16 - UF_V_FDLO)) & ~0xFFFF))
95 #define UF_GETGHI(x) (int32) ((((x) >> (16 + UF_V_GHI)) & G_FRACW) | \
96 (((x) >> (UF_V_GHI - 16)) & ~0xFFFF))
97 #define UF_GETGLO(x) (int32) ((((x) >> (16 + UF_V_GLO)) & 0xFFFF) | \
98 (((x) << (16 - UF_V_GLO)) & ~0xFFFF))
99
100 void unpackf (int32 hi, UFP *a);
101 void unpackd (int32 hi, int32 lo, UFP *a);
102 void unpackg (int32 hi, int32 lo, UFP *a);
103 void norm (UFP *a);
104 int32 rpackfd (UFP *a, int32 *rh);
105 int32 rpackg (UFP *a, int32 *rh);
106 void vax_fadd (UFP *a, UFP *b);
107 void vax_fmul (UFP *a, UFP *b, t_bool qd, int32 bias, uint32 mhi, uint32 mlo);
108 void vax_fdiv (UFP *b, UFP *a, int32 prec, int32 bias);
109 void vax_fmod (UFP *a, int32 bias, int32 *intgr, int32 *flg);
110
111 /* Quadword arithmetic shift
112
113 opnd[0] = shift count (cnt.rb)
114 opnd[1:2] = source (src.rq)
115 opnd[3:4] = destination (dst.wq)
116 */
117
118 int32 op_ashq (int32 *opnd, int32 *rh, int32 *flg)
119 {
120 t_int64 src, r;
121 int32 sc = opnd[0];
122
123 src = CONCAT (opnd[2], opnd[1]); /* build src */
124 if (sc & BSIGN) { /* right shift? */
125 *flg = 0; /* no ovflo */
126 sc = 0x100 - sc; /* |shift| */
127 if (sc > 63) r = (opnd[2] & LSIGN)? -1: 0; /* sc > 63? */
128 else r = src >> sc;
129 }
130 else {
131 if (sc > 63) { /* left shift */
132 r = 0; /* sc > 63? */
133 *flg = (src != 0); /* ovflo test */
134 }
135 else {
136 r = src << sc; /* do shift */
137 *flg = (src != (r >> sc)); /* ovflo test */
138 }
139 }
140 *rh = (int32) ((r >> 32) & LMASK); /* hi result */
141 return ((int32) (r & LMASK)); /* lo result */
142 }
143
144 /* Extended multiply subroutine */
145
146 int32 op_emul (int32 mpy, int32 mpc, int32 *rh)
147 {
148 t_int64 lmpy = mpy;
149 t_int64 lmpc = mpc;
150
151 lmpy = lmpy * lmpc;
152 *rh = (int32) ((lmpy >> 32) & LMASK);
153 return ((int32) (lmpy & LMASK));
154 }
155
156 /* Extended divide
157
158 opnd[0] = divisor (non-zero)
159 opnd[1:2] = dividend
160 */
161
162 int32 op_ediv (int32 *opnd, int32 *rh, int32 *flg)
163 {
164 t_int64 ldvd, ldvr;
165 int32 quo, rem;
166
167 *flg = CC_V; /* assume error */
168 *rh = 0;
169 ldvr = ((opnd[0] & LSIGN)? -opnd[0]: opnd[0]) & LMASK; /* |divisor| */
170 ldvd = CONCAT (opnd[2], opnd[1]); /* 64b dividend */
171 if (opnd[2] & LSIGN) ldvd = -ldvd; /* |dividend| */
172 if (((ldvd >> 32) & LMASK) >= ldvr) return opnd[1]; /* divide work? */
173 quo = (int32) (ldvd / ldvr); /* do divide */
174 rem = (int32) (ldvd % ldvr);
175 if ((opnd[0] ^ opnd[2]) & LSIGN) { /* result -? */
176 quo = -quo; /* negate */
177 if (quo && ((quo & LSIGN) == 0)) return opnd[1]; /* right sign? */
178 }
179 else if (quo & LSIGN) return opnd[1];
180 if (opnd[2] & LSIGN) rem = -rem; /* sign of rem */
181 *flg = 0; /* no overflow */
182 *rh = rem & LMASK; /* set rem */
183 return (quo & LMASK); /* return quo */
184 }
185
186 /* Compare floating */
187
188 int32 op_cmpfd (int32 h1, int32 l1, int32 h2, int32 l2)
189 {
190 t_uint64 n1, n2;
191
192 if ((h1 & FD_EXP) == 0) {
193 if (h1 & FPSIGN) RSVD_OPND_FAULT;
194 h1 = l1 = 0;
195 }
196 if ((h2 & FD_EXP) == 0) {
197 if (h2 & FPSIGN) RSVD_OPND_FAULT;
198 h2 = l2 = 0;
199 }
200 if ((h1 ^ h2) & FPSIGN) return ((h1 & FPSIGN)? CC_N: 0);
201 n1 = UNSCRAM (h1, l1);
202 n2 = UNSCRAM (h2, l2);
203 if (n1 == n2) return CC_Z;
204 return (((n1 < n2) ^ ((h1 & FPSIGN) != 0))? CC_N: 0);
205 }
206
207 int32 op_cmpg (int32 h1, int32 l1, int32 h2, int32 l2)
208 {
209 t_uint64 n1, n2;
210
211 if ((h1 & G_EXP) == 0) {
212 if (h1 & FPSIGN) RSVD_OPND_FAULT;
213 h1 = l1 = 0;
214 }
215 if ((h2 & G_EXP) == 0) {
216 if (h2 & FPSIGN) RSVD_OPND_FAULT;
217 h2 = l2 = 0;
218 }
219 if ((h1 ^ h2) & FPSIGN) return ((h1 & FPSIGN)? CC_N: 0);
220 n1 = UNSCRAM (h1, l1);
221 n2 = UNSCRAM (h2, l2);
222 if (n1 == n2) return CC_Z;
223 return (((n1 < n2) ^ ((h1 & FPSIGN) != 0))? CC_N: 0);
224 }
225
226 /* Integer to floating convert */
227
228 int32 op_cvtifdg (int32 val, int32 *rh, int32 opc)
229 {
230 UFP a;
231
232 if (val == 0) {
233 if (rh) *rh = 0;
234 return 0;
235 }
236 if (val < 0) {
237 a.sign = FPSIGN;
238 val = - val;
239 }
240 else a.sign = 0;
241 a.exp = 32 + ((opc & 0x100)? G_BIAS: FD_BIAS);
242 a.frac = ((t_uint64) val) << (UF_V_NM - 31);
243 norm (&a);
244 if (opc & 0x100) return rpackg (&a, rh);
245 return rpackfd (&a, rh);
246 }
247
248 /* Floating to integer convert */
249
250 int32 op_cvtfdgi (int32 *opnd, int32 *flg, int32 opc)
251 {
252 UFP a;
253 int32 lnt = opc & 03;
254 int32 ubexp;
255 static t_uint64 maxv[4] = { 0x7F, 0x7FFF, 0x7FFFFFFF, 0x7FFFFFFF };
256
257 *flg = 0;
258 if (opc & 0x100) {
259 unpackg (opnd[0], opnd[1], &a);
260 ubexp = a.exp - G_BIAS;
261 }
262 else {
263 if (opc & 0x20) unpackd (opnd[0], opnd[1], &a);
264 else unpackf (opnd[0], &a);
265 ubexp = a.exp - FD_BIAS;
266 }
267 if ((a.exp == 0) || (ubexp < 0)) return 0;
268 if (ubexp <= UF_V_NM) {
269 a.frac = a.frac >> (UF_V_NM - ubexp); /* leave rnd bit */
270 if ((opc & 03) == 03) a.frac = a.frac + 1; /* if CVTR, round */
271 a.frac = a.frac >> 1; /* now justified */
272 if (a.frac > (maxv[lnt] + (a.sign? 1: 0))) *flg = CC_V;
273 }
274 else {
275 *flg = CC_V; /* set overflow */
276 if (ubexp > (UF_V_NM + 32)) return 0;
277 a.frac = a.frac << (ubexp - UF_V_NM - 1); /* no rnd bit */
278 }
279 return ((int32) ((a.sign? (a.frac ^ LMASK) + 1: a.frac) & LMASK));
280 }
281
282 /* Extended modularize
283
284 One of three floating point instructions dropped from the architecture,
285 EMOD presents two sets of complications. First, it requires an extended
286 fraction multiply, with precise (and unusual) truncation conditions.
287 Second, it has two write operands, a dubious distinction it shares
288 with EDIV.
289 */
290
291 int32 op_emodf (int32 *opnd, int32 *intgr, int32 *flg)
292 {
293 UFP a, b;
294
295 unpackf (opnd[0], &a); /* unpack operands */
296 unpackf (opnd[2], &b);
297 a.frac = a.frac | (((t_uint64) opnd[1]) << 32); /* extend src1 */
298 vax_fmul (&a, &b, 0, FD_BIAS, 0, LMASK); /* multiply */
299 vax_fmod (&a, FD_BIAS, intgr, flg); /* sep int & frac */
300 return rpackfd (&a, NULL); /* return frac */
301 }
302
303 int32 op_emodd (int32 *opnd, int32 *flo, int32 *intgr, int32 *flg)
304 {
305 UFP a, b;
306
307 unpackd (opnd[0], opnd[1], &a); /* unpack operands */
308 unpackd (opnd[3], opnd[4], &b);
309 a.frac = a.frac | opnd[2]; /* extend src1 */
310 vax_fmul (&a, &b, 1, FD_BIAS, 0, 0); /* multiply */
311 vax_fmod (&a, FD_BIAS, intgr, flg); /* sep int & frac */
312 return rpackfd (&a, flo); /* return frac */
313 }
314
315 int32 op_emodg (int32 *opnd, int32 *flo, int32 *intgr, int32 *flg)
316 {
317 UFP a, b;
318
319 unpackg (opnd[0], opnd[1], &a); /* unpack operands */
320 unpackg (opnd[3], opnd[4], &b);
321 a.frac = a.frac | (opnd[2] >> 5); /* extend src1 */
322 vax_fmul (&a, &b, 1, G_BIAS, 0, 0); /* multiply */
323 vax_fmod (&a, G_BIAS, intgr, flg); /* sep int & frac */
324 return rpackg (&a, flo); /* return frac */
325 }
326
327 /* Unpacked floating point routines */
328
329 void vax_fadd (UFP *a, UFP *b)
330 {
331 int32 ediff;
332 UFP t;
333
334 if (a->frac == 0) { /* s1 = 0? */
335 *a = *b;
336 return;
337 }
338 if (b->frac == 0) return; /* s2 = 0? */
339 if ((a->exp < b->exp) || /* |s1| < |s2|? swap */
340 ((a->exp == b->exp) && (a->frac < b->frac))) {
341 t = *a;
342 *a = *b;
343 *b = t;
344 }
345 ediff = a->exp - b->exp; /* exp diff */
346 if (a->sign ^ b->sign) { /* eff sub? */
347 if (ediff) { /* exp diff? */
348 if (ediff > 63) b->frac = M64; /* retain sticky */
349 else b->frac = ((-((t_int64) b->frac) >> ediff) | /* denormalize */
350 (M64 << (64 - ediff))); /* preserve sign */
351 a->frac = a->frac + b->frac; /* add frac */
352 }
353 else a->frac = a->frac - b->frac; /* sub frac */
354 norm (a); /* normalize */
355 }
356 else {
357 if (ediff > 63) b->frac = 0; /* add */
358 else if (ediff) b->frac = b->frac >> ediff; /* denormalize */
359 a->frac = a->frac + b->frac; /* add frac */
360 if (a->frac < b->frac) { /* chk for carry */
361 a->frac = UF_NM | (a->frac >> 1); /* shift in carry */
362 a->exp = a->exp + 1; /* skip norm */
363 }
364 }
365 return;
366 }
367
368 /* Floating multiply - 64b * 64b with cross products */
369
370 void vax_fmul (UFP *a, UFP *b, t_bool qd, int32 bias, uint32 mhi, uint32 mlo)
371 {
372 t_uint64 ah, bh, al, bl, rhi, rlo, rmid1, rmid2;
373 t_uint64 mask = (((t_uint64) mhi) << 32) | ((t_uint64) mlo);
374
375 if ((a->exp == 0) || (b->exp == 0)) { /* zero argument? */
376 a->frac = a->sign = a->exp = 0; /* result is zero */
377 return;
378 }
379 a->sign = a->sign ^ b->sign; /* sign of result */
380 a->exp = a->exp + b->exp - bias; /* add exponents */
381 ah = (a->frac >> 32) & LMASK; /* split operands */
382 bh = (b->frac >> 32) & LMASK; /* into 32b chunks */
383 rhi = ah * bh; /* high result */
384 if (qd) { /* 64b needed? */
385 al = a->frac & LMASK;
386 bl = b->frac & LMASK;
387 rmid1 = ah * bl;
388 rmid2 = al * bh;
389 rlo = al * bl;
390 rhi = rhi + ((rmid1 >> 32) & LMASK) + ((rmid2 >> 32) & LMASK);
391 rmid1 = rlo + (rmid1 << 32); /* add mid1 to lo */
392 if (rmid1 < rlo) rhi = rhi + 1; /* carry? incr hi */
393 rmid2 = rmid1 + (rmid2 << 32); /* add mid2 to lo */
394 if (rmid2 < rmid1) rhi = rhi + 1; /* carry? incr hi */
395 }
396 a->frac = rhi & ~mask;
397 norm (a); /* normalize */
398 return;
399 }
400
401 /* Floating modulus - there are three cases
402
403 exp <= bias - integer is 0, fraction is input,
404 no overflow
405 bias < exp <= bias+64 - separate integer and fraction,
406 integer overflow may occur
407 bias+64 < exp - result is integer, fraction is 0
408 integer overflow
409 */
410
411 void vax_fmod (UFP *a, int32 bias, int32 *intgr, int32 *flg)
412 {
413 if (a->exp <= bias) *intgr = *flg = 0; /* 0 or <1? int = 0 */
414 else if (a->exp <= (bias + 64)) { /* in range [1,64]? */
415 *intgr = (int32) (a->frac >> (64 - (a->exp - bias)));
416 if ((a->exp > (bias + 32)) || /* test ovflo */
417 ((a->exp == (bias + 32)) &&
418 (((uint32) *intgr) > (a->sign? 0x80000000: 0x7FFFFFFF))))
419 *flg = CC_V;
420 else *flg = 0;
421 if (a->sign) *intgr = -*intgr; /* -? comp int */
422 if (a->exp == (bias + 64)) a->frac = 0; /* special case 64 */
423 else a->frac = a->frac << (a->exp - bias);
424 a->exp = bias;
425 }
426 else {
427 *intgr = 0; /* out of range */
428 a->frac = a->sign = a->exp = 0; /* result 0 */
429 *flg = CC_V; /* overflow */
430 }
431 norm (a); /* normalize */
432 return;
433 }
434
435 /* Floating divide
436 Needs to develop at least one rounding bit. Since the first
437 divide step can fail, caller should specify 2 more bits than
438 the precision of the fraction.
439 */
440
441 void vax_fdiv (UFP *a, UFP *b, int32 prec, int32 bias)
442 {
443 int32 i;
444 t_uint64 quo = 0;
445
446 if (a->exp == 0) FLT_DZRO_FAULT; /* divr = 0? */
447 if (b->exp == 0) return; /* divd = 0? */
448 b->sign = b->sign ^ a->sign; /* result sign */
449 b->exp = b->exp - a->exp + bias + 1; /* unbiased exp */
450 a->frac = a->frac >> 1; /* allow 1 bit left */
451 b->frac = b->frac >> 1;
452 for (i = 0; (i < prec) && b->frac; i++) { /* divide loop */
453 quo = quo << 1; /* shift quo */
454 if (b->frac >= a->frac) { /* div step ok? */
455 b->frac = b->frac - a->frac; /* subtract */
456 quo = quo + 1; /* quo bit = 1 */
457 }
458 b->frac = b->frac << 1; /* shift divd */
459 }
460 b->frac = quo << (UF_V_NM - i + 1); /* shift quo */
461 norm (b); /* normalize */
462 return;
463 }
464
465 /* Support routines */
466
467 void unpackf (int32 hi, UFP *r)
468 {
469 r->sign = hi & FPSIGN; /* get sign */
470 r->exp = FD_GETEXP (hi); /* get exponent */
471 if (r->exp == 0) { /* exp = 0? */
472 if (r->sign) RSVD_OPND_FAULT; /* if -, rsvd op */
473 r->frac = 0; /* else 0 */
474 return;
475 }
476 hi = (((hi & FD_FRACW) | FD_HB) << 16) | ((hi >> 16) & 0xFFFF);
477 r->frac = ((t_uint64) hi) << (32 + UF_V_FDLO);
478 return;
479 }
480
481 void unpackd (int32 hi, int32 lo, UFP *r)
482 {
483 r->sign = hi & FPSIGN; /* get sign */
484 r->exp = FD_GETEXP (hi); /* get exponent */
485 if (r->exp == 0) { /* exp = 0? */
486 if (r->sign) RSVD_OPND_FAULT; /* if -, rsvd op */
487 r->frac = 0; /* else 0 */
488 return;
489 }
490 hi = (hi & FD_FRACL) | FD_HB; /* canonical form */
491 r->frac = UNSCRAM (hi, lo) << UF_V_FDLO; /* guard bits */
492 return;
493 }
494
495 void unpackg (int32 hi, int32 lo, UFP *r)
496 {
497 r->sign = hi & FPSIGN; /* get sign */
498 r->exp = G_GETEXP (hi); /* get exponent */
499 if (r->exp == 0) { /* exp = 0? */
500 if (r->sign) RSVD_OPND_FAULT; /* if -, rsvd op */
501 r->frac = 0; /* else 0 */
502 return;
503 }
504 hi = (hi & G_FRACL) | G_HB; /* canonical form */
505 r->frac = UNSCRAM (hi, lo) << UF_V_GLO; /* guard bits */
506 return;
507 }
508
509 void norm (UFP *r)
510 {
511 int32 i;
512 static t_uint64 normmask[5] = {
513 0xc000000000000000, 0xf000000000000000, 0xff00000000000000,
514 0xffff000000000000, 0xffffffff00000000
515 };
516 static int32 normtab[6] = { 1, 2, 4, 8, 16, 32};
517
518 if (r->frac == 0) { /* if fraction = 0 */
519 r->sign = r->exp = 0; /* result is 0 */
520 return;
521 }
522 while ((r->frac & UF_NM) == 0) { /* normalized? */
523 for (i = 0; i < 5; i++) { /* find first 1 */
524 if (r->frac & normmask[i]) break;
525 }
526 r->frac = r->frac << normtab[i]; /* shift frac */
527 r->exp = r->exp - normtab[i]; /* decr exp */
528 }
529 return;
530 }
531
532 int32 rpackfd (UFP *r, int32 *rh)
533 {
534 if (rh) *rh = 0; /* assume 0 */
535 if (r->frac == 0) return 0; /* result 0? */
536 r->frac = r->frac + (rh? UF_DRND: UF_FRND); /* round */
537 if ((r->frac & UF_NM) == 0) { /* carry out? */
538 r->frac = r->frac >> 1; /* renormalize */
539 r->exp = r->exp + 1;
540 }
541 if (r->exp > (int32) FD_M_EXP) FLT_OVFL_FAULT; /* ovflo? fault */
542 if (r->exp <= 0) { /* underflow? */
543 if (PSL & PSW_FU) FLT_UNFL_FAULT; /* fault if fu */
544 return 0; /* else 0 */
545 }
546 if (rh) *rh = UF_GETFDLO (r->frac); /* get low */
547 return r->sign | (r->exp << FD_V_EXP) | UF_GETFDHI (r->frac);
548 }
549
550 int32 rpackg (UFP *r, int32 *rh)
551 {
552 *rh = 0; /* assume 0 */
553 if (r->frac == 0) return 0; /* result 0? */
554 r->frac = r->frac + UF_GRND; /* round */
555 if ((r->frac & UF_NM) == 0) { /* carry out? */
556 r->frac = r->frac >> 1; /* renormalize */
557 r->exp = r->exp + 1;
558 }
559 if (r->exp > (int32) G_M_EXP) FLT_OVFL_FAULT; /* ovflo? fault */
560 if (r->exp <= 0) { /* underflow? */
561 if (PSL & PSW_FU) FLT_UNFL_FAULT; /* fault if fu */
562 return 0; /* else 0 */
563 }
564 if (rh) *rh = UF_GETGLO (r->frac); /* get low */
565 return r->sign | (r->exp << G_V_EXP) | UF_GETGHI (r->frac);
566 }
567
568 #else /* 32b code */
569
570 #define WORDSWAP(x) ((((x) & WMASK) << 16) | (((x) >> 16) & WMASK))
571
572 typedef struct {
573 uint32 lo;
574 uint32 hi;
575 } UDP;
576
577 typedef struct {
578 int32 sign;
579 int32 exp;
580 UDP frac;
581 } UFP;
582
583 #define UF_NM_H 0x80000000 /* normalized */
584 #define UF_FRND_H 0x00000080 /* F round */
585 #define UF_FRND_L 0x00000000
586 #define UF_DRND_H 0x00000000 /* D round */
587 #define UF_DRND_L 0x00000080
588 #define UF_GRND_H 0x00000000 /* G round */
589 #define UF_GRND_L 0x00000400
590 #define UF_V_NM 63
591
592 void unpackf (uint32 hi, UFP *a);
593 void unpackd (uint32 hi, uint32 lo, UFP *a);
594 void unpackg (uint32 hi, uint32 lo, UFP *a);
595 void norm (UFP *a);
596 int32 rpackfd (UFP *a, int32 *rh);
597 int32 rpackg (UFP *a, int32 *rh);
598 void vax_fadd (UFP *a, UFP *b);
599 void vax_fmul (UFP *a, UFP *b, t_bool qd, int32 bias, uint32 mhi, uint32 mlo);
600 void vax_fmod (UFP *a, int32 bias, int32 *intgr, int32 *flg);
601 void vax_fdiv (UFP *b, UFP *a, int32 prec, int32 bias);
602 void dp_add (UDP *a, UDP *b);
603 void dp_inc (UDP *a);
604 void dp_sub (UDP *a, UDP *b);
605 void dp_imul (uint32 a, uint32 b, UDP *r);
606 void dp_lsh (UDP *a, uint32 sc);
607 void dp_rsh (UDP *a, uint32 sc);
608 void dp_rsh_s (UDP *a, uint32 sc, uint32 neg);
609 void dp_neg (UDP *a);
610 int32 dp_cmp (UDP *a, UDP *b);
611
612 /* Quadword arithmetic shift
613
614 opnd[0] = shift count (cnt.rb)
615 opnd[1:2] = source (src.rq)
616 opnd[3:4] = destination (dst.wq)
617 */
618
619 int32 op_ashq (int32 *opnd, int32 *rh, int32 *flg)
620 {
621 UDP r, sr;
622 uint32 sc = opnd[0];
623
624 r.lo = opnd[1]; /* get source */
625 r.hi = opnd[2];
626 *flg = 0; /* assume no ovflo */
627 if (sc & BSIGN) /* right shift? */
628 dp_rsh_s (&r, 0x100 - sc, r.hi & LSIGN); /* signed right */
629 else {
630 dp_lsh (&r, sc); /* left shift */
631 sr = r; /* copy result */
632 dp_rsh_s (&sr, sc, sr.hi & LSIGN); /* signed right */
633 if ((sr.hi != ((uint32) opnd[2])) || /* reshift != orig? */
634 (sr.lo != ((uint32) opnd[1]))) *flg = 1; /* overflow */
635 }
636 *rh = r.hi; /* hi result */
637 return r.lo; /* lo result */
638 }
639
640 /* Extended multiply subroutine */
641
642 int32 op_emul (int32 mpy, int32 mpc, int32 *rh)
643 {
644 UDP r;
645 int32 sign = mpy ^ mpc; /* sign of result */
646
647 if (mpy & LSIGN) mpy = -mpy; /* abs value */
648 if (mpc & LSIGN) mpc = -mpc;
649 dp_imul (mpy & LMASK, mpc & LMASK, &r); /* 32b * 32b -> 64b */
650 if (sign & LSIGN) dp_neg (&r); /* negative result? */
651 *rh = r.hi;
652 return r.lo;
653 }
654
655 /* Extended divide
656
657 opnd[0] = divisor (non-zero)
658 opnd[1:2] = dividend
659 */
660
661 int32 op_ediv (int32 *opnd, int32 *rh, int32 *flg)
662 {
663 UDP dvd;
664 uint32 i, dvr, quo;
665
666 dvr = opnd[0]; /* get divisor */
667 dvd.lo = opnd[1]; /* get dividend */
668 dvd.hi = opnd[2];
669 *flg = CC_V; /* assume error */
670 *rh = 0;
671 if (dvd.hi & LSIGN) dp_neg (&dvd); /* |dividend| */
672 if (dvr & LSIGN) dvr = NEG (dvr); /* |divisor| */
673 if (dvd.hi >= dvr) return opnd[1]; /* divide work? */
674 for (i = quo = 0; i < 32; i++) { /* 32 iterations */
675 quo = quo << 1; /* shift quotient */
676 dp_lsh (&dvd, 1); /* shift dividend */
677 if (dvd.hi >= dvr) { /* step work? */
678 dvd.hi = (dvd.hi - dvr) & LMASK; /* subtract dvr */
679 quo = quo + 1;
680 }
681 }
682 if ((opnd[0] ^ opnd[2]) & LSIGN) { /* result -? */
683 quo = NEG (quo); /* negate */
684 if (quo && ((quo & LSIGN) == 0)) return opnd[1]; /* right sign? */
685 }
686 else if (quo & LSIGN) return opnd[1];
687 if (opnd[2] & LSIGN) *rh = NEG (dvd.hi); /* sign of rem */
688 else *rh = dvd.hi;
689 *flg = 0; /* no overflow */
690 return quo; /* return quo */
691 }
692
693 /* Compare floating */
694
695 int32 op_cmpfd (int32 h1, int32 l1, int32 h2, int32 l2)
696 {
697 UFP a, b;
698 int32 r;
699
700 unpackd (h1, l1, &a);
701 unpackd (h2, l2, &b);
702 if (a.sign != b.sign) return (a.sign? CC_N: 0);
703 r = a.exp - b.exp;
704 if (r == 0) r = dp_cmp (&a.frac, &b.frac);
705 if (r < 0) return (a.sign? 0: CC_N);
706 if (r > 0) return (a.sign? CC_N: 0);
707 return CC_Z;
708 }
709
710 int32 op_cmpg (int32 h1, int32 l1, int32 h2, int32 l2)
711 {
712 UFP a, b;
713 int32 r;
714
715 unpackg (h1, l1, &a);
716 unpackg (h2, l2, &b);
717 if (a.sign != b.sign) return (a.sign? CC_N: 0);
718 r = a.exp - b.exp;
719 if (r == 0) r = dp_cmp (&a.frac, &b.frac);
720 if (r < 0) return (a.sign? 0: CC_N);
721 if (r > 0) return (a.sign? CC_N: 0);
722 return CC_Z;
723 }
724
725 /* Integer to floating convert */
726
727 int32 op_cvtifdg (int32 val, int32 *rh, int32 opc)
728 {
729 UFP a;
730
731 if (val == 0) { /* zero? */
732 if (rh) *rh = 0; /* return true 0 */
733 return 0;
734 }
735 if (val < 0) { /* negative? */
736 a.sign = FPSIGN; /* sign = - */
737 val = -val;
738 }
739 else a.sign = 0; /* else sign = + */
740 a.exp = 32 + ((opc & 0x100)? G_BIAS: FD_BIAS); /* initial exp */
741 a.frac.hi = val & LMASK; /* fraction */
742 a.frac.lo = 0;
743 norm (&a); /* normalize */
744 if (opc & 0x100) return rpackg (&a, rh); /* pack and return */
745 return rpackfd (&a, rh);
746 }
747
748 /* Floating to integer convert */
749
750 int32 op_cvtfdgi (int32 *opnd, int32 *flg, int32 opc)
751 {
752 UFP a;
753 int32 lnt = opc & 03;
754 int32 ubexp;
755 static uint32 maxv[4] = { 0x7F, 0x7FFF, 0x7FFFFFFF, 0x7FFFFFFF };
756
757 *flg = 0;
758 if (opc & 0x100) { /* G? */
759 unpackg (opnd[0], opnd[1], &a); /* unpack */
760 ubexp = a.exp - G_BIAS; /* unbiased exp */
761 }
762 else {
763 if (opc & 0x20) unpackd (opnd[0], opnd[1], &a); /* F or D */
764 else unpackf (opnd[0], &a); /* unpack */
765 ubexp = a.exp - FD_BIAS; /* unbiased exp */
766 }
767 if ((a.exp == 0) || (ubexp < 0)) return 0; /* true zero or frac? */
768 if (ubexp <= UF_V_NM) { /* exp in range? */
769 dp_rsh (&a.frac, UF_V_NM - ubexp); /* leave rnd bit */
770 if (lnt == 03) dp_inc (&a.frac); /* if CVTR, round */
771 dp_rsh (&a.frac, 1); /* now justified */
772 if ((a.frac.hi != 0) ||
773 (a.frac.lo > (maxv[lnt] + (a.sign? 1: 0)))) *flg = CC_V;
774 }
775 else {
776 *flg = CC_V; /* always ovflo */
777 if (ubexp > (UF_V_NM + 32)) return 0; /* in ext range? */
778 dp_lsh (&a.frac, ubexp - UF_V_NM - 1); /* no rnd bit */
779 }
780 return (a.sign? NEG (a.frac.lo): a.frac.lo); /* return lo frac */
781 }
782
783 /* Extended modularize
784
785 One of three floating point instructions dropped from the architecture,
786 EMOD presents two sets of complications. First, it requires an extended
787 fraction multiply, with precise (and unusual) truncation conditions.
788 Second, it has two write operands, a dubious distinction it shares
789 with EDIV.
790 */
791
792 int32 op_emodf (int32 *opnd, int32 *intgr, int32 *flg)
793 {
794 UFP a, b;
795
796 unpackf (opnd[0], &a); /* unpack operands */
797 unpackf (opnd[2], &b);
798 a.frac.hi = a.frac.hi | opnd[1]; /* extend src1 */
799 vax_fmul (&a, &b, 0, FD_BIAS, 0, LMASK); /* multiply */
800 vax_fmod (&a, FD_BIAS, intgr, flg); /* sep int & frac */
801 return rpackfd (&a, NULL); /* return frac */
802 }
803
804 int32 op_emodd (int32 *opnd, int32 *flo, int32 *intgr, int32 *flg)
805 {
806 UFP a, b;
807
808 unpackd (opnd[0], opnd[1], &a); /* unpack operands */
809 unpackd (opnd[3], opnd[4], &b);
810 a.frac.lo = a.frac.lo | opnd[2]; /* extend src1 */
811 vax_fmul (&a, &b, 1, FD_BIAS, 0, 0); /* multiply */
812 vax_fmod (&a, FD_BIAS, intgr, flg); /* sep int & frac */
813 return rpackfd (&a, flo); /* return frac */
814 }
815
816 int32 op_emodg (int32 *opnd, int32 *flo, int32 *intgr, int32 *flg)
817 {
818 UFP a, b;
819
820 unpackg (opnd[0], opnd[1], &a); /* unpack operands */
821 unpackg (opnd[3], opnd[4], &b);
822 a.frac.lo = a.frac.lo | (opnd[2] >> 5); /* extend src1 */
823 vax_fmul (&a, &b, 1, G_BIAS, 0, 0); /* multiply */
824 vax_fmod (&a, G_BIAS, intgr, flg); /* sep int & frac */
825 return rpackg (&a, flo); /* return frac */
826 }
827
828 /* Unpacked floating point routines */
829
830 /* Floating add */
831
832 void vax_fadd (UFP *a, UFP *b)
833 {
834 int32 ediff;
835 UFP t;
836
837 if ((a->frac.hi == 0) && (a->frac.lo == 0)) { /* s1 = 0? */
838 *a = *b;
839 return;
840 }
841 if ((b->frac.hi == 0) && (b->frac.lo == 0)) return; /* s2 = 0? */
842 if ((a->exp < b->exp) || /* |s1| < |s2|? swap */
843 ((a->exp == b->exp) && (dp_cmp (&a->frac, &b->frac) < 0))) {
844 t = *a;
845 *a = *b;
846 *b = t;
847 }
848 ediff = a->exp - b->exp; /* exp diff */
849 if (a->sign ^ b->sign) { /* eff sub? */
850 if (ediff) { /* exp diff? */
851 dp_neg (&b->frac); /* negate fraction */
852 dp_rsh_s (&b->frac, ediff, 1); /* signed right */
853 dp_add (&a->frac, &b->frac); /* "add" frac */
854 }
855 else dp_sub (&a->frac, &b->frac); /* a >= b */
856 norm (a); /* normalize */
857 }
858 else {
859 if (ediff) dp_rsh (&b->frac, ediff); /* add, denormalize */
860 dp_add (&a->frac, &b->frac); /* add frac */
861 if (dp_cmp (&a->frac, &b->frac) < 0) { /* chk for carry */
862 dp_rsh (&a->frac, 1); /* renormalize */
863 a->frac.hi = a->frac.hi | UF_NM_H; /* add norm bit */
864 a->exp = a->exp + 1; /* skip norm */
865 }
866 }
867 return;
868 }
869
870 /* Floating multiply - 64b * 64b with cross products */
871
872 void vax_fmul (UFP *a, UFP *b, t_bool qd, int32 bias, uint32 mhi, uint32 mlo)
873 {
874 UDP rhi, rlo, rmid1, rmid2;
875
876 if ((a->exp == 0) || (b->exp == 0)) { /* zero argument? */
877 a->frac.hi = a->frac.lo = 0; /* result is zero */
878 a->sign = a->exp = 0;
879 return;
880 }
881 a->sign = a->sign ^ b->sign; /* sign of result */
882 a->exp = a->exp + b->exp - bias; /* add exponents */
883 dp_imul (a->frac.hi, b->frac.hi, &rhi); /* high result */
884 if (qd) { /* 64b needed? */
885 dp_imul (a->frac.hi, b->frac.lo, &rmid1); /* cross products */
886 dp_imul (a->frac.lo, b->frac.hi, &rmid2);
887 dp_imul (a->frac.lo, b->frac.lo, &rlo); /* low result */
888 rhi.lo = (rhi.lo + rmid1.hi) & LMASK; /* add hi cross */
889 if (rhi.lo < rmid1.hi) /* to low high res */
890 rhi.hi = (rhi.hi + 1) & LMASK;
891 rhi.lo = (rhi.lo + rmid2.hi) & LMASK;
892 if (rhi.lo < rmid2.hi)
893 rhi.hi = (rhi.hi + 1) & LMASK;
894 rlo.hi = (rlo.hi + rmid1.lo) & LMASK; /* add mid1 to low res */
895 if (rlo.hi < rmid1.lo) dp_inc (&rhi); /* carry? incr high res */
896 rlo.hi = (rlo.hi + rmid2.lo) & LMASK; /* add mid2 to low res */
897 if (rlo.hi < rmid2.lo) dp_inc (&rhi); /* carry? incr high res */
898 }
899 a->frac.hi = rhi.hi & ~mhi; /* mask fraction */
900 a->frac.lo = rhi.lo & ~mlo;
901 norm (a); /* normalize */
902 return;
903 }
904
905 /* Floating modulus - there are three cases
906
907 exp <= bias - integer is 0, fraction is input,
908 no overflow
909 bias < exp <= bias+64 - separate integer and fraction,
910 integer overflow may occur
911 bias+64 < exp - result is integer, fraction is 0
912 integer overflow
913 */
914
915 void vax_fmod (UFP *a, int32 bias, int32 *intgr, int32 *flg)
916 {
917 UDP ifr;
918
919 if (a->exp <= bias) *intgr = *flg = 0; /* 0 or <1? int = 0 */
920 else if (a->exp <= (bias + 64)) { /* in range [1,64]? */
921 ifr = a->frac;
922 dp_rsh (&ifr, 64 - (a->exp - bias)); /* separate integer */
923 if ((a->exp > (bias + 32)) || /* test ovflo */
924 ((a->exp == (bias + 32)) &&
925 (ifr.lo > (a->sign? 0x80000000: 0x7FFFFFFF))))
926 *flg = CC_V;
927 else *flg = 0;
928 *intgr = ifr.lo;
929 if (a->sign) *intgr = -*intgr; /* -? comp int */
930 dp_lsh (&a->frac, a->exp - bias); /* excise integer */
931 a->exp = bias;
932 }
933 else {
934 *intgr = 0; /* out of range */
935 a->frac.hi = a->frac.lo = a->sign = a->exp = 0; /* result 0 */
936 *flg = CC_V; /* overflow */
937 }
938 norm (a); /* normalize */
939 return;
940 }
941
942 /* Floating divide
943 Needs to develop at least one rounding bit. Since the first
944 divide step can fail, caller should specify 2 more bits than
945 the precision of the fraction.
946 */
947
948 void vax_fdiv (UFP *a, UFP *b, int32 prec, int32 bias)
949 {
950 int32 i;
951 UDP quo = { 0, 0 };
952
953 if (a->exp == 0) FLT_DZRO_FAULT; /* divr = 0? */
954 if (b->exp == 0) return; /* divd = 0? */
955 b->sign = b->sign ^ a->sign; /* result sign */
956 b->exp = b->exp - a->exp + bias + 1; /* unbiased exp */
957 dp_rsh (&a->frac, 1); /* allow 1 bit left */
958 dp_rsh (&b->frac, 1);
959 for (i = 0; i < prec; i++) { /* divide loop */
960 dp_lsh (&quo, 1); /* shift quo */
961 if (dp_cmp (&b->frac, &a->frac) >= 0) { /* div step ok? */
962 dp_sub (&b->frac, &a->frac); /* subtract */
963 quo.lo = quo.lo + 1; /* quo bit = 1 */
964 }
965 dp_lsh (&b->frac, 1); /* shift divd */
966 }
967 dp_lsh (&quo, UF_V_NM - prec + 1); /* put in position */
968 b->frac = quo;
969 norm (b); /* normalize */
970 return;
971 }
972
973 /* Double precision integer routines */
974
975 int32 dp_cmp (UDP *a, UDP *b)
976 {
977 if (a->hi < b->hi) return -1; /* compare hi */
978 if (a->hi > b->hi) return +1;
979 if (a->lo < b->lo) return -1; /* hi =, compare lo */
980 if (a->lo > b->lo) return +1;
981 return 0; /* hi, lo equal */
982 }
983
984 void dp_add (UDP *a, UDP *b)
985 {
986 a->lo = (a->lo + b->lo) & LMASK; /* add lo */
987 if (a->lo < b->lo) a->hi = a->hi + 1; /* carry? */
988 a->hi = (a->hi + b->hi) & LMASK; /* add hi */
989 return;
990 }
991
992 void dp_inc (UDP *a)
993 {
994 a->lo = (a->lo + 1) & LMASK; /* inc lo */
995 if (a->lo == 0) a->hi = (a->hi + 1) & LMASK; /* carry? inc hi */
996 return;
997 }
998
999 void dp_sub (UDP *a, UDP *b)
1000 {
1001 if (a->lo < b->lo) a->hi = a->hi - 1; /* borrow? decr hi */
1002 a->lo = (a->lo - b->lo) & LMASK; /* sub lo */
1003 a->hi = (a->hi - b->hi) & LMASK; /* sub hi */
1004 return;
1005 }
1006
1007 void dp_lsh (UDP *r, uint32 sc)
1008 {
1009 if (sc > 63) r->hi = r->lo = 0; /* > 63? result 0 */
1010 else if (sc > 31) { /* [32,63]? */
1011 r->hi = (r->lo << (sc - 32)) & LMASK;
1012 r->lo = 0;
1013 }
1014 else if (sc != 0) {
1015 r->hi = ((r->hi << sc) | (r->lo >> (32 - sc))) & LMASK;
1016 r->lo = (r->lo << sc) & LMASK;
1017 }
1018 return;
1019 }
1020
1021 void dp_rsh (UDP *r, uint32 sc)
1022 {
1023 if (sc > 63) r->hi = r->lo = 0; /* > 63? result 0 */
1024 else if (sc > 31) { /* [32,63]? */
1025 r->lo = (r->hi >> (sc - 32)) & LMASK;
1026 r->hi = 0;
1027 }
1028 else if (sc != 0) {
1029 r->lo = ((r->lo >> sc) | (r->hi << (32 - sc))) & LMASK;
1030 r->hi = (r->hi >> sc) & LMASK;
1031 }
1032 return;
1033 }
1034
1035 void dp_rsh_s (UDP *r, uint32 sc, uint32 neg)
1036 {
1037 dp_rsh (r, sc); /* do unsigned right */
1038 if (neg && sc) { /* negative? */
1039 if (sc > 63) r->hi = r->lo = LMASK; /* > 63? result -1 */
1040 else {
1041 UDP ones = { LMASK, LMASK };
1042 dp_lsh (&ones, 64 - sc); /* shift ones */
1043 r->hi = r->hi | ones.hi; /* or into result */
1044 r->lo = r->lo | ones.lo;
1045 }
1046 }
1047 return;
1048 }
1049
1050 void dp_imul (uint32 a, uint32 b, UDP *r)
1051 {
1052 uint32 ah, bh, al, bl, rhi, rlo, rmid1, rmid2;
1053
1054 if ((a == 0) || (b == 0)) { /* zero argument? */
1055 r->hi = r->lo = 0; /* result is zero */
1056 return;
1057 }
1058 ah = (a >> 16) & WMASK; /* split operands */
1059 bh = (b >> 16) & WMASK; /* into 16b chunks */
1060 al = a & WMASK;
1061 bl = b & WMASK;
1062 rhi = ah * bh; /* high result */
1063 rmid1 = ah * bl;
1064 rmid2 = al * bh;
1065 rlo = al * bl;
1066 rhi = rhi + ((rmid1 >> 16) & WMASK) + ((rmid2 >> 16) & WMASK);
1067 rmid1 = (rlo + (rmid1 << 16)) & LMASK; /* add mid1 to lo */
1068 if (rmid1 < rlo) rhi = rhi + 1; /* carry? incr hi */
1069 rmid2 = (rmid1 + (rmid2 << 16)) & LMASK; /* add mid2 to to */
1070 if (rmid2 < rmid1) rhi = rhi + 1; /* carry? incr hi */
1071 r->hi = rhi & LMASK; /* mask result */
1072 r->lo = rmid2;
1073 return;
1074 }
1075
1076 void dp_neg (UDP *r)
1077 {
1078 r->lo = NEG (r->lo);
1079 r->hi = (~r->hi + (r->lo == 0)) & LMASK;
1080 return;
1081 }
1082
1083 /* Support routines */
1084
1085 void unpackf (uint32 hi, UFP *r)
1086 {
1087 r->sign = hi & FPSIGN; /* get sign */
1088 r->exp = FD_GETEXP (hi); /* get exponent */
1089 if (r->exp == 0) { /* exp = 0? */
1090 if (r->sign) RSVD_OPND_FAULT; /* if -, rsvd op */
1091 r->frac.hi = r->frac.lo = 0; /* else 0 */
1092 return;
1093 }
1094 r->frac.hi = WORDSWAP ((hi & ~(FPSIGN | FD_EXP)) | FD_HB);
1095 r->frac.lo = 0;
1096 dp_lsh (&r->frac, FD_GUARD);
1097 return;
1098 }
1099
1100 void unpackd (uint32 hi, uint32 lo, UFP *r)
1101 {
1102 r->sign = hi & FPSIGN; /* get sign */
1103 r->exp = FD_GETEXP (hi); /* get exponent */
1104 if (r->exp == 0) { /* exp = 0? */
1105 if (r->sign) RSVD_OPND_FAULT; /* if -, rsvd op */
1106 r->frac.hi = r->frac.lo = 0; /* else 0 */
1107 return;
1108 }
1109 r->frac.hi = WORDSWAP ((hi & ~(FPSIGN | FD_EXP)) | FD_HB);
1110 r->frac.lo = WORDSWAP (lo);
1111 dp_lsh (&r->frac, FD_GUARD);
1112 return;
1113 }
1114
1115 void unpackg (uint32 hi, uint32 lo, UFP *r)
1116 {
1117 r->sign = hi & FPSIGN; /* get sign */
1118 r->exp = G_GETEXP (hi); /* get exponent */
1119 if (r->exp == 0) { /* exp = 0? */
1120 if (r->sign) RSVD_OPND_FAULT; /* if -, rsvd op */
1121 r->frac.hi = r->frac.lo = 0; /* else 0 */
1122 return;
1123 }
1124 r->frac.hi = WORDSWAP ((hi & ~(FPSIGN | G_EXP)) | G_HB);
1125 r->frac.lo = WORDSWAP (lo);
1126 dp_lsh (&r->frac, G_GUARD);
1127 return;
1128 }
1129
1130 void norm (UFP *r)
1131 {
1132 int32 i;
1133 static uint32 normmask[5] = {
1134 0xc0000000, 0xf0000000, 0xff000000, 0xffff0000, 0xffffffff
1135 };
1136 static int32 normtab[6] = { 1, 2, 4, 8, 16, 32};
1137
1138 if ((r->frac.hi == 0) && (r->frac.lo == 0)) { /* if fraction = 0 */
1139 r->sign = r->exp = 0; /* result is 0 */
1140 return;
1141 }
1142 while ((r->frac.hi & UF_NM_H) == 0) { /* normalized? */
1143 for (i = 0; i < 5; i++) { /* find first 1 */
1144 if (r->frac.hi & normmask[i]) break;
1145 }
1146 dp_lsh (&r->frac, normtab[i]); /* shift frac */
1147 r->exp = r->exp - normtab[i]; /* decr exp */
1148 }
1149 return;
1150 }
1151
1152 int32 rpackfd (UFP *r, int32 *rh)
1153 {
1154 static UDP f_round = { UF_FRND_L, UF_FRND_H };
1155 static UDP d_round = { UF_DRND_L, UF_DRND_H };
1156
1157 if (rh) *rh = 0; /* assume 0 */
1158 if ((r->frac.hi == 0) && (r->frac.lo == 0)) return 0; /* result 0? */
1159 if (rh) dp_add (&r->frac, &d_round); /* round */
1160 else dp_add (&r->frac, &f_round);
1161 if ((r->frac.hi & UF_NM_H) == 0) { /* carry out? */
1162 dp_rsh (&r->frac, 1); /* renormalize */
1163 r->exp = r->exp + 1;
1164 }
1165 if (r->exp > (int32) FD_M_EXP) FLT_OVFL_FAULT; /* ovflo? fault */
1166 if (r->exp <= 0) { /* underflow? */
1167 if (PSL & PSW_FU) FLT_UNFL_FAULT; /* fault if fu */
1168 return 0; /* else 0 */
1169 }
1170 dp_rsh (&r->frac, FD_GUARD); /* remove guard */
1171 if (rh) *rh = WORDSWAP (r->frac.lo); /* get low */
1172 return r->sign | (r->exp << FD_V_EXP) |
1173 (WORDSWAP (r->frac.hi) & ~(FD_HB | FPSIGN | FD_EXP));
1174 }
1175
1176 int32 rpackg (UFP *r, int32 *rh)
1177 {
1178 static UDP g_round = { UF_GRND_L, UF_GRND_H };
1179
1180 *rh = 0; /* assume 0 */
1181 if ((r->frac.hi == 0) && (r->frac.lo == 0)) return 0; /* result 0? */
1182 dp_add (&r->frac, &g_round); /* round */
1183 if ((r->frac.hi & UF_NM_H) == 0) { /* carry out? */
1184 dp_rsh (&r->frac, 1); /* renormalize */
1185 r->exp = r->exp + 1;
1186 }
1187 if (r->exp > (int32) G_M_EXP) FLT_OVFL_FAULT; /* ovflo? fault */
1188 if (r->exp <= 0) { /* underflow? */
1189 if (PSL & PSW_FU) FLT_UNFL_FAULT; /* fault if fu */
1190 return 0; /* else 0 */
1191 }
1192 dp_rsh (&r->frac, G_GUARD); /* remove guard */
1193 *rh = WORDSWAP (r->frac.lo); /* get low */
1194 return r->sign | (r->exp << G_V_EXP) |
1195 (WORDSWAP (r->frac.hi) & ~(G_HB | FPSIGN | G_EXP));
1196 }
1197
1198 #endif
1199
1200 /* Floating point instructions */
1201
1202 /* Move/test/move negated floating
1203
1204 Note that only the high 32b is processed.
1205 If the high 32b is not zero, it is unchanged.
1206 */
1207
1208 int32 op_movfd (int32 val)
1209 {
1210 if (val & FD_EXP) return val;
1211 if (val & FPSIGN) RSVD_OPND_FAULT;
1212 return 0;
1213 }
1214
1215 int32 op_mnegfd (int32 val)
1216 {
1217 if (val & FD_EXP) return (val ^ FPSIGN);
1218 if (val & FPSIGN) RSVD_OPND_FAULT;
1219 return 0;
1220 }
1221
1222 int32 op_movg (int32 val)
1223 {
1224 if (val & G_EXP) return val;
1225 if (val & FPSIGN) RSVD_OPND_FAULT;
1226 return 0;
1227 }
1228
1229 int32 op_mnegg (int32 val)
1230 {
1231 if (val & G_EXP) return (val ^ FPSIGN);
1232 if (val & FPSIGN) RSVD_OPND_FAULT;
1233 return 0;
1234 }
1235
1236 /* Floating to floating convert - F to D is essentially done with MOVFD */
1237
1238 int32 op_cvtdf (int32 *opnd)
1239 {
1240 UFP a;
1241
1242 unpackd (opnd[0], opnd[1], &a);
1243 return rpackfd (&a, NULL);
1244 }
1245
1246 int32 op_cvtfg (int32 *opnd, int32 *rh)
1247 {
1248 UFP a;
1249
1250 unpackf (opnd[0], &a);
1251 a.exp = a.exp - FD_BIAS + G_BIAS;
1252 return rpackg (&a, rh);
1253 }
1254
1255 int32 op_cvtgf (int32 *opnd)
1256 {
1257 UFP a;
1258
1259 unpackg (opnd[0], opnd[1], &a);
1260 a.exp = a.exp - G_BIAS + FD_BIAS;
1261 return rpackfd (&a, NULL);
1262 }
1263
1264 /* Floating add and subtract */
1265
1266 int32 op_addf (int32 *opnd, t_bool sub)
1267 {
1268 UFP a, b;
1269
1270 unpackf (opnd[0], &a); /* F format */
1271 unpackf (opnd[1], &b);
1272 if (sub) a.sign = a.sign ^ FPSIGN; /* sub? -s1 */
1273 vax_fadd (&a, &b); /* add fractions */
1274 return rpackfd (&a, NULL);
1275 }
1276
1277 int32 op_addd (int32 *opnd, int32 *rh, t_bool sub)
1278 {
1279 UFP a, b;
1280
1281 unpackd (opnd[0], opnd[1], &a);
1282 unpackd (opnd[2], opnd[3], &b);
1283 if (sub) a.sign = a.sign ^ FPSIGN; /* sub? -s1 */
1284 vax_fadd (&a, &b); /* add fractions */
1285 return rpackfd (&a, rh);
1286 }
1287
1288 int32 op_addg (int32 *opnd, int32 *rh, t_bool sub)
1289 {
1290 UFP a, b;
1291
1292 unpackg (opnd[0], opnd[1], &a);
1293 unpackg (opnd[2], opnd[3], &b);
1294 if (sub) a.sign = a.sign ^ FPSIGN; /* sub? -s1 */
1295 vax_fadd (&a, &b); /* add fractions */
1296 return rpackg (&a, rh); /* round and pack */
1297 }
1298
1299 /* Floating multiply */
1300
1301 int32 op_mulf (int32 *opnd)
1302 {
1303 UFP a, b;
1304
1305 unpackf (opnd[0], &a); /* F format */
1306 unpackf (opnd[1], &b);
1307 vax_fmul (&a, &b, 0, FD_BIAS, 0, 0); /* do multiply */
1308 return rpackfd (&a, NULL); /* round and pack */
1309 }
1310
1311 int32 op_muld (int32 *opnd, int32 *rh)
1312 {
1313 UFP a, b;
1314
1315 unpackd (opnd[0], opnd[1], &a); /* D format */
1316 unpackd (opnd[2], opnd[3], &b);
1317 vax_fmul (&a, &b, 1, FD_BIAS, 0, 0); /* do multiply */
1318 return rpackfd (&a, rh); /* round and pack */
1319 }
1320
1321 int32 op_mulg (int32 *opnd, int32 *rh)
1322 {
1323 UFP a, b;
1324
1325 unpackg (opnd[0], opnd[1], &a); /* G format */
1326 unpackg (opnd[2], opnd[3], &b);
1327 vax_fmul (&a, &b, 1, G_BIAS, 0, 0); /* do multiply */
1328 return rpackg (&a, rh); /* round and pack */
1329 }
1330
1331 /* Floating divide */
1332
1333 int32 op_divf (int32 *opnd)
1334 {
1335 UFP a, b;
1336
1337 unpackf (opnd[0], &a); /* F format */
1338 unpackf (opnd[1], &b);
1339 vax_fdiv (&a, &b, 26, FD_BIAS); /* do divide */
1340 return rpackfd (&b, NULL); /* round and pack */
1341 }
1342
1343 int32 op_divd (int32 *opnd, int32 *rh)
1344 {
1345 UFP a, b;
1346
1347 unpackd (opnd[0], opnd[1], &a); /* D format */
1348 unpackd (opnd[2], opnd[3], &b);
1349 vax_fdiv (&a, &b, 58, FD_BIAS); /* do divide */
1350 return rpackfd (&b, rh); /* round and pack */
1351 }
1352
1353 int32 op_divg (int32 *opnd, int32 *rh)
1354 {
1355 UFP a, b;
1356
1357 unpackg (opnd[0], opnd[1], &a); /* G format */
1358 unpackg (opnd[2], opnd[3], &b);
1359 vax_fdiv (&a, &b, 55, G_BIAS); /* do divide */
1360 return rpackg (&b, rh); /* round and pack */
1361 }
1362
1363 /* Polynomial evaluation
1364 The most mis-implemented instruction in the VAX (probably here too).
1365 POLY requires a precise combination of masking versus normalizing
1366 to achieve the desired answer. In particular, the multiply step
1367 is masked prior to normalization. In addition, negative small
1368 fractions must not be treated as 0 during denorm.
1369 */
1370
1371 void op_polyf (int32 *opnd, int32 acc)
1372 {
1373 UFP r, a, c;
1374 int32 deg = opnd[1];
1375 int32 ptr = opnd[2];
1376 int32 i, wd, res;
1377
1378 if (deg > 31) RSVD_OPND_FAULT; /* degree > 31? fault */
1379 unpackf (opnd[0], &a); /* unpack arg */
1380 wd = Read (ptr, L_LONG, RD); /* get C0 */
1381 ptr = ptr + 4;
1382 unpackf (wd, &r); /* unpack C0 */
1383 res = rpackfd (&r, NULL); /* first result */
1384 for (i = 0; i < deg; i++) { /* loop */
1385 unpackf (res, &r); /* unpack result */
1386 vax_fmul (&r, &a, 0, FD_BIAS, 1, LMASK); /* r = r * arg, mask */
1387 wd = Read (ptr, L_LONG, RD); /* get Cnext */
1388 ptr = ptr + 4;
1389 unpackf (wd, &c); /* unpack Cnext */
1390 vax_fadd (&r, &c); /* r = r + Cnext */
1391 res = rpackfd (&r, NULL); /* round and pack */
1392 }
1393 R[0] = res;
1394 R[1] = R[2] = 0;
1395 R[3] = ptr;
1396 return;
1397 }
1398
1399 void op_polyd (int32 *opnd, int32 acc)
1400 {
1401 UFP r, a, c;
1402 int32 deg = opnd[2];
1403 int32 ptr = opnd[3];
1404 int32 i, wd, wd1, res, resh;
1405
1406 if (deg > 31) RSVD_OPND_FAULT; /* degree > 31? fault */
1407 unpackd (opnd[0], opnd[1], &a); /* unpack arg */
1408 wd = Read (ptr, L_LONG, RD); /* get C0 */
1409 wd1 = Read (ptr + 4, L_LONG, RD);
1410 ptr = ptr + 8;
1411 unpackd (wd, wd1, &r); /* unpack C0 */
1412 res = rpackfd (&r, &resh); /* first result */
1413 for (i = 0; i < deg; i++) { /* loop */
1414 unpackd (res, resh, &r); /* unpack result */
1415 vax_fmul (&r, &a, 1, FD_BIAS, 0, 1); /* r = r * arg, mask */
1416 wd = Read (ptr, L_LONG, RD); /* get Cnext */
1417 wd1 = Read (ptr + 4, L_LONG, RD);
1418 ptr = ptr + 8;
1419 unpackd (wd, wd1, &c); /* unpack Cnext */
1420 vax_fadd (&r, &c); /* r = r + Cnext */
1421 res = rpackfd (&r, &resh); /* round and pack */
1422 }
1423 R[0] = res;
1424 R[1] = resh;
1425 R[2] = 0;
1426 R[3] = ptr;
1427 R[4] = 0;
1428 R[5] = 0;
1429 return;
1430 }
1431
1432 void op_polyg (int32 *opnd, int32 acc)
1433 {
1434 UFP r, a, c;
1435 int32 deg = opnd[2];
1436 int32 ptr = opnd[3];
1437 int32 i, wd, wd1, res, resh;
1438
1439 if (deg > 31) RSVD_OPND_FAULT; /* degree > 31? fault */
1440 unpackg (opnd[0], opnd[1], &a); /* unpack arg */
1441 wd = Read (ptr, L_LONG, RD); /* get C0 */
1442 wd1 = Read (ptr + 4, L_LONG, RD);
1443 ptr = ptr + 8;
1444 unpackg (wd, wd1, &r); /* unpack C0 */
1445 res = rpackg (&r, &resh); /* first result */
1446 for (i = 0; i < deg; i++) { /* loop */
1447 unpackg (res, resh, &r); /* unpack result */
1448 vax_fmul (&r, &a, 1, G_BIAS, 0, 1); /* r = r * arg */
1449 wd = Read (ptr, L_LONG, RD); /* get Cnext */
1450 wd1 = Read (ptr + 4, L_LONG, RD);
1451 ptr = ptr + 8;
1452 unpackg (wd, wd1, &c); /* unpack Cnext */
1453 vax_fadd (&r, &c); /* r = r + Cnext */
1454 res = rpackg (&r, &resh); /* round and pack */
1455 }
1456 R[0] = res;
1457 R[1] = resh;
1458 R[2] = 0;
1459 R[3] = ptr;
1460 R[4] = 0;
1461 R[5] = 0;
1462 return;
1463 }
Magically enhanced SIMH