1 /* hp2100_fp.c: HP 2100 floating point instructions
3 Copyright (c) 2002-2008, Robert M. Supnik
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:
12 The above copyright notice and this permission notice shall be included in
13 all copies or substantial portions of the Software.
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.
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.
26 21-Jan-08 JDB Corrected fp_unpack mantissa high-word return
28 01-Dec-06 JDB Reworked FFP helpers for 1000-F support, deleted f_pwr2
29 22-Jul-05 RMS Fixed compiler warning in Solaris (from Doug Glyn)
30 25-Feb-05 JDB Added FFP helpers f_pack, f_unpack, f_pwr2
31 11-Feb-05 JDB Fixed missing negative overflow renorm in StoreFP
32 26-Dec-04 RMS Separated A/B from M[0/1] for DMA IO (from Dave Bryan)
33 15-Jul-03 RMS Fixed signed/unsigned warning
34 21-Oct-02 RMS Recoded for compatibility with 21MX microcode algorithms
37 Implementation note: The 2100/1000-M/E Fast FORTRAN Processor (FFP) and 1000
38 F-Series Floating Point Processor (FPP) simulations require that the host
39 compiler support 64-bit integers and the HAVE_INT64 symbol be defined during
40 compilation. If this symbol is defined, two-word floating-point operations
41 are handled in the FPP code, and this module is not used. If it is not
42 defined, then FFP and FPP operations are not available, and this module
43 provides the floating-point support.
46 The HP2100 uses a unique binary floating point format:
49 +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
50 |S | fraction high | : A
51 +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
52 | fraction low | exponent |XS| : A + 1
53 +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+
56 where S = 0 for plus fraction, 1 for minus fraction
57 fraction = s.bbbbb..., 24 binary digits
59 XS = 0 for plus exponent, 1 for minus exponent
61 Numbers can be normalized or unnormalized but are always normalized
64 Unpacked floating point numbers are stored in structure ufp
66 exp = exponent, 2's complement
67 h'l = fraction, 2's comp, left justified
69 This routine tries to reproduce the algorithms of the 2100/21MX
70 microcode in order to achieve 'bug-for-bug' compatibility. In
73 - The FIX code produces various results in B.
74 - The fraction multiply code uses 16b x 16b multiplies to produce
75 a 31b result. It always loses the low order bit of the product.
76 - The fraction divide code is an approximation that may produce
78 - Signs are tracked implicitly as part of the fraction. Unnormalized
79 inputs may cause the packup code to produce the wrong sign.
80 - "Unclean" zeros (zero fraction, non-zero exponent) are processed
84 #include "hp2100_defs.h"
85 #include "hp2100_cpu1.h"
86 #include "hp2100_fp.h"
88 #if !defined (HAVE_INT64) /* int64 support unavailable */
90 struct ufp
{ /* unpacked fp */
95 #define FP_V_SIGN 31 /* sign */
97 #define FP_V_FR 8 /* fraction */
98 #define FP_M_FR 077777777
99 #define FP_V_EXP 1 /* exponent */
100 #define FP_M_EXP 0177
101 #define FP_V_EXPS 0 /* exp sign */
103 #define FP_SIGN (FP_M_SIGN << FP_V_SIGN)
104 #define FP_FR (FP_M_FR << FP_V_FR)
105 #define FP_EXP (FP_M_EXP << FP_V_EXP)
106 #define FP_EXPS (FP_M_EXPS << FP_V_EXPS)
107 #define FP_GETSIGN(x) (((x) >> FP_V_SIGN) & FP_M_SIGN)
108 #define FP_GETEXP(x) (((x) >> FP_V_EXP) & FP_M_EXP)
109 #define FP_GETEXPS(x) (((x) >> FP_V_EXPS) & FP_M_EXPS)
111 #define FP_NORM (1 << (FP_V_SIGN - 1)) /* normalized */
112 #define FP_LOW (1 << FP_V_FR)
113 #define FP_RNDP (1 << (FP_V_FR - 1)) /* round for plus */
114 #define FP_RNDM (FP_RNDP - 1) /* round for minus */
116 #define FPAB ((((uint32) AR) << 16) | ((uint32) BR))
118 /* Fraction shift; 0 < shift < 32 */
120 #define FR_ARS(v,s) (((v) >> (s)) | (((v) & FP_SIGN)? \
121 (((uint32) DMASK32) << (32 - (s))): 0)) & DMASK32
123 #define FR_NEG(v) ((~(v) + 1) & DMASK32)
125 extern uint16 ABREG
[2];
127 uint32
UnpackFP (struct ufp
*fop
, uint32 opnd
);
128 void NegFP (struct ufp
*fop
);
129 void NormFP (struct ufp
*fop
);
130 uint32
PackFP (struct ufp
*fop
);
131 uint32
StoreFP (struct ufp
*fop
);
133 /* Floating to integer conversion */
140 UnpackFP (&fop
, FPAB
); /* unpack op */
141 if (fop
.exp
< 0) { /* exp < 0? */
142 AR
= 0; /* result = 0 */
143 return 0; /* B unchanged */
145 if (fop
.exp
> 15) { /* exp > 15? */
146 BR
= AR
; /* B has high bits */
147 AR
= 077777; /* result = 77777 */
148 return 1; /* overflow */
150 if (fop
.exp
< 15) { /* if not aligned */
151 res
= FR_ARS (fop
.fr
, 15 - fop
.exp
); /* shift right */
152 AR
= (res
>> 16) & DMASK
; /* AR gets result */
155 if ((AR
& SIGN
) && ((fop
.fr
| res
) & DMASK
)) /* any low bits lost? */
156 AR
= (AR
+ 1) & DMASK
; /* round up */
160 /* Integer to floating conversion */
164 struct ufp res
= { 15, 0 }; /* +, 2**15 */
166 res
.fr
= ((uint32
) AR
) << 16; /* left justify */
167 StoreFP (&res
); /* store result */
168 return 0; /* clr overflow */
171 /* Floating point add/subtract */
173 uint32
f_as (uint32 opnd
, t_bool sub
)
175 struct ufp fop1
, fop2
, t
;
178 UnpackFP (&fop1
, FPAB
); /* unpack A-B */
179 UnpackFP (&fop2
, opnd
); /* get op */
180 if (sub
) { /* subtract? */
181 fop2
.fr
= FR_NEG (fop2
.fr
); /* negate frac */
182 if (fop2
.fr
== ((uint32
) FP_SIGN
)) { /* -1/2? */
183 fop2
.fr
= fop2
.fr
>> 1; /* special case */
184 fop2
.exp
= fop2
.exp
+ 1;
187 if (fop1
.fr
== 0) fop1
= fop2
; /* op1 = 0? res = op2 */
188 else if (fop2
.fr
!= 0) { /* op2 = 0? no add */
189 if (fop1
.exp
< fop2
.exp
) { /* |op1| < |op2|? */
190 t
= fop2
; /* swap operands */
194 ediff
= fop1
.exp
- fop2
.exp
; /* get exp diff */
196 if (ediff
) fop2
.fr
= FR_ARS (fop2
.fr
, ediff
); /* denorm, signed */
197 if ((fop1
.fr
^ fop2
.fr
) & FP_SIGN
) /* unlike signs? */
198 fop1
.fr
= fop1
.fr
+ fop2
.fr
; /* eff subtract */
199 else { /* like signs */
200 fop1
.fr
= fop1
.fr
+ fop2
.fr
; /* eff add */
201 if (fop2
.fr
& FP_SIGN
) { /* both -? */
202 if ((fop1
.fr
& FP_SIGN
) == 0) { /* overflow? */
203 fop1
.fr
= FP_SIGN
| (fop1
.fr
>> 1); /* renormalize */
204 fop1
.exp
= fop1
.exp
+ 1; /* incr exp */
207 else if (fop1
.fr
& FP_SIGN
) { /* both +, cry out? */
208 fop1
.fr
= fop1
.fr
>> 1; /* renormalize */
209 fop1
.exp
= fop1
.exp
+ 1; /* incr exp */
211 } /* end else like */
214 return StoreFP (&fop1
); /* store result */
217 /* Floating point multiply - passes diagnostic */
219 uint32
f_mul (uint32 opnd
)
221 struct ufp fop1
, fop2
;
222 struct ufp res
= { 0, 0 };
223 int32 shi1
, shi2
, t1
, t2
, t3
, t4
, t5
;
225 UnpackFP (&fop1
, FPAB
); /* unpack A-B */
226 UnpackFP (&fop2
, opnd
); /* unpack op */
227 if (fop1
.fr
&& fop2
.fr
) { /* if both != 0 */
228 res
.exp
= fop1
.exp
+ fop2
.exp
+ 1; /* exp = sum */
229 shi1
= SEXT (fop1
.fr
>> 16); /* mpy hi */
230 shi2
= SEXT (fop2
.fr
>> 16); /* mpc hi */
231 t1
= shi2
* ((int32
) ((fop1
.fr
>> 1) & 077600)); /* mpc hi * (mpy lo/2) */
232 t2
= shi1
* ((int32
) ((fop2
.fr
>> 1) & 077600)); /* mpc lo * (mpy hi/2) */
233 t3
= t1
+ t2
; /* cross product */
234 t4
= (shi1
* shi2
) & ~1; /* mpy hi * mpc hi */
235 t5
= (SEXT (t3
>> 16)) << 1; /* add in cross */
236 res
.fr
= (t4
+ t5
) & DMASK32
; /* bit<0> is lost */
238 return StoreFP (&res
); /* store */
241 /* Floating point divide - reverse engineered from diagnostic */
243 uint32
divx (uint32 ba
, uint32 dvr
, uint32
*rem
)
245 int32 sdvd
= 0, sdvr
= 0;
248 if (ba
& FP_SIGN
) sdvd
= 1; /* 32b/16b signed dvd */
249 if (dvr
& SIGN
) sdvr
= 1; /* use old-fashioned */
250 if (sdvd
) ba
= (~ba
+ 1) & DMASK32
; /* unsigned divides, */
251 if (sdvr
) dvr
= (~dvr
+ 1) & DMASK
; /* as results may ovflo */
254 if (sdvd
^ sdvr
) q
= (~q
+ 1) & DMASK
;
255 if (sdvd
) r
= (~r
+ 1) & DMASK
;
260 uint32
f_div (uint32 opnd
)
262 struct ufp fop1
, fop2
;
263 struct ufp quo
= { 0, 0 };
264 uint32 ba
, q0
, q1
, q2
, dvrh
;
266 UnpackFP (&fop1
, FPAB
); /* unpack A-B */
267 UnpackFP (&fop2
, opnd
); /* unpack op */
268 dvrh
= (fop2
.fr
>> 16) & DMASK
; /* high divisor */
269 if (dvrh
== 0) { /* div by zero? */
270 AR
= 0077777; /* return most pos */
274 if (fop1
.fr
) { /* dvd != 0? */
275 quo
.exp
= fop1
.exp
- fop2
.exp
+ 1; /* exp = diff */
276 ba
= FR_ARS (fop1
.fr
, 2); /* prevent ovflo */
277 q0
= divx (ba
, dvrh
, &ba
); /* Q0 = dvd / dvrh */
278 ba
= (ba
& ~1) << 16; /* remainder */
279 ba
= FR_ARS (ba
, 1); /* prevent ovflo */
280 q1
= divx (ba
, dvrh
, NULL
); /* Q1 = rem / dvrh */
281 ba
= (fop2
.fr
& 0xFF00) << 13; /* dvrl / 8 */
282 q2
= divx (ba
, dvrh
, NULL
); /* dvrl / dvrh */
283 ba
= -(SEXT (q2
)) * (SEXT (q0
)); /* -Q0 * Q2 */
284 ba
= (ba
>> 16) & 0xFFFF; /* save ms half */
285 if (q1
& SIGN
) quo
.fr
= quo
.fr
- 0x00010000; /* Q1 < 0? -1 */
286 if (ba
& SIGN
) quo
.fr
= quo
.fr
- 0x00010000; /* -Q0*Q2 < 0? */
287 quo
.fr
= quo
.fr
+ ((ba
<< 2) & 0xFFFF) + q1
; /* rest prod, add Q1 */
288 quo
.fr
= quo
.fr
<< 1; /* shift result */
289 quo
.fr
= quo
.fr
+ (q0
<< 16); /* add Q0 */
290 } /* end if fop1.h */
291 return StoreFP (&quo
); /* store result */
294 /* Utility routines */
298 uint32
UnpackFP (struct ufp
*fop
, uint32 opnd
)
300 fop
->fr
= opnd
& FP_FR
; /* get frac */
301 fop
->exp
= FP_GETEXP (opnd
); /* get exp */
302 if (FP_GETEXPS (opnd
)) fop
->exp
= fop
->exp
| ~FP_M_EXP
; /* < 0? sext */
303 return FP_GETSIGN (opnd
); /* return sign */
306 /* Normalize unpacked floating point number */
308 void NormFP (struct ufp
*fop
)
310 if (fop
->fr
) { /* any fraction? */
311 uint32 test
= (fop
->fr
>> 1) & FP_NORM
;
312 while ((fop
->fr
& FP_NORM
) == test
) { /* until norm */
313 fop
->exp
= fop
->exp
- 1;
314 fop
->fr
= (fop
->fr
<< 1);
317 else fop
->exp
= 0; /* clean 0 */
323 uint32
PackFP (struct ufp
*fop
)
325 return (fop
->fr
& FP_FR
) | /* merge frac */
326 ((fop
->exp
& FP_M_EXP
) << FP_V_EXP
) | /* and exp */
327 ((fop
->exp
< 0)? (1 << FP_V_EXPS
): 0); /* add exp sign */
330 /* Round fp number, store, generate overflow */
332 uint32
StoreFP (struct ufp
*fop
)
334 uint32 sign
, svfr
, hi
, ov
= 0;
336 NormFP (fop
); /* normalize */
337 svfr
= fop
->fr
; /* save fraction */
338 sign
= FP_GETSIGN (fop
->fr
); /* save sign */
339 fop
->fr
= (fop
->fr
+ (sign
? FP_RNDM
: FP_RNDP
)) & FP_FR
; /* round */
340 if ((fop
->fr
^ svfr
) & FP_SIGN
) { /* sign change? */
341 fop
->fr
= fop
->fr
>> 1; /* renormalize */
342 fop
->exp
= fop
->exp
+ 1;
344 else NormFP (fop
); /* check for norm */
345 if (fop
->fr
== 0) hi
= 0; /* result 0? */
346 else if (fop
->exp
< -(FP_M_EXP
+ 1)) { /* underflow? */
347 hi
= 0; /* store clean 0 */
350 else if (fop
->exp
> FP_M_EXP
) { /* overflow? */
351 hi
= 0x7FFFFFFE; /* all 1's */
354 else hi
= PackFP (fop
); /* pack mant and exp */
355 AR
= (hi
>> 16) & DMASK
;
361 /* Single-precision Fast FORTRAN Processor helpers. */
363 /* Pack mantissa and exponent and return fp value. */
365 uint32
fp_pack (OP
*result
, OP mantissa
, int32 exponent
, OPSIZE precision
)
370 fop
.fr
= ((uint32
) mantissa
.fpk
[0] << 16) | mantissa
.fpk
[1];
373 result
->fpk
[0] = (int16
) (val
>> 16);
374 result
->fpk
[1] = (int16
) val
;
378 /* Normalize, round, and pack mantissa and exponent and return fp value. */
380 uint32
fp_nrpack (OP
*result
, OP mantissa
, int32 exponent
, OPSIZE precision
)
385 fop
.fr
= ((uint32
) mantissa
.fpk
[0] << 16) | mantissa
.fpk
[1];
387 ovf
= StoreFP (&fop
);
393 /* Unpack fp number in into mantissa and exponent. */
395 uint32
fp_unpack (OP
*mantissa
, int32
*exponent
, OP packed
, OPSIZE precision
)
400 operand
= ((uint32
) packed
.fpk
[0] << 16) | packed
.fpk
[1];
401 UnpackFP (&fop
, operand
);
402 mantissa
->fpk
[0] = (uint16
) (fop
.fr
>> 16);
403 mantissa
->fpk
[1] = (uint16
) fop
.fr
;
408 #endif /* int64 support unavailable */