Commit | Line | Data |
---|---|---|
196ba1fc PH |
1 | /* hp2100_fp.c: HP 2100 floating point instructions\r |
2 | \r | |
3 | Copyright (c) 2002-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 | 21-Jan-08 JDB Corrected fp_unpack mantissa high-word return\r | |
27 | (from Mark Pizzolato)\r | |
28 | 01-Dec-06 JDB Reworked FFP helpers for 1000-F support, deleted f_pwr2\r | |
29 | 22-Jul-05 RMS Fixed compiler warning in Solaris (from Doug Glyn)\r | |
30 | 25-Feb-05 JDB Added FFP helpers f_pack, f_unpack, f_pwr2\r | |
31 | 11-Feb-05 JDB Fixed missing negative overflow renorm in StoreFP\r | |
32 | 26-Dec-04 RMS Separated A/B from M[0/1] for DMA IO (from Dave Bryan)\r | |
33 | 15-Jul-03 RMS Fixed signed/unsigned warning\r | |
34 | 21-Oct-02 RMS Recoded for compatibility with 21MX microcode algorithms\r | |
35 | \r | |
36 | \r | |
37 | Implementation note: The 2100/1000-M/E Fast FORTRAN Processor (FFP) and 1000\r | |
38 | F-Series Floating Point Processor (FPP) simulations require that the host\r | |
39 | compiler support 64-bit integers and the HAVE_INT64 symbol be defined during\r | |
40 | compilation. If this symbol is defined, two-word floating-point operations\r | |
41 | are handled in the FPP code, and this module is not used. If it is not\r | |
42 | defined, then FFP and FPP operations are not available, and this module\r | |
43 | provides the floating-point support.\r | |
44 | \r | |
45 | \r | |
46 | The HP2100 uses a unique binary floating point format:\r | |
47 | \r | |
48 | 15 14 0\r | |
49 | +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+\r | |
50 | |S | fraction high | : A\r | |
51 | +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+\r | |
52 | | fraction low | exponent |XS| : A + 1\r | |
53 | +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+\r | |
54 | 15 8 7 1 0\r | |
55 | \r | |
56 | where S = 0 for plus fraction, 1 for minus fraction\r | |
57 | fraction = s.bbbbb..., 24 binary digits\r | |
58 | exponent = 2**+/-n\r | |
59 | XS = 0 for plus exponent, 1 for minus exponent\r | |
60 | \r | |
61 | Numbers can be normalized or unnormalized but are always normalized\r | |
62 | when loaded.\r | |
63 | \r | |
64 | Unpacked floating point numbers are stored in structure ufp\r | |
65 | \r | |
66 | exp = exponent, 2's complement\r | |
67 | h'l = fraction, 2's comp, left justified\r | |
68 | \r | |
69 | This routine tries to reproduce the algorithms of the 2100/21MX\r | |
70 | microcode in order to achieve 'bug-for-bug' compatibility. In\r | |
71 | particular,\r | |
72 | \r | |
73 | - The FIX code produces various results in B.\r | |
74 | - The fraction multiply code uses 16b x 16b multiplies to produce\r | |
75 | a 31b result. It always loses the low order bit of the product.\r | |
76 | - The fraction divide code is an approximation that may produce\r | |
77 | an error of 1 LSB.\r | |
78 | - Signs are tracked implicitly as part of the fraction. Unnormalized\r | |
79 | inputs may cause the packup code to produce the wrong sign.\r | |
80 | - "Unclean" zeros (zero fraction, non-zero exponent) are processed\r | |
81 | like normal operands.\r | |
82 | */\r | |
83 | \r | |
84 | #include "hp2100_defs.h"\r | |
85 | #include "hp2100_cpu1.h"\r | |
86 | #include "hp2100_fp.h"\r | |
87 | \r | |
88 | #if !defined (HAVE_INT64) /* int64 support unavailable */\r | |
89 | \r | |
90 | struct ufp { /* unpacked fp */\r | |
91 | int32 exp; /* exp */\r | |
92 | uint32 fr; /* frac */\r | |
93 | };\r | |
94 | \r | |
95 | #define FP_V_SIGN 31 /* sign */\r | |
96 | #define FP_M_SIGN 01\r | |
97 | #define FP_V_FR 8 /* fraction */\r | |
98 | #define FP_M_FR 077777777\r | |
99 | #define FP_V_EXP 1 /* exponent */\r | |
100 | #define FP_M_EXP 0177\r | |
101 | #define FP_V_EXPS 0 /* exp sign */\r | |
102 | #define FP_M_EXPS 01\r | |
103 | #define FP_SIGN (FP_M_SIGN << FP_V_SIGN)\r | |
104 | #define FP_FR (FP_M_FR << FP_V_FR)\r | |
105 | #define FP_EXP (FP_M_EXP << FP_V_EXP)\r | |
106 | #define FP_EXPS (FP_M_EXPS << FP_V_EXPS)\r | |
107 | #define FP_GETSIGN(x) (((x) >> FP_V_SIGN) & FP_M_SIGN)\r | |
108 | #define FP_GETEXP(x) (((x) >> FP_V_EXP) & FP_M_EXP)\r | |
109 | #define FP_GETEXPS(x) (((x) >> FP_V_EXPS) & FP_M_EXPS)\r | |
110 | \r | |
111 | #define FP_NORM (1 << (FP_V_SIGN - 1)) /* normalized */\r | |
112 | #define FP_LOW (1 << FP_V_FR)\r | |
113 | #define FP_RNDP (1 << (FP_V_FR - 1)) /* round for plus */\r | |
114 | #define FP_RNDM (FP_RNDP - 1) /* round for minus */\r | |
115 | \r | |
116 | #define FPAB ((((uint32) AR) << 16) | ((uint32) BR))\r | |
117 | \r | |
118 | /* Fraction shift; 0 < shift < 32 */\r | |
119 | \r | |
120 | #define FR_ARS(v,s) (((v) >> (s)) | (((v) & FP_SIGN)? \\r | |
121 | (((uint32) DMASK32) << (32 - (s))): 0)) & DMASK32\r | |
122 | \r | |
123 | #define FR_NEG(v) ((~(v) + 1) & DMASK32)\r | |
124 | \r | |
125 | extern uint16 ABREG[2];\r | |
126 | \r | |
127 | uint32 UnpackFP (struct ufp *fop, uint32 opnd);\r | |
128 | void NegFP (struct ufp *fop);\r | |
129 | void NormFP (struct ufp *fop);\r | |
130 | uint32 PackFP (struct ufp *fop);\r | |
131 | uint32 StoreFP (struct ufp *fop);\r | |
132 | \r | |
133 | /* Floating to integer conversion */\r | |
134 | \r | |
135 | uint32 f_fix (void)\r | |
136 | {\r | |
137 | struct ufp fop;\r | |
138 | uint32 res = 0;\r | |
139 | \r | |
140 | UnpackFP (&fop, FPAB); /* unpack op */\r | |
141 | if (fop.exp < 0) { /* exp < 0? */\r | |
142 | AR = 0; /* result = 0 */\r | |
143 | return 0; /* B unchanged */\r | |
144 | }\r | |
145 | if (fop.exp > 15) { /* exp > 15? */\r | |
146 | BR = AR; /* B has high bits */\r | |
147 | AR = 077777; /* result = 77777 */\r | |
148 | return 1; /* overflow */\r | |
149 | }\r | |
150 | if (fop.exp < 15) { /* if not aligned */\r | |
151 | res = FR_ARS (fop.fr, 15 - fop.exp); /* shift right */\r | |
152 | AR = (res >> 16) & DMASK; /* AR gets result */\r | |
153 | }\r | |
154 | BR = AR;\r | |
155 | if ((AR & SIGN) && ((fop.fr | res) & DMASK)) /* any low bits lost? */\r | |
156 | AR = (AR + 1) & DMASK; /* round up */\r | |
157 | return 0;\r | |
158 | }\r | |
159 | \r | |
160 | /* Integer to floating conversion */\r | |
161 | \r | |
162 | uint32 f_flt (void)\r | |
163 | {\r | |
164 | struct ufp res = { 15, 0 }; /* +, 2**15 */\r | |
165 | \r | |
166 | res.fr = ((uint32) AR) << 16; /* left justify */\r | |
167 | StoreFP (&res); /* store result */\r | |
168 | return 0; /* clr overflow */\r | |
169 | }\r | |
170 | \r | |
171 | /* Floating point add/subtract */\r | |
172 | \r | |
173 | uint32 f_as (uint32 opnd, t_bool sub)\r | |
174 | {\r | |
175 | struct ufp fop1, fop2, t;\r | |
176 | int32 ediff;\r | |
177 | \r | |
178 | UnpackFP (&fop1, FPAB); /* unpack A-B */\r | |
179 | UnpackFP (&fop2, opnd); /* get op */\r | |
180 | if (sub) { /* subtract? */\r | |
181 | fop2.fr = FR_NEG (fop2.fr); /* negate frac */\r | |
182 | if (fop2.fr == ((uint32) FP_SIGN)) { /* -1/2? */\r | |
183 | fop2.fr = fop2.fr >> 1; /* special case */\r | |
184 | fop2.exp = fop2.exp + 1;\r | |
185 | }\r | |
186 | }\r | |
187 | if (fop1.fr == 0) fop1 = fop2; /* op1 = 0? res = op2 */\r | |
188 | else if (fop2.fr != 0) { /* op2 = 0? no add */\r | |
189 | if (fop1.exp < fop2.exp) { /* |op1| < |op2|? */\r | |
190 | t = fop2; /* swap operands */\r | |
191 | fop2 = fop1;\r | |
192 | fop1 = t;\r | |
193 | }\r | |
194 | ediff = fop1.exp - fop2.exp; /* get exp diff */\r | |
195 | if (ediff <= 24) {\r | |
196 | if (ediff) fop2.fr = FR_ARS (fop2.fr, ediff); /* denorm, signed */\r | |
197 | if ((fop1.fr ^ fop2.fr) & FP_SIGN) /* unlike signs? */\r | |
198 | fop1.fr = fop1.fr + fop2.fr; /* eff subtract */\r | |
199 | else { /* like signs */\r | |
200 | fop1.fr = fop1.fr + fop2.fr; /* eff add */\r | |
201 | if (fop2.fr & FP_SIGN) { /* both -? */\r | |
202 | if ((fop1.fr & FP_SIGN) == 0) { /* overflow? */\r | |
203 | fop1.fr = FP_SIGN | (fop1.fr >> 1); /* renormalize */\r | |
204 | fop1.exp = fop1.exp + 1; /* incr exp */\r | |
205 | }\r | |
206 | }\r | |
207 | else if (fop1.fr & FP_SIGN) { /* both +, cry out? */\r | |
208 | fop1.fr = fop1.fr >> 1; /* renormalize */\r | |
209 | fop1.exp = fop1.exp + 1; /* incr exp */\r | |
210 | }\r | |
211 | } /* end else like */\r | |
212 | } /* end if ediff */\r | |
213 | } /* end if fop2 */\r | |
214 | return StoreFP (&fop1); /* store result */\r | |
215 | }\r | |
216 | \r | |
217 | /* Floating point multiply - passes diagnostic */\r | |
218 | \r | |
219 | uint32 f_mul (uint32 opnd)\r | |
220 | {\r | |
221 | struct ufp fop1, fop2;\r | |
222 | struct ufp res = { 0, 0 };\r | |
223 | int32 shi1, shi2, t1, t2, t3, t4, t5;\r | |
224 | \r | |
225 | UnpackFP (&fop1, FPAB); /* unpack A-B */\r | |
226 | UnpackFP (&fop2, opnd); /* unpack op */\r | |
227 | if (fop1.fr && fop2.fr) { /* if both != 0 */\r | |
228 | res.exp = fop1.exp + fop2.exp + 1; /* exp = sum */\r | |
229 | shi1 = SEXT (fop1.fr >> 16); /* mpy hi */\r | |
230 | shi2 = SEXT (fop2.fr >> 16); /* mpc hi */\r | |
231 | t1 = shi2 * ((int32) ((fop1.fr >> 1) & 077600)); /* mpc hi * (mpy lo/2) */\r | |
232 | t2 = shi1 * ((int32) ((fop2.fr >> 1) & 077600)); /* mpc lo * (mpy hi/2) */\r | |
233 | t3 = t1 + t2; /* cross product */\r | |
234 | t4 = (shi1 * shi2) & ~1; /* mpy hi * mpc hi */\r | |
235 | t5 = (SEXT (t3 >> 16)) << 1; /* add in cross */\r | |
236 | res.fr = (t4 + t5) & DMASK32; /* bit<0> is lost */\r | |
237 | }\r | |
238 | return StoreFP (&res); /* store */\r | |
239 | }\r | |
240 | \r | |
241 | /* Floating point divide - reverse engineered from diagnostic */\r | |
242 | \r | |
243 | uint32 divx (uint32 ba, uint32 dvr, uint32 *rem)\r | |
244 | {\r | |
245 | int32 sdvd = 0, sdvr = 0;\r | |
246 | uint32 q, r;\r | |
247 | \r | |
248 | if (ba & FP_SIGN) sdvd = 1; /* 32b/16b signed dvd */\r | |
249 | if (dvr & SIGN) sdvr = 1; /* use old-fashioned */\r | |
250 | if (sdvd) ba = (~ba + 1) & DMASK32; /* unsigned divides, */\r | |
251 | if (sdvr) dvr = (~dvr + 1) & DMASK; /* as results may ovflo */\r | |
252 | q = ba / dvr;\r | |
253 | r = ba % dvr;\r | |
254 | if (sdvd ^ sdvr) q = (~q + 1) & DMASK;\r | |
255 | if (sdvd) r = (~r + 1) & DMASK;\r | |
256 | if (rem) *rem = r;\r | |
257 | return q;\r | |
258 | }\r | |
259 | \r | |
260 | uint32 f_div (uint32 opnd)\r | |
261 | {\r | |
262 | struct ufp fop1, fop2;\r | |
263 | struct ufp quo = { 0, 0 };\r | |
264 | uint32 ba, q0, q1, q2, dvrh;\r | |
265 | \r | |
266 | UnpackFP (&fop1, FPAB); /* unpack A-B */\r | |
267 | UnpackFP (&fop2, opnd); /* unpack op */\r | |
268 | dvrh = (fop2.fr >> 16) & DMASK; /* high divisor */\r | |
269 | if (dvrh == 0) { /* div by zero? */\r | |
270 | AR = 0077777; /* return most pos */\r | |
271 | BR = 0177776;\r | |
272 | return 1;\r | |
273 | }\r | |
274 | if (fop1.fr) { /* dvd != 0? */\r | |
275 | quo.exp = fop1.exp - fop2.exp + 1; /* exp = diff */\r | |
276 | ba = FR_ARS (fop1.fr, 2); /* prevent ovflo */\r | |
277 | q0 = divx (ba, dvrh, &ba); /* Q0 = dvd / dvrh */\r | |
278 | ba = (ba & ~1) << 16; /* remainder */\r | |
279 | ba = FR_ARS (ba, 1); /* prevent ovflo */\r | |
280 | q1 = divx (ba, dvrh, NULL); /* Q1 = rem / dvrh */\r | |
281 | ba = (fop2.fr & 0xFF00) << 13; /* dvrl / 8 */\r | |
282 | q2 = divx (ba, dvrh, NULL); /* dvrl / dvrh */\r | |
283 | ba = -(SEXT (q2)) * (SEXT (q0)); /* -Q0 * Q2 */\r | |
284 | ba = (ba >> 16) & 0xFFFF; /* save ms half */\r | |
285 | if (q1 & SIGN) quo.fr = quo.fr - 0x00010000; /* Q1 < 0? -1 */\r | |
286 | if (ba & SIGN) quo.fr = quo.fr - 0x00010000; /* -Q0*Q2 < 0? */\r | |
287 | quo.fr = quo.fr + ((ba << 2) & 0xFFFF) + q1; /* rest prod, add Q1 */\r | |
288 | quo.fr = quo.fr << 1; /* shift result */\r | |
289 | quo.fr = quo.fr + (q0 << 16); /* add Q0 */\r | |
290 | } /* end if fop1.h */\r | |
291 | return StoreFP (&quo); /* store result */\r | |
292 | }\r | |
293 | \r | |
294 | /* Utility routines */\r | |
295 | \r | |
296 | /* Unpack operand */\r | |
297 | \r | |
298 | uint32 UnpackFP (struct ufp *fop, uint32 opnd)\r | |
299 | {\r | |
300 | fop->fr = opnd & FP_FR; /* get frac */\r | |
301 | fop->exp = FP_GETEXP (opnd); /* get exp */\r | |
302 | if (FP_GETEXPS (opnd)) fop->exp = fop->exp | ~FP_M_EXP; /* < 0? sext */\r | |
303 | return FP_GETSIGN (opnd); /* return sign */\r | |
304 | }\r | |
305 | \r | |
306 | /* Normalize unpacked floating point number */\r | |
307 | \r | |
308 | void NormFP (struct ufp *fop)\r | |
309 | {\r | |
310 | if (fop->fr) { /* any fraction? */\r | |
311 | uint32 test = (fop->fr >> 1) & FP_NORM;\r | |
312 | while ((fop->fr & FP_NORM) == test) { /* until norm */\r | |
313 | fop->exp = fop->exp - 1;\r | |
314 | fop->fr = (fop->fr << 1);\r | |
315 | }\r | |
316 | }\r | |
317 | else fop->exp = 0; /* clean 0 */\r | |
318 | return;\r | |
319 | }\r | |
320 | \r | |
321 | /* Pack fp number */\r | |
322 | \r | |
323 | uint32 PackFP (struct ufp *fop)\r | |
324 | {\r | |
325 | return (fop->fr & FP_FR) | /* merge frac */\r | |
326 | ((fop->exp & FP_M_EXP) << FP_V_EXP) | /* and exp */\r | |
327 | ((fop->exp < 0)? (1 << FP_V_EXPS): 0); /* add exp sign */\r | |
328 | }\r | |
329 | \r | |
330 | /* Round fp number, store, generate overflow */\r | |
331 | \r | |
332 | uint32 StoreFP (struct ufp *fop)\r | |
333 | {\r | |
334 | uint32 sign, svfr, hi, ov = 0;\r | |
335 | \r | |
336 | NormFP (fop); /* normalize */\r | |
337 | svfr = fop->fr; /* save fraction */\r | |
338 | sign = FP_GETSIGN (fop->fr); /* save sign */\r | |
339 | fop->fr = (fop->fr + (sign? FP_RNDM: FP_RNDP)) & FP_FR; /* round */\r | |
340 | if ((fop->fr ^ svfr) & FP_SIGN) { /* sign change? */\r | |
341 | fop->fr = fop->fr >> 1; /* renormalize */\r | |
342 | fop->exp = fop->exp + 1;\r | |
343 | }\r | |
344 | else NormFP (fop); /* check for norm */\r | |
345 | if (fop->fr == 0) hi = 0; /* result 0? */\r | |
346 | else if (fop->exp < -(FP_M_EXP + 1)) { /* underflow? */\r | |
347 | hi = 0; /* store clean 0 */\r | |
348 | ov = 1;\r | |
349 | }\r | |
350 | else if (fop->exp > FP_M_EXP) { /* overflow? */\r | |
351 | hi = 0x7FFFFFFE; /* all 1's */\r | |
352 | ov = 1;\r | |
353 | }\r | |
354 | else hi = PackFP (fop); /* pack mant and exp */\r | |
355 | AR = (hi >> 16) & DMASK;\r | |
356 | BR = hi & DMASK;\r | |
357 | return ov;\r | |
358 | }\r | |
359 | \r | |
360 | \r | |
361 | /* Single-precision Fast FORTRAN Processor helpers. */\r | |
362 | \r | |
363 | /* Pack mantissa and exponent and return fp value. */\r | |
364 | \r | |
365 | uint32 fp_pack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)\r | |
366 | {\r | |
367 | struct ufp fop;\r | |
368 | uint32 val;\r | |
369 | \r | |
370 | fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];\r | |
371 | fop.exp = exponent;\r | |
372 | val = PackFP (&fop);\r | |
373 | result->fpk[0] = (int16) (val >> 16);\r | |
374 | result->fpk[1] = (int16) val;\r | |
375 | return 0;\r | |
376 | }\r | |
377 | \r | |
378 | /* Normalize, round, and pack mantissa and exponent and return fp value. */\r | |
379 | \r | |
380 | uint32 fp_nrpack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)\r | |
381 | {\r | |
382 | struct ufp fop;\r | |
383 | uint32 ovf;\r | |
384 | \r | |
385 | fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];\r | |
386 | fop.exp = exponent;\r | |
387 | ovf = StoreFP (&fop);\r | |
388 | result->fpk[0] = AR;\r | |
389 | result->fpk[1] = BR;\r | |
390 | return ovf;\r | |
391 | }\r | |
392 | \r | |
393 | /* Unpack fp number in into mantissa and exponent. */\r | |
394 | \r | |
395 | uint32 fp_unpack (OP *mantissa, int32 *exponent, OP packed, OPSIZE precision)\r | |
396 | {\r | |
397 | struct ufp fop;\r | |
398 | uint32 operand;\r | |
399 | \r | |
400 | operand = ((uint32) packed.fpk[0] << 16) | packed.fpk[1];\r | |
401 | UnpackFP (&fop, operand);\r | |
402 | mantissa->fpk[0] = (uint16) (fop.fr >> 16);\r | |
403 | mantissa->fpk[1] = (uint16) fop.fr;\r | |
404 | *exponent = fop.exp;\r | |
405 | return 0;\r | |
406 | }\r | |
407 | \r | |
408 | #endif /* int64 support unavailable */\r |