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