1 /* pdp10_mdfp.c: PDP-10 multiply/divide and floating point simulator
2
3 Copyright (c) 1993-2008, Robert M Supnik
4
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:
11
12 The above copyright notice and this permission notice shall be included in
13 all copies or substantial portions of the Software.
14
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.
21
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.
25
26 2-Apr-04 RMS Fixed bug in floating point unpack
27 Fixed bug in FIXR (found by Phil Stone, fixed by
28 Chris Smith)
29 31-Aug-01 RMS Changed int64 to t_int64 for Windoze
30 10-Aug-01 RMS Removed register in declarations
31
32 Instructions handled in this module:
33 imul integer multiply
34 idiv integer divide
35 mul multiply
36 div divide
37 dmul double precision multiply
38 ddiv double precision divide
39 fad(r) floating add (and round)
40 fsb(r) floating subtract (and round)
41 fmp(r) floating multiply (and round)
42 fdv(r) floating divide and round
43 fsc floating scale
44 fix(r) floating to fixed (and round)
45 fltr fixed to floating and round
46 dfad double precision floating add/subtract
47 dfmp double precision floating multiply
48 dfdv double precision floating divide
49
50 The PDP-10 stores double (quad) precision integers in sequential
51 AC's or memory locations. Integers are stored in 2's complement
52 form. Only the sign of the high order word matters; the signs
53 in low order words are ignored on input and set to the sign of
54 the result on output. Quad precision integers exist only in the
55 AC's as the result of a DMUL or the dividend of a DDIV.
56
57 0 00000000011111111112222222222333333
58 0 12345678901234567890123456789012345
59 +-+-----------------------------------+
60 |S| high order integer | AC(n), A
61 +-+-----------------------------------+
62 |S| low order integer | AC(n + 1), A + 1
63 +-+-----------------------------------+
64 |S| low order integer | AC(n + 2)
65 +-+-----------------------------------+
66 |S| low order integer | AC(n + 3)
67 +-+-----------------------------------+
68
69 The PDP-10 supports two floating point formats: single and double
70 precision. In both, the exponent is 8 bits, stored in excess
71 128 notation. The fraction is expected to be normalized. A
72 single precision floating point number has 27 bits of fraction;
73 a double precision number has 62 bits of fraction (the sign
74 bit of the second word is ignored and is set to zero).
75
76 In a negative floating point number, the exponent is stored in
77 one's complement form, the fraction in two's complement form.
78
79 0 00000000 011111111112222222222333333
80 0 12345678 901234567890123456789012345
81 +-+--------+---------------------------+
82 |S|exponent| high order fraction | AC(n), A
83 +-+--------+---------------------------+
84 |0| low order fraction | AC(n + 1), A + 1
85 +-+------------------------------------+
86
87 Note that treatment of the sign is different for double precision
88 integers and double precision floating point. DMOVN (implemented
89 as an inline macro) follows floating point conventions.
90
91 The original PDP-10 CPU (KA10) used a different format for double
92 precision numbers and included certain instructions to make
93 software support easier. These instructions were phased out in
94 the KL10 and KS10 and are treated as MUUO's.
95
96 The KL10 added extended precision (11-bit exponent) floating point
97 format (so-called G floating). These instructions were not
98 implemented in the KS10 and are treated as MUUO's.
99 */
100
101 #include "pdp10_defs.h"
102 #include <setjmp.h>
103
104 typedef struct { /* unpacked fp number */
105 int32 sign; /* sign */
106 int32 exp; /* exponent */
107 t_uint64 fhi; /* fraction high */
108 t_uint64 flo; /* for double prec */
109 } UFP;
110
111 #define MSK32 0xFFFFFFFF
112 #define FIT27 (DMASK - 0x07FFFFFF)
113 #define FIT32 (DMASK - MSK32)
114 #define SFRC TRUE /* frac 2's comp */
115 #define AFRC FALSE /* frac abs value */
116
117 /* In packed floating point number */
118
119 #define FP_BIAS 0200 /* exponent bias */
120 #define FP_N_FHI 27 /* # of hi frac bits */
121 #define FP_V_FHI 0 /* must be zero */
122 #define FP_M_FHI 0000777777777
123 #define FP_N_EXP 8 /* # of exp bits */
124 #define FP_V_EXP (FP_V_FHI + FP_N_FHI)
125 #define FP_M_EXP 0377
126 #define FP_V_SIGN (FP_V_EXP + FP_N_EXP) /* sign */
127 #define FP_N_FLO 35 /* # of lo frac bits */
128 #define FP_V_FLO 0 /* must be zero */
129 #define FP_M_FLO 0377777777777
130 #define GET_FPSIGN(x) ((int32) (((x) >> FP_V_SIGN) & 1))
131 #define GET_FPEXP(x) ((int32) (((x) >> FP_V_EXP) & FP_M_EXP))
132 #define GET_FPHI(x) ((x) & FP_M_FHI)
133 #define GET_FPLO(x) ((x) & FP_M_FLO)
134
135 /* In unpacked floating point number */
136
137 #define FP_N_GUARD 1 /* # of guard bits */
138 #define FP_V_UFLO FP_N_GUARD /* <35:1> */
139 #define FP_V_URNDD (FP_V_UFLO - 1) /* dp round bit */
140 #define FP_V_UFHI (FP_V_UFLO + FP_N_FLO) /* <62:36> */
141 #define FP_V_URNDS (FP_V_UFHI - 1) /* sp round bit */
142 #define FP_V_UCRY (FP_V_UFHI + FP_N_FHI) /* <63> */
143 #define FP_V_UNORM (FP_V_UCRY - 1) /* normalized bit */
144 #define FP_UFHI 0x7FFFFFF000000000
145 #define FP_UFLO 0x0000000FFFFFFFFE
146 #define FP_UFRAC 0x7FFFFFFFFFFFFFFE
147 #define FP_URNDD 0x0000000000000001
148 #define FP_URNDS 0x0000000800000000
149 #define FP_UNORM 0x4000000000000000
150 #define FP_UCRY 0x8000000000000000
151 #define FP_ONES 0xFFFFFFFFFFFFFFFF
152
153 #define UNEG(x) ((~x) + 1)
154 #define DUNEG(x) x.flo = UNEG (x.flo); x.fhi = ~x.fhi + (x.flo == 0)
155
156 extern d10 *ac_cur; /* current AC block */
157 extern int32 flags; /* flags */
158 void mul (d10 a, d10 b, d10 *rs);
159 void funpack (d10 h, d10 l, UFP *r, t_bool sgn);
160 void fnorm (UFP *r, t_int64 rnd);
161 d10 fpack (UFP *r, d10 *lo, t_bool fdvneg);
162
163 /* Integer multiply - checked against KS-10 ucode */
164
imul(d10 a,d10 b)165 d10 imul (d10 a, d10 b)
166 {
167 d10 rs[2];
168
169 if ((a == SIGN) && (b == SIGN)) { /* KS10 hack */
170 SETF (F_AOV | F_T1); /* -2**35 squared */
171 return SIGN;
172 }
173 mul (a, b, rs); /* mpy, dprec result */
174 if (rs[0] && (rs[0] != ONES)) { /* high not all sign? */
175 rs[1] = TSTS (a ^ b)? SETS (rs[1]): CLRS (rs[1]); /* set sign */
176 SETF (F_AOV | F_T1); /* overflow */
177 }
178 return rs[1];
179 }
180
181 /* Integer divide, return quotient, remainder - checked against KS10 ucode
182 The KS10 does not recognize -2^35/-1 as an error. Instead, it produces
183 2^35 (that is, -2^35) as the incorrect result.
184 */
185
idiv(d10 a,d10 b,d10 * rs)186 t_bool idiv (d10 a, d10 b, d10 *rs)
187 {
188 d10 dvd = ABS (a); /* make ops positive */
189 d10 dvr = ABS (b);
190
191 if (dvr == 0) { /* divide by 0? */
192 SETF (F_DCK | F_AOV | F_T1); /* set flags, return */
193 return FALSE;
194 }
195 rs[0] = dvd / dvr; /* get quotient */
196 rs[1] = dvd % dvr; /* get remainder */
197 if (TSTS (a ^ b)) /* sign of result */
198 rs[0] = NEG (rs[0]);
199 if (TSTS (a)) /* sign of remainder */
200 rs[1] = NEG (rs[1]);
201 return TRUE;
202 }
203
204 /* Multiply, return double precision result - checked against KS10 ucode */
205
mul(d10 s1,d10 s2,d10 * rs)206 void mul (d10 s1, d10 s2, d10 *rs)
207 {
208 t_uint64 a = ABS (s1);
209 t_uint64 b = ABS (s2);
210 t_uint64 t, u, r;
211
212 if ((a == 0) || (b == 0)) { /* operand = 0? */
213 rs[0] = rs[1] = 0; /* result 0 */
214 return;
215 }
216 if ((a & FIT32) || (b & FIT32)) { /* fit in 64b? */
217 t = a >> 18; /* no, split in half */
218 a = a & RMASK; /* "dp" multiply */
219 u = b >> 18;
220 b = b & RMASK;
221 r = (a * b) + (((a * u) + (b * t)) << 18); /* low is only 35b */
222 rs[0] = ((t * u) << 1) + (r >> 35); /* so lsh hi 1 */
223 rs[1] = r & MMASK;
224 }
225 else {
226 r = a * b; /* fits, native mpy */
227 rs[0] = r >> 35; /* split at bit 35 */
228 rs[1] = r & MMASK;
229 }
230
231 if (TSTS (s1 ^ s2)) { /* result -? */
232 MKDNEG (rs);
233 }
234 else if (TSTS (rs[0])) { /* result +, 2**70? */
235 SETF (F_AOV | F_T1); /* overflow */
236 rs[1] = SETS (rs[1]); /* consistent - */
237 }
238 return;
239 }
240
241 /* Divide, return quotient and remainder - checked against KS10 ucode
242 Note that the initial divide check catches the case -2^70/-2^35;
243 thus, the quotient can have at most 35 bits.
244 */
245
divi(int32 ac,d10 b,d10 * rs)246 t_bool divi (int32 ac, d10 b, d10 *rs)
247 {
248 int32 p1 = ADDAC (ac, 1);
249 d10 dvr = ABS (b); /* make divr positive */
250 t_int64 t;
251 int32 i;
252 d10 dvd[2];
253
254 dvd[0] = AC(ac); /* divd high */
255 dvd[1] = CLRS (AC(p1)); /* divd lo, clr sgn */
256 if (TSTS (AC(ac))) { /* make divd positive */
257 DMOVN (dvd);
258 }
259 if (dvd[0] >= dvr) { /* divide fail? */
260 SETF (F_AOV | F_DCK | F_T1); /* set flags, return */
261 return FALSE;
262 }
263 if (dvd[0] & FIT27) { /* fit in 63b? */
264 for (i = 0, rs[0] = 0; i < 35; i++) { /* 35 quotient bits */
265 dvd[0] = (dvd[0] << 1) | ((dvd[1] >> 34) & 1);
266 dvd[1] = (dvd[1] << 1) & MMASK; /* shift dividend */
267 rs[0] = rs[0] << 1; /* shift quotient */
268 if (dvd[0] >= dvr) { /* subtract work? */
269 dvd[0] = dvd[0] - dvr; /* quo bit is 1 */
270 rs[0] = rs[0] + 1;
271 }
272 }
273 rs[1] = dvd[0]; /* store remainder */
274 }
275 else {
276 t = (dvd[0] << 35) | dvd[1]; /* concatenate */
277 rs[0] = t / dvr; /* quotient */
278 rs[1] = t % dvr; /* remainder */
279 }
280 if (TSTS (AC(ac) ^ b)) /* sign of result */
281 rs[0] = NEG (rs[0]);
282 if (TSTS (AC(ac))) /* sign of remainder */
283 rs[1] = NEG (rs[1]);
284 return TRUE;
285 }
286
287 /* Double precision multiply. This is done the old fashioned way. Cross
288 product multiplies would be a lot faster but would require more code.
289 */
290
dmul(int32 ac,d10 * mpy)291 void dmul (int32 ac, d10 *mpy)
292 {
293 int32 p1 = ADDAC (ac, 1);
294 int32 p2 = ADDAC (ac, 2);
295 int32 p3 = ADDAC (ac, 3);
296 int32 i;
297 d10 mpc[2], sign;
298
299 mpc[0] = AC(ac); /* mplcnd hi */
300 mpc[1] = CLRS (AC(p1)); /* mplcnd lo, clr sgn */
301 sign = mpc[0] ^ mpy[0]; /* sign of result */
302 if (TSTS (mpc[0])) { /* get abs (mpcnd) */
303 DMOVN (mpc);
304 }
305 if (TSTS (mpy[0])) { /* get abs (mpyer) */
306 DMOVN (mpy);
307 }
308 else mpy[1] = CLRS (mpy[1]); /* clear mpy lo sign */
309 AC(ac) = AC(p1) = AC(p2) = AC(p3) = 0; /* clear AC's */
310 if (((mpy[0] | mpy[1]) == 0) || ((mpc[0] | mpc[1]) == 0))
311 return;
312 for (i = 0; i < 71; i++) { /* 71 mpyer bits */
313 if (i) { /* shift res, mpy */
314 AC(p3) = (AC(p3) >> 1) | ((AC(p2) & 1) << 34);
315 AC(p2) = (AC(p2) >> 1) | ((AC(p1) & 1) << 34);
316 AC(p1) = (AC(p1) >> 1) | ((AC(ac) & 1) << 34);
317 AC(ac) = AC(ac) >> 1;
318 mpy[1] = (mpy[1] >> 1) | ((mpy[0] & 1) << 34);
319 mpy[0] = mpy[0] >> 1;
320 }
321 if (mpy[1] & 1) { /* if mpy lo bit = 1 */
322 AC(p1) = AC(p1) + mpc[1];
323 AC(ac) = AC(ac) + mpc[0] + (TSTS (AC(p1) != 0));
324 AC(p1) = CLRS (AC(p1));
325 }
326 }
327 if (TSTS (sign)) { /* result minus? */
328 AC(p3) = (-AC(p3)) & MMASK; /* quad negate */
329 AC(p2) = (~AC(p2) + (AC(p3) == 0)) & MMASK;
330 AC(p1) = (~AC(p1) + (AC(p2) == 0)) & MMASK;
331 AC(ac) = (~AC(ac) + (AC(p1) == 0)) & DMASK;
332 }
333 else if (TSTS (AC(ac))) /* wrong sign */
334 SETF (F_AOV | F_T1);
335 if (TSTS (AC(ac))) { /* if result - */
336 AC(p1) = SETS (AC(p1)); /* make signs consistent */
337 AC(p2) = SETS (AC(p2));
338 AC(p3) = SETS (AC(p3));
339 }
340 return;
341 }
342
343 /* Double precision divide - checked against KS10 ucode */
344
ddiv(int32 ac,d10 * dvr)345 void ddiv (int32 ac, d10 *dvr)
346 {
347 int32 i, cryin;
348 d10 sign, qu[2], dvd[4];
349
350 dvd[0] = AC(ac); /* save dividend */
351 for (i = 1; i < 4; i++)
352 dvd[i] = CLRS (AC(ADDAC (ac, i)));
353 sign = AC(ac) ^ dvr[0]; /* sign of result */
354 if (TSTS (AC(ac))) { /* get abs (dividend) */
355 for (i = 3, cryin = 1; i > 0; i--) { /* negate quad */
356 dvd[i] = (~dvd[i] + cryin) & MMASK; /* comp + carry in */
357 if (dvd[i]) /* next carry in */
358 cryin = 0;
359 }
360 dvd[0] = (~dvd[0] + cryin) & DMASK;
361 }
362 if (TSTS (dvr[0])) { /* get abs (divisor) */
363 DMOVN (dvr);
364 }
365 else dvr[1] = CLRS (dvr[1]);
366 if (DCMPGE (dvd, dvr)) { /* will divide work? */
367 SETF (F_AOV | F_DCK | F_T1); /* no, set flags */
368 return;
369 }
370 qu[0] = qu[1] = 0; /* clear quotient */
371 for (i = 0; i < 70; i++) { /* 70 quotient bits */
372 dvd[0] = ((dvd[0] << 1) | ((dvd[1] >> 34) & 1)) & DMASK;;
373 dvd[1] = ((dvd[1] << 1) | ((dvd[2] >> 34) & 1)) & MMASK;
374 dvd[2] = ((dvd[2] << 1) | ((dvd[3] >> 34) & 1)) & MMASK;
375 dvd[3] = (dvd[3] << 1) & MMASK; /* shift dividend */
376 qu[0] = (qu[0] << 1) | ((qu[1] >> 34) & 1); /* shift quotient */
377 qu[1] = (qu[1] << 1) & MMASK;
378 if (DCMPGE (dvd, dvr)) { /* subtract work? */
379 dvd[0] = dvd[0] - dvr[0] - (dvd[1] < dvr[1]);
380 dvd[1] = (dvd[1] - dvr[1]) & MMASK; /* do subtract */
381 qu[1] = qu[1] + 1; /* set quotient bit */
382 }
383 }
384 if (TSTS (sign) && (qu[0] | qu[1])) {
385 MKDNEG (qu);
386 }
387 if (TSTS (AC(ac)) && (dvd[0] | dvd[1])) {
388 MKDNEG (dvd);
389 }
390 AC(ac) = qu[0]; /* quotient */
391 AC(ADDAC(ac, 1)) = qu[1];
392 AC(ADDAC(ac, 2)) = dvd[0]; /* remainder */
393 AC(ADDAC(ac, 3)) = dvd[1];
394 return;
395 }
396
397 /* Single precision floating add/subtract - checked against KS10 ucode
398 The KS10 shifts the smaller operand regardless of the exponent diff.
399 This code will not shift more than 63 places; shifts beyond that
400 cannot change the value of the smaller operand.
401
402 If the signs of the operands are the same, the result sign is the
403 same as the source sign; the sign of the result fraction is actually
404 part of the data. If the signs of the operands are different, the
405 result sign is determined by the fraction sign.
406 */
407
fad(d10 op1,d10 op2,t_bool rnd,int32 inv)408 d10 fad (d10 op1, d10 op2, t_bool rnd, int32 inv)
409 {
410 int32 ediff;
411 UFP a, b, t;
412
413 if (inv) /* subtract? -b */
414 op2 = NEG (op2);
415 if (op1 == 0) /* a = 0? result is b */
416 funpack (op2, 0, &a, AFRC);
417 else if (op2 == 0) /* b = 0? result is a */
418 funpack (op1, 0, &a, AFRC);
419 else {
420 funpack (op1, 0, &a, SFRC); /* unpack operands */
421 funpack (op2, 0, &b, SFRC); /* fracs are 2's comp */
422 ediff = a.exp - b.exp; /* get exp diff */
423 if (ediff < 0) { /* a < b? switch */
424 t = a;
425 a = b;
426 b = t;
427 ediff = -ediff;
428 }
429 if (ediff > 63) /* cap diff at 63 */
430 ediff = 63;
431 if (ediff) /* shift b (signed) */
432 b.fhi = (t_int64) b.fhi >> ediff;
433 a.fhi = a.fhi + b.fhi; /* add fractions */
434 if (a.sign ^ b.sign) { /* add or subtract? */
435 if (a.fhi & FP_UCRY) { /* subtract, frac -? */
436 a.fhi = UNEG (a.fhi); /* complement result */
437 a.sign = 1; /* result is - */
438 }
439 else a.sign = 0; /* result is + */
440 }
441 else {
442 if (a.sign) /* add, src -? comp */
443 a.fhi = UNEG (a.fhi);
444 if (a.fhi & FP_UCRY) { /* check for carry */
445 a.fhi = a.fhi >> 1; /* flo won't be used */
446 a.exp = a.exp + 1;
447 }
448 }
449 }
450 fnorm (&a, (rnd? FP_URNDS: 0)); /* normalize, round */
451 return fpack (&a, NULL, FALSE);
452 }
453
454 /* Single precision floating multiply. Because the fractions are 27b,
455 a 64b multiply can be used for the fraction multiply. The 27b
456 fractions are positioned 0'frac'0000, resulting in 00'hifrac'0..0.
457 The extra 0 is accounted for by biasing the result exponent.
458 */
459
460 #define FP_V_SPM (FP_V_UFHI - (32 - FP_N_FHI - 1))
fmp(d10 op1,d10 op2,t_bool rnd)461 d10 fmp (d10 op1, d10 op2, t_bool rnd)
462 {
463 UFP a, b;
464
465 funpack (op1, 0, &a, AFRC); /* unpack operands */
466 funpack (op2, 0, &b, AFRC); /* fracs are abs val */
467 if ((a.fhi == 0) || (b.fhi == 0)) /* either 0? */
468 return 0;
469 a.sign = a.sign ^ b.sign; /* result sign */
470 a.exp = a.exp + b.exp - FP_BIAS + 1; /* result exponent */
471 a.fhi = (a.fhi >> FP_V_SPM) * (b.fhi >> FP_V_SPM); /* high 27b of result */
472 fnorm (&a, (rnd? FP_URNDS: 0)); /* normalize, round */
473 return fpack (&a, NULL, FALSE);
474 }
475
476 /* Single precision floating divide. Because the fractions are 27b, a
477 64b divide can be used for the fraction divide. Note that 28b-29b
478 of fraction are developed; the code will do one special normalize to
479 make sure that the 28th bit is not lost. Also note the special
480 treatment of negative quotients with non-zero remainders; this
481 implements the note on p2-23 of the Processor Reference Manual.
482 */
483
fdv(d10 op1,d10 op2,d10 * rs,t_bool rnd)484 t_bool fdv (d10 op1, d10 op2, d10 *rs, t_bool rnd)
485 {
486 UFP a, b;
487 t_uint64 savhi;
488 t_bool rem = FALSE;
489
490 funpack (op1, 0, &a, AFRC); /* unpack operands */
491 funpack (op2, 0, &b, AFRC); /* fracs are abs val */
492 if (a.fhi >= 2 * b.fhi) { /* will divide work? */
493 SETF (F_AOV | F_DCK | F_FOV | F_T1);
494 return FALSE;
495 }
496 if ((savhi = a.fhi)) { /* dvd = 0? quo = 0 */
497 a.sign = a.sign ^ b.sign; /* result sign */
498 a.exp = a.exp - b.exp + FP_BIAS + 1; /* result exponent */
499 a.fhi = a.fhi / (b.fhi >> (FP_N_FHI + 1)); /* do divide */
500 if (a.sign && (savhi != (a.fhi * (b.fhi >> (FP_N_FHI + 1)))))
501 rem = TRUE; /* KL/KS hack */
502 a.fhi = a.fhi << (FP_V_UNORM - FP_N_FHI - 1); /* put quo in place */
503 if ((a.fhi & FP_UNORM) == 0) { /* normalize 1b */
504 a.fhi = a.fhi << 1; /* before masking */
505 a.exp = a.exp - 1;
506 }
507 a.fhi = a.fhi & (FP_UFHI | FP_URNDS); /* mask quo to 28b */
508 }
509 fnorm (&a, (rnd? FP_URNDS: 0)); /* normalize, round */
510 *rs = fpack (&a, NULL, rem); /* pack result */
511 return TRUE;
512 }
513
514 /* Single precision floating scale. */
515
fsc(d10 val,a10 ea)516 d10 fsc (d10 val, a10 ea)
517 {
518 int32 sc = LIT8 (ea);
519 UFP a;
520
521 if (val == 0)
522 return 0;
523 funpack (val, 0, &a, AFRC); /* unpack operand */
524 if (ea & RSIGN) /* adjust exponent */
525 a.exp = a.exp - sc;
526 else a.exp = a.exp + sc;
527 fnorm (&a, 0); /* renormalize */
528 return fpack (&a, NULL, FALSE); /* pack result */
529 }
530
531 /* Float integer operand and round */
532
fltr(d10 mb)533 d10 fltr (d10 mb)
534 {
535 UFP a;
536 d10 val = ABS (mb);
537
538 a.sign = GET_FPSIGN (mb); /* get sign */
539 a.exp = FP_BIAS + 36; /* initial exponent */
540 a.fhi = val << (FP_V_UNORM - 35); /* left justify op */
541 a.flo = 0;
542 fnorm (&a, FP_URNDS); /* normalize, round */
543 return fpack (&a, NULL, FALSE); /* pack result */
544 }
545
546 /* Fix and truncate/round floating operand */
547
fix(int32 ac,d10 mb,t_bool rnd)548 void fix (int32 ac, d10 mb, t_bool rnd)
549 {
550 int32 sc;
551 t_uint64 so;
552 UFP a;
553
554 funpack (mb, 0, &a, AFRC); /* unpack operand */
555 if (a.exp > (FP_BIAS + FP_N_FHI + FP_N_EXP))
556 SETF (F_AOV | F_T1);
557 else if (a.exp < FP_BIAS) /* < 1/2? */
558 AC(ac) = 0;
559 else {
560 sc = FP_V_UNORM - (a.exp - FP_BIAS) + 1;
561 AC(ac) = a.fhi >> sc;
562 if (rnd) {
563 so = a.fhi << (64 - sc);
564 if (so >= (0x8000000000000000 + a.sign))
565 AC(ac) = AC(ac) + 1;
566 }
567 if (a.sign)
568 AC(ac) = NEG (AC(ac));
569 }
570 return;
571 }
572
573 /* Double precision floating add/subtract
574 Since a.flo is 0, adding b.flo is just a copy - this is incorporated into
575 the denormalization step. If there's no denormalization, bflo is zero too.
576 */
577
dfad(int32 ac,d10 * rs,int32 inv)578 void dfad (int32 ac, d10 *rs, int32 inv)
579 {
580 int32 p1 = ADDAC (ac, 1);
581 int32 ediff;
582 UFP a, b, t;
583
584 if (inv) { /* subtract? -b */
585 DMOVN (rs);
586 }
587 if ((AC(ac) | AC(p1)) == 0) /* a == 0? sum = b */
588 funpack (rs[0], rs[1], &a, AFRC);
589 else if ((rs[0] | rs[1]) == 0) /* b == 0? sum = a */
590 funpack (AC(ac), AC(p1), &a, AFRC);
591 else {
592 funpack (AC(ac), AC(p1), &a, SFRC); /* unpack operands */
593 funpack (rs[0], rs[1], &b, SFRC);
594 ediff = a.exp - b.exp; /* get exp diff */
595 if (ediff < 0) { /* a < b? switch */
596 t = a;
597 a = b;
598 b = t;
599 ediff = -ediff;
600 }
601 if (ediff > 127) /* cap diff at 127 */
602 ediff = 127;
603 if (ediff > 63) { /* diff > 63? */
604 a.flo = (t_int64) b.fhi >> (ediff - 64); /* b hi to a lo */
605 b.fhi = b.sign? FP_ONES: 0; /* hi = all sign */
606 }
607 else if (ediff) { /* diff <= 63 */
608 a.flo = (b.flo >> ediff) | (b.fhi << (64 - ediff));
609 b.fhi = (t_int64) b.fhi >> ediff; /* shift b (signed) */
610 }
611 a.fhi = a.fhi + b.fhi; /* do add */
612 if (a.sign ^ b.sign) { /* add or subtract? */
613 if (a.fhi & FP_UCRY) { /* subtract, frac -? */
614 DUNEG (a); /* complement result */
615 a.sign = 1; /* result is - */
616 }
617 else a.sign = 0; /* result is + */
618 }
619 else {
620 if (a.sign) { /* add, src -? comp */
621 DUNEG (a);
622 };
623 if (a.fhi & FP_UCRY) { /* check for carry */
624 a.fhi = a.fhi >> 1; /* flo won't be used */
625 a.exp = a.exp + 1;
626 }
627 }
628 }
629 fnorm (&a, FP_URNDD); /* normalize, round */
630 AC(ac) = fpack (&a, &AC(p1), FALSE); /* pack result */
631 return;
632 }
633
634 /* Double precision floating multiply
635 The 62b fractions are multiplied, with cross products, to produce a
636 124b fraction with two leading and two trailing 0's. Because the
637 product has 2 leading 0's, instead of the normal 1, an extra
638 normalization step is needed. Accordingly, the exponent calculation
639 increments the result exponent, to compensate for normalization.
640 */
641
dfmp(int32 ac,d10 * rs)642 void dfmp (int32 ac, d10 *rs)
643 {
644 int32 p1 = ADDAC (ac, 1);
645 t_uint64 xh, xl, yh, yl, mid;
646 UFP a, b;
647
648 funpack (AC(ac), AC(p1), &a, AFRC); /* unpack operands */
649 funpack (rs[0], rs[1], &b, AFRC);
650 if ((a.fhi == 0) || (b.fhi == 0)) { /* either 0? result 0 */
651 AC(ac) = AC(p1) = 0;
652 return;
653 }
654 a.sign = a.sign ^ b.sign; /* result sign */
655 a.exp = a.exp + b.exp - FP_BIAS + 1; /* result exponent */
656 xh = a.fhi >> 32; /* split 62b fracs */
657 xl = a.fhi & MSK32; /* into 32b halves */
658 yh = b.fhi >> 32;
659 yl = b.fhi & MSK32;
660 a.fhi = xh * yh; /* hi xproduct */
661 a.flo = xl * yl; /* low xproduct */
662 mid = (xh * yl) + (yh * xl); /* fits in 64b */
663 a.flo = a.flo + (mid << 32); /* add mid lo to lo */
664 a.fhi = a.fhi + ((mid >> 32) & MSK32) + (a.flo < (mid << 32));
665 fnorm (&a, FP_URNDD); /* normalize, round */
666 AC(ac) = fpack (&a, &AC(p1), FALSE); /* pack result */
667 return;
668 }
669
670 /* Double precision floating divide
671 This algorithm develops a full 62 bits of quotient, plus one rounding
672 bit, in the low order 63b of a 64b number. To do this, we must assure
673 that the initial divide step generates a 1. If it would fail, shift
674 the dividend left and decrement the result exponent accordingly.
675 */
676
dfdv(int32 ac,d10 * rs)677 void dfdv (int32 ac, d10 *rs)
678 {
679 int32 p1 = ADDAC (ac, 1);
680 int32 i;
681 t_uint64 qu = 0;
682 UFP a, b;
683
684 funpack (AC(ac), AC(p1), &a, AFRC); /* unpack operands */
685 funpack (rs[0], rs[1], &b, AFRC);
686 if (a.fhi >= 2 * b.fhi) { /* will divide work? */
687 SETF (F_AOV | F_DCK | F_FOV | F_T1);
688 return;
689 }
690 if (a.fhi) { /* dvd = 0? quo = 0 */
691 a.sign = a.sign ^ b.sign; /* result sign */
692 a.exp = a.exp - b.exp + FP_BIAS + 1; /* result exponent */
693 if (a.fhi < b.fhi) { /* make sure initial */
694 a.fhi = a.fhi << 1; /* divide step will work */
695 a.exp = a.exp - 1;
696 }
697 for (i = 0; i < 63; i++) { /* 63b of quotient */
698 qu = qu << 1; /* shift quotient */
699 if (a.fhi >= b.fhi) { /* will div work? */
700 a.fhi = a.fhi - b.fhi; /* sub, quo = 1 */
701 qu = qu + 1;
702 }
703 a.fhi = a.fhi << 1; /* shift dividend */
704 }
705 a.fhi = qu;
706 }
707 fnorm (&a, FP_URNDD); /* normalize, round */
708 AC(ac) = fpack (&a, &AC(p1), FALSE); /* pack result */
709 return;
710 }
711
712 /* Unpack floating point operand */
713
funpack(d10 h,d10 l,UFP * r,t_bool sgn)714 void funpack (d10 h, d10 l, UFP *r, t_bool sgn)
715 {
716 d10 fphi, fplo;
717
718 r->sign = GET_FPSIGN (h);
719 r->exp = GET_FPEXP (h);
720 fphi = GET_FPHI (h);
721 fplo = GET_FPLO (l);
722 r->fhi = (fphi << FP_V_UFHI) | (fplo << FP_V_UFLO);
723 r->flo = 0;
724 if (r->sign) {
725 r->exp = r->exp ^ FP_M_EXP; /* 1s comp exp */
726 if (sgn) { /* signed frac? */
727 if (r->fhi) /* extend sign */
728 r->fhi = r->fhi | FP_UCRY;
729 else {
730 r->exp = r->exp + 1;
731 r->fhi = FP_UCRY | FP_UNORM;
732 }
733 }
734 else { /* abs frac */
735 if (r->fhi)
736 r->fhi = UNEG (r->fhi) & FP_UFRAC;
737 else {
738 r->exp = r->exp + 1;
739 r->fhi = FP_UNORM;
740 }
741 }
742 }
743 return;
744 }
745
746 /* Normalize and optionally round floating point operand */
747
fnorm(UFP * a,t_int64 rnd)748 void fnorm (UFP *a, t_int64 rnd)
749 {
750 int32 i;
751 static t_uint64 normmask[6] = {
752 0x6000000000000000, 0x7800000000000000, 0x7F80000000000000,
753 0x7FFF800000000000, 0x7FFFFFFF80000000, 0x7FFFFFFFFFFFFFFF
754 };
755 static int32 normtab[7] = { 1, 2, 4, 8, 16, 32, 63 };
756 extern a10 pager_PC;
757
758 if (a->fhi & FP_UCRY) { /* carry set? */
759 printf ("%%PDP-10 FP: carry bit set at normalization, PC = %o\n", pager_PC);
760 a->flo = (a->flo >> 1) | ((a->fhi & 1) << 63); /* try to recover */
761 a->fhi = a->fhi >> 1; /* but root cause */
762 a->exp = a->exp + 1; /* should be fixed! */
763 }
764 if ((a->fhi | a->flo) == 0) { /* if fraction = 0 */
765 a->sign = a->exp = 0; /* result is 0 */
766 return;
767 }
768 while ((a->fhi & FP_UNORM) == 0) { /* normalized? */
769 for (i = 0; i < 6; i++) {
770 if (a->fhi & normmask[i])
771 break;
772 }
773 a->fhi = (a->fhi << normtab[i]) | (a->flo >> (64 - normtab[i]));
774 a->flo = a->flo << normtab[i];
775 a->exp = a->exp - normtab[i];
776 }
777 if (rnd) { /* rounding? */
778 a->fhi = a->fhi + rnd; /* add round const */
779 if (a->fhi & FP_UCRY) { /* if carry out, */
780 a->fhi = a->fhi >> 1; /* renormalize */
781 a->exp = a->exp + 1;
782 }
783 }
784 return;
785 }
786
787 /* Pack floating point result */
788
fpack(UFP * r,d10 * lo,t_bool fdvneg)789 d10 fpack (UFP *r, d10 *lo, t_bool fdvneg)
790 {
791 d10 val[2];
792
793 if (r->exp < 0)
794 SETF (F_AOV | F_FOV | F_FXU | F_T1);
795 else if (r->exp > FP_M_EXP)
796 SETF (F_AOV | F_FOV | F_T1);
797 val[0] = (((((d10) r->exp) & FP_M_EXP) << FP_V_EXP) |
798 ((r->fhi & FP_UFHI) >> FP_V_UFHI)) & DMASK;
799 if (lo)
800 val[1] = ((r->fhi & FP_UFLO) >> FP_V_UFLO) & MMASK;
801 else val[1] = 0;
802 if (r->sign) { /* negate? */
803 if (fdvneg) { /* fdvr special? */
804 val[1] = ~val[1] & MMASK; /* 1's comp */
805 val[0] = ~val[0] & DMASK;
806 }
807 else { /* 2's comp */
808 DMOVN (val);
809 }
810 }
811 if (lo)
812 *lo = val[1];
813 return val[0];
814 }
815