projects / simh.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
First Commit of my working state
[simh.git] / HP2100 / hp2100_fp.c
CommitLineData
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
90struct 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
125extern uint16 ABREG[2];\r
126\r
127uint32 UnpackFP (struct ufp *fop, uint32 opnd);\r
128void NegFP (struct ufp *fop);\r
129void NormFP (struct ufp *fop);\r
130uint32 PackFP (struct ufp *fop);\r
131uint32 StoreFP (struct ufp *fop);\r
132\r
133/* Floating to integer conversion */\r
134\r
135uint32 f_fix (void)\r
136{\r
137struct ufp fop;\r
138uint32 res = 0;\r
139\r
140UnpackFP (&fop, FPAB); /* unpack op */\r
141if (fop.exp < 0) { /* exp < 0? */\r
142 AR = 0; /* result = 0 */\r
143 return 0; /* B unchanged */\r
144 }\r
145if (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
150if (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
154BR = AR;\r
155if ((AR & SIGN) && ((fop.fr | res) & DMASK)) /* any low bits lost? */\r
156 AR = (AR + 1) & DMASK; /* round up */\r
157return 0;\r
158}\r
159\r
160/* Integer to floating conversion */\r
161\r
162uint32 f_flt (void)\r
163{\r
164struct ufp res = { 15, 0 }; /* +, 2**15 */\r
165\r
166res.fr = ((uint32) AR) << 16; /* left justify */\r
167StoreFP (&res); /* store result */\r
168return 0; /* clr overflow */\r
169}\r
170\r
171/* Floating point add/subtract */\r
172\r
173uint32 f_as (uint32 opnd, t_bool sub)\r
174{\r
175struct ufp fop1, fop2, t;\r
176int32 ediff;\r
177\r
178UnpackFP (&fop1, FPAB); /* unpack A-B */\r
179UnpackFP (&fop2, opnd); /* get op */\r
180if (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
187if (fop1.fr == 0) fop1 = fop2; /* op1 = 0? res = op2 */\r
188else 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
214return StoreFP (&fop1); /* store result */\r
215}\r
216\r
217/* Floating point multiply - passes diagnostic */\r
218\r
219uint32 f_mul (uint32 opnd)\r
220{\r
221struct ufp fop1, fop2;\r
222struct ufp res = { 0, 0 };\r
223int32 shi1, shi2, t1, t2, t3, t4, t5;\r
224\r
225UnpackFP (&fop1, FPAB); /* unpack A-B */\r
226UnpackFP (&fop2, opnd); /* unpack op */\r
227if (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
238return StoreFP (&res); /* store */\r
239}\r
240\r
241/* Floating point divide - reverse engineered from diagnostic */\r
242\r
243uint32 divx (uint32 ba, uint32 dvr, uint32 *rem)\r
244{\r
245int32 sdvd = 0, sdvr = 0;\r
246uint32 q, r;\r
247\r
248if (ba & FP_SIGN) sdvd = 1; /* 32b/16b signed dvd */\r
249if (dvr & SIGN) sdvr = 1; /* use old-fashioned */\r
250if (sdvd) ba = (~ba + 1) & DMASK32; /* unsigned divides, */\r
251if (sdvr) dvr = (~dvr + 1) & DMASK; /* as results may ovflo */\r
252q = ba / dvr;\r
253r = ba % dvr;\r
254if (sdvd ^ sdvr) q = (~q + 1) & DMASK;\r
255if (sdvd) r = (~r + 1) & DMASK;\r
256if (rem) *rem = r;\r
257return q;\r
258}\r
259\r
260uint32 f_div (uint32 opnd)\r
261{\r
262struct ufp fop1, fop2;\r
263struct ufp quo = { 0, 0 };\r
264uint32 ba, q0, q1, q2, dvrh;\r
265\r
266UnpackFP (&fop1, FPAB); /* unpack A-B */\r
267UnpackFP (&fop2, opnd); /* unpack op */\r
268dvrh = (fop2.fr >> 16) & DMASK; /* high divisor */\r
269if (dvrh == 0) { /* div by zero? */\r
270 AR = 0077777; /* return most pos */\r
271 BR = 0177776;\r
272 return 1;\r
273 }\r
274if (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
291return StoreFP (&quo); /* store result */\r
292}\r
293\r
294/* Utility routines */\r
295\r
296/* Unpack operand */\r
297\r
298uint32 UnpackFP (struct ufp *fop, uint32 opnd)\r
299{\r
300fop->fr = opnd & FP_FR; /* get frac */\r
301fop->exp = FP_GETEXP (opnd); /* get exp */\r
302if (FP_GETEXPS (opnd)) fop->exp = fop->exp | ~FP_M_EXP; /* < 0? sext */\r
303return FP_GETSIGN (opnd); /* return sign */\r
304}\r
305\r
306/* Normalize unpacked floating point number */\r
307\r
308void NormFP (struct ufp *fop)\r
309{\r
310if (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
317else fop->exp = 0; /* clean 0 */\r
318return;\r
319}\r
320\r
321/* Pack fp number */\r
322\r
323uint32 PackFP (struct ufp *fop)\r
324{\r
325return (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
332uint32 StoreFP (struct ufp *fop)\r
333{\r
334uint32 sign, svfr, hi, ov = 0;\r
335\r
336NormFP (fop); /* normalize */\r
337svfr = fop->fr; /* save fraction */\r
338sign = FP_GETSIGN (fop->fr); /* save sign */\r
339fop->fr = (fop->fr + (sign? FP_RNDM: FP_RNDP)) & FP_FR; /* round */\r
340if ((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
344else NormFP (fop); /* check for norm */\r
345if (fop->fr == 0) hi = 0; /* result 0? */\r
346else if (fop->exp < -(FP_M_EXP + 1)) { /* underflow? */\r
347 hi = 0; /* store clean 0 */\r
348 ov = 1;\r
349 }\r
350else if (fop->exp > FP_M_EXP) { /* overflow? */\r
351 hi = 0x7FFFFFFE; /* all 1's */\r
352 ov = 1;\r
353 }\r
354else hi = PackFP (fop); /* pack mant and exp */\r
355AR = (hi >> 16) & DMASK;\r
356BR = hi & DMASK;\r
357return 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
365uint32 fp_pack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)\r
366{\r
367struct ufp fop;\r
368uint32 val;\r
369\r
370fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];\r
371fop.exp = exponent;\r
372val = PackFP (&fop);\r
373result->fpk[0] = (int16) (val >> 16);\r
374result->fpk[1] = (int16) val;\r
375return 0;\r
376}\r
377\r
378/* Normalize, round, and pack mantissa and exponent and return fp value. */\r
379\r
380uint32 fp_nrpack (OP *result, OP mantissa, int32 exponent, OPSIZE precision)\r
381{\r
382struct ufp fop;\r
383uint32 ovf;\r
384\r
385fop.fr = ((uint32) mantissa.fpk[0] << 16) | mantissa.fpk[1];\r
386fop.exp = exponent;\r
387ovf = StoreFP (&fop);\r
388result->fpk[0] = AR;\r
389result->fpk[1] = BR;\r
390return ovf;\r
391}\r
392\r
393/* Unpack fp number in into mantissa and exponent. */\r
394\r
395uint32 fp_unpack (OP *mantissa, int32 *exponent, OP packed, OPSIZE precision)\r
396{\r
397struct ufp fop;\r
398uint32 operand;\r
399\r
400operand = ((uint32) packed.fpk[0] << 16) | packed.fpk[1];\r
401UnpackFP (&fop, operand);\r
402mantissa->fpk[0] = (uint16) (fop.fr >> 16);\r
403mantissa->fpk[1] = (uint16) fop.fr;\r
404*exponent = fop.exp;\r
405return 0;\r
406}\r
407\r
408#endif /* int64 support unavailable */\r
Magically enhanced SIMH