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