1 /* pdp10_mdfp.c: PDP-10 multiply/divide and floating point simulator
3 Copyright (c) 1993-2005, 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 2-Apr-04 RMS Fixed bug in floating point unpack
27 Fixed bug in FIXR (found by Phil Stone, fixed by
29 31-Aug-01 RMS Changed int64 to t_int64 for Windoze
30 10-Aug-01 RMS Removed register in declarations
32 Instructions handled in this module:
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
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
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.
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 +-+-----------------------------------+
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).
76 In a negative floating point number, the exponent is stored in
77 one's complement form, the fraction in two's complement form.
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 +-+------------------------------------+
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.
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.
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.
101 #include "pdp10_defs.h"
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 */
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 */
117 /* In packed floating point number */
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)
135 /* In unpacked floating point number */
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
153 #define UNEG(x) ((~x) + 1)
154 #define DUNEG(x) x.flo = UNEG (x.flo); x.fhi = ~x.fhi + (x.flo == 0)
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
);
163 /* Integer multiply - checked against KS-10 ucode */
165 d10
imul (d10 a
, d10 b
)
169 if ((a
== SIGN
) && (b
== SIGN
)) { /* KS10 hack */
170 SETF (F_AOV
| F_T1
); /* -2**35 squared */
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 */
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.
186 t_bool
idiv (d10 a
, d10 b
, d10
*rs
)
188 d10 dvd
= ABS (a
); /* make ops positive */
191 if (dvr
== 0) { /* divide by 0? */
192 SETF (F_DCK
| F_AOV
| F_T1
); /* set flags, return */
195 rs
[0] = dvd
/ dvr
; /* get quotient */
196 rs
[1] = dvd
% dvr
; /* get remainder */
197 if (TSTS (a
^ b
)) rs
[0] = NEG (rs
[0]); /* sign of result */
198 if (TSTS (a
)) rs
[1] = NEG (rs
[1]); /* sign of remainder */
202 /* Multiply, return double precision result - checked against KS10 ucode */
204 void mul (d10 s1
, d10 s2
, d10
*rs
)
206 t_uint64 a
= ABS (s1
);
207 t_uint64 b
= ABS (s2
);
210 if ((a
== 0) || (b
== 0)) { /* operand = 0? */
211 rs
[0] = rs
[1] = 0; /* result 0 */
214 if ((a
& FIT32
) || (b
& FIT32
)) { /* fit in 64b? */
215 t
= a
>> 18; /* no, split in half */
216 a
= a
& RMASK
; /* "dp" multiply */
219 r
= (a
* b
) + (((a
* u
) + (b
* t
)) << 18); /* low is only 35b */
220 rs
[0] = ((t
* u
) << 1) + (r
>> 35); /* so lsh hi 1 */
224 r
= a
* b
; /* fits, native mpy */
225 rs
[0] = r
>> 35; /* split at bit 35 */
229 if (TSTS (s1
^ s2
)) { MKDNEG (rs
); } /* result -? */
230 else if (TSTS (rs
[0])) { /* result +, 2**70? */
231 SETF (F_AOV
| F_T1
); /* overflow */
232 rs
[1] = SETS (rs
[1]); /* consistent - */
237 /* Divide, return quotient and remainder - checked against KS10 ucode
238 Note that the initial divide check catches the case -2^70/-2^35;
239 thus, the quotient can have at most 35 bits.
242 t_bool
divi (int32 ac
, d10 b
, d10
*rs
)
244 int32 p1
= ADDAC (ac
, 1);
245 d10 dvr
= ABS (b
); /* make divr positive */
250 dvd
[0] = AC(ac
); /* divd high */
251 dvd
[1] = CLRS (AC(p1
)); /* divd lo, clr sgn */
252 if (TSTS (AC(ac
))) { DMOVN (dvd
); } /* make divd positive */
253 if (dvd
[0] >= dvr
) { /* divide fail? */
254 SETF (F_AOV
| F_DCK
| F_T1
); /* set flags, return */
257 if (dvd
[0] & FIT27
) { /* fit in 63b? */
258 for (i
= 0, rs
[0] = 0; i
< 35; i
++) { /* 35 quotient bits */
259 dvd
[0] = (dvd
[0] << 1) | ((dvd
[1] >> 34) & 1);
260 dvd
[1] = (dvd
[1] << 1) & MMASK
; /* shift dividend */
261 rs
[0] = rs
[0] << 1; /* shift quotient */
262 if (dvd
[0] >= dvr
) { /* subtract work? */
263 dvd
[0] = dvd
[0] - dvr
; /* quo bit is 1 */
267 rs
[1] = dvd
[0]; /* store remainder */
270 t
= (dvd
[0] << 35) | dvd
[1]; /* concatenate */
271 rs
[0] = t
/ dvr
; /* quotient */
272 rs
[1] = t
% dvr
; /* remainder */
274 if (TSTS (AC(ac
) ^ b
)) rs
[0] = NEG (rs
[0]); /* sign of result */
275 if (TSTS (AC(ac
))) rs
[1] = NEG (rs
[1]); /* sign of remainder */
279 /* Double precision multiply. This is done the old fashioned way. Cross
280 product multiplies would be a lot faster but would require more code.
283 void dmul (int32 ac
, d10
*mpy
)
285 int32 p1
= ADDAC (ac
, 1);
286 int32 p2
= ADDAC (ac
, 2);
287 int32 p3
= ADDAC (ac
, 3);
291 mpc
[0] = AC(ac
); /* mplcnd hi */
292 mpc
[1] = CLRS (AC(p1
)); /* mplcnd lo, clr sgn */
293 sign
= mpc
[0] ^ mpy
[0]; /* sign of result */
294 if (TSTS (mpc
[0])) { DMOVN (mpc
); } /* get abs (mpcnd) */
295 if (TSTS (mpy
[0])) { DMOVN (mpy
); } /* get abs (mpyer) */
296 else mpy
[1] = CLRS (mpy
[1]); /* clear mpy lo sign */
297 AC(ac
) = AC(p1
) = AC(p2
) = AC(p3
) = 0; /* clear AC's */
298 if (((mpy
[0] | mpy
[1]) == 0) || ((mpc
[0] | mpc
[1]) == 0)) return;
299 for (i
= 0; i
< 71; i
++) { /* 71 mpyer bits */
300 if (i
) { /* shift res, mpy */
301 AC(p3
) = (AC(p3
) >> 1) | ((AC(p2
) & 1) << 34);
302 AC(p2
) = (AC(p2
) >> 1) | ((AC(p1
) & 1) << 34);
303 AC(p1
) = (AC(p1
) >> 1) | ((AC(ac
) & 1) << 34);
304 AC(ac
) = AC(ac
) >> 1;
305 mpy
[1] = (mpy
[1] >> 1) | ((mpy
[0] & 1) << 34);
306 mpy
[0] = mpy
[0] >> 1;
308 if (mpy
[1] & 1) { /* if mpy lo bit = 1 */
309 AC(p1
) = AC(p1
) + mpc
[1];
310 AC(ac
) = AC(ac
) + mpc
[0] + (TSTS (AC(p1
) != 0));
311 AC(p1
) = CLRS (AC(p1
));
314 if (TSTS (sign
)) { /* result minus? */
315 AC(p3
) = (-AC(p3
)) & MMASK
; /* quad negate */
316 AC(p2
) = (~AC(p2
) + (AC(p3
) == 0)) & MMASK
;
317 AC(p1
) = (~AC(p1
) + (AC(p2
) == 0)) & MMASK
;
318 AC(ac
) = (~AC(ac
) + (AC(p1
) == 0)) & DMASK
;
320 else if (TSTS (AC(ac
))) SETF (F_AOV
| F_T1
); /* wrong sign */
321 if (TSTS (AC(ac
))) { /* if result - */
322 AC(p1
) = SETS (AC(p1
)); /* make signs consistent */
323 AC(p2
) = SETS (AC(p2
));
324 AC(p3
) = SETS (AC(p3
));
329 /* Double precision divide - checked against KS10 ucode */
331 void ddiv (int32 ac
, d10
*dvr
)
334 d10 sign
, qu
[2], dvd
[4];
336 dvd
[0] = AC(ac
); /* save dividend */
337 for (i
= 1; i
< 4; i
++) dvd
[i
] = CLRS (AC(ADDAC (ac
, i
)));
338 sign
= AC(ac
) ^ dvr
[0]; /* sign of result */
339 if (TSTS (AC(ac
))) { /* get abs (dividend) */
340 for (i
= 3, cryin
= 1; i
> 0; i
--) { /* negate quad */
341 dvd
[i
] = (~dvd
[i
] + cryin
) & MMASK
; /* comp + carry in */
342 if (dvd
[i
]) cryin
= 0; /* next carry in */
344 dvd
[0] = (~dvd
[0] + cryin
) & DMASK
;
346 if (TSTS (dvr
[0])) { DMOVN (dvr
); } /* get abs (divisor) */
347 else dvr
[1] = CLRS (dvr
[1]);
348 if (DCMPGE (dvd
, dvr
)) { /* will divide work? */
349 SETF (F_AOV
| F_DCK
| F_T1
); /* no, set flags */
352 qu
[0] = qu
[1] = 0; /* clear quotient */
353 for (i
= 0; i
< 70; i
++) { /* 70 quotient bits */
354 dvd
[0] = ((dvd
[0] << 1) | ((dvd
[1] >> 34) & 1)) & DMASK
;;
355 dvd
[1] = ((dvd
[1] << 1) | ((dvd
[2] >> 34) & 1)) & MMASK
;
356 dvd
[2] = ((dvd
[2] << 1) | ((dvd
[3] >> 34) & 1)) & MMASK
;
357 dvd
[3] = (dvd
[3] << 1) & MMASK
; /* shift dividend */
358 qu
[0] = (qu
[0] << 1) | ((qu
[1] >> 34) & 1); /* shift quotient */
359 qu
[1] = (qu
[1] << 1) & MMASK
;
360 if (DCMPGE (dvd
, dvr
)) { /* subtract work? */
361 dvd
[0] = dvd
[0] - dvr
[0] - (dvd
[1] < dvr
[1]);
362 dvd
[1] = (dvd
[1] - dvr
[1]) & MMASK
; /* do subtract */
363 qu
[1] = qu
[1] + 1; /* set quotient bit */
366 if (TSTS (sign
) && (qu
[0] | qu
[1])) { MKDNEG (qu
); }
367 if (TSTS (AC(ac
)) && (dvd
[0] | dvd
[1])) { MKDNEG (dvd
); }
368 AC(ac
) = qu
[0]; /* quotient */
369 AC(ADDAC(ac
, 1)) = qu
[1];
370 AC(ADDAC(ac
, 2)) = dvd
[0]; /* remainder */
371 AC(ADDAC(ac
, 3)) = dvd
[1];
375 /* Single precision floating add/subtract - checked against KS10 ucode
376 The KS10 shifts the smaller operand regardless of the exponent diff.
377 This code will not shift more than 63 places; shifts beyond that
378 cannot change the value of the smaller operand.
380 If the signs of the operands are the same, the result sign is the
381 same as the source sign; the sign of the result fraction is actually
382 part of the data. If the signs of the operands are different, the
383 result sign is determined by the fraction sign.
386 d10
fad (d10 op1
, d10 op2
, t_bool rnd
, int32 inv
)
391 if (inv
) op2
= NEG (op2
); /* subtract? -b */
392 if (op1
== 0) funpack (op2
, 0, &a
, AFRC
); /* a = 0? result is b */
393 else if (op2
== 0) funpack (op1
, 0, &a
, AFRC
); /* b = 0? result is a */
395 funpack (op1
, 0, &a
, SFRC
); /* unpack operands */
396 funpack (op2
, 0, &b
, SFRC
); /* fracs are 2's comp */
397 ediff
= a
.exp
- b
.exp
; /* get exp diff */
398 if (ediff
< 0) { /* a < b? switch */
404 if (ediff
> 63) ediff
= 63; /* cap diff at 63 */
405 if (ediff
) b
.fhi
= (t_int64
) b
.fhi
>> ediff
; /* shift b (signed) */
406 a
.fhi
= a
.fhi
+ b
.fhi
; /* add fractions */
407 if (a
.sign
^ b
.sign
) { /* add or subtract? */
408 if (a
.fhi
& FP_UCRY
) { /* subtract, frac -? */
409 a
.fhi
= UNEG (a
.fhi
); /* complement result */
410 a
.sign
= 1; /* result is - */
412 else a
.sign
= 0; /* result is + */
415 if (a
.sign
) a
.fhi
= UNEG (a
.fhi
); /* add, src -? comp */
416 if (a
.fhi
& FP_UCRY
) { /* check for carry */
417 a
.fhi
= a
.fhi
>> 1; /* flo won't be used */
422 fnorm (&a
, (rnd
? FP_URNDS
: 0)); /* normalize, round */
423 return fpack (&a
, NULL
, FALSE
);
426 /* Single precision floating multiply. Because the fractions are 27b,
427 a 64b multiply can be used for the fraction multiply. The 27b
428 fractions are positioned 0'frac'0000, resulting in 00'hifrac'0..0.
429 The extra 0 is accounted for by biasing the result exponent.
432 #define FP_V_SPM (FP_V_UFHI - (32 - FP_N_FHI - 1))
433 d10
fmp (d10 op1
, d10 op2
, t_bool rnd
)
437 funpack (op1
, 0, &a
, AFRC
); /* unpack operands */
438 funpack (op2
, 0, &b
, AFRC
); /* fracs are abs val */
439 if ((a
.fhi
== 0) || (b
.fhi
== 0)) return 0; /* either 0? */
440 a
.sign
= a
.sign
^ b
.sign
; /* result sign */
441 a
.exp
= a
.exp
+ b
.exp
- FP_BIAS
+ 1; /* result exponent */
442 a
.fhi
= (a
.fhi
>> FP_V_SPM
) * (b
.fhi
>> FP_V_SPM
); /* high 27b of result */
443 fnorm (&a
, (rnd
? FP_URNDS
: 0)); /* normalize, round */
444 return fpack (&a
, NULL
, FALSE
);
447 /* Single precision floating divide. Because the fractions are 27b, a
448 64b divide can be used for the fraction divide. Note that 28b-29b
449 of fraction are developed; the code will do one special normalize to
450 make sure that the 28th bit is not lost. Also note the special
451 treatment of negative quotients with non-zero remainders; this
452 implements the note on p2-23 of the Processor Reference Manual.
455 t_bool
fdv (d10 op1
, d10 op2
, d10
*rs
, t_bool rnd
)
461 funpack (op1
, 0, &a
, AFRC
); /* unpack operands */
462 funpack (op2
, 0, &b
, AFRC
); /* fracs are abs val */
463 if (a
.fhi
>= 2 * b
.fhi
) { /* will divide work? */
464 SETF (F_AOV
| F_DCK
| F_FOV
| F_T1
);
467 if (savhi
= a
.fhi
) { /* dvd = 0? quo = 0 */
468 a
.sign
= a
.sign
^ b
.sign
; /* result sign */
469 a
.exp
= a
.exp
- b
.exp
+ FP_BIAS
+ 1; /* result exponent */
470 a
.fhi
= a
.fhi
/ (b
.fhi
>> (FP_N_FHI
+ 1)); /* do divide */
471 if (a
.sign
&& (savhi
!= (a
.fhi
* (b
.fhi
>> (FP_N_FHI
+ 1)))))
472 rem
= TRUE
; /* KL/KS hack */
473 a
.fhi
= a
.fhi
<< (FP_V_UNORM
- FP_N_FHI
- 1); /* put quo in place */
474 if ((a
.fhi
& FP_UNORM
) == 0) { /* normalize 1b */
475 a
.fhi
= a
.fhi
<< 1; /* before masking */
478 a
.fhi
= a
.fhi
& (FP_UFHI
| FP_URNDS
); /* mask quo to 28b */
480 fnorm (&a
, (rnd
? FP_URNDS
: 0)); /* normalize, round */
481 *rs
= fpack (&a
, NULL
, rem
); /* pack result */
485 /* Single precision floating scale. */
487 d10
fsc (d10 val
, a10 ea
)
489 int32 sc
= LIT8 (ea
);
492 if (val
== 0) return 0;
493 funpack (val
, 0, &a
, AFRC
); /* unpack operand */
494 if (ea
& RSIGN
) a
.exp
= a
.exp
- sc
; /* adjust exponent */
495 else a
.exp
= a
.exp
+ sc
;
496 fnorm (&a
, 0); /* renormalize */
497 return fpack (&a
, NULL
, FALSE
); /* pack result */
500 /* Float integer operand and round */
507 a
.sign
= GET_FPSIGN (mb
); /* get sign */
508 a
.exp
= FP_BIAS
+ 36; /* initial exponent */
509 a
.fhi
= val
<< (FP_V_UNORM
- 35); /* left justify op */
511 fnorm (&a
, FP_URNDS
); /* normalize, round */
512 return fpack (&a
, NULL
, FALSE
); /* pack result */
515 /* Fix and truncate/round floating operand */
517 void fix (int32 ac
, d10 mb
, t_bool rnd
)
523 funpack (mb
, 0, &a
, AFRC
); /* unpack operand */
524 if (a
.exp
> (FP_BIAS
+ FP_N_FHI
+ FP_N_EXP
)) SETF (F_AOV
| F_T1
);
525 else if (a
.exp
< FP_BIAS
) AC(ac
) = 0; /* < 1/2? */
527 sc
= FP_V_UNORM
- (a
.exp
- FP_BIAS
) + 1;
528 AC(ac
) = a
.fhi
>> sc
;
530 so
= a
.fhi
<< (64 - sc
);
531 if (so
>= (0x8000000000000000 + a
.sign
)) AC(ac
) = AC(ac
) + 1;
533 if (a
.sign
) AC(ac
) = NEG (AC(ac
));
538 /* Double precision floating add/subtract
539 Since a.flo is 0, adding b.flo is just a copy - this is incorporated into
540 the denormalization step. If there's no denormalization, bflo is zero too.
543 void dfad (int32 ac
, d10
*rs
, int32 inv
)
545 int32 p1
= ADDAC (ac
, 1);
549 if (inv
) { DMOVN (rs
); } /* subtract? -b */
550 if ((AC(ac
) | AC(p1
)) == 0) funpack (rs
[0], rs
[1], &a
, AFRC
);
551 /* a == 0? sum = b */
552 else if ((rs
[0] | rs
[1]) == 0) funpack (AC(ac
), AC(p1
), &a
, AFRC
);
553 /* b == 0? sum = a */
555 funpack (AC(ac
), AC(p1
), &a
, SFRC
); /* unpack operands */
556 funpack (rs
[0], rs
[1], &b
, SFRC
);
557 ediff
= a
.exp
- b
.exp
; /* get exp diff */
558 if (ediff
< 0) { /* a < b? switch */
564 if (ediff
> 127) ediff
= 127; /* cap diff at 127 */
565 if (ediff
> 63) { /* diff > 63? */
566 a
.flo
= (t_int64
) b
.fhi
>> (ediff
- 64); /* b hi to a lo */
567 b
.fhi
= b
.sign
? FP_ONES
: 0; /* hi = all sign */
569 else if (ediff
) { /* diff <= 63 */
570 a
.flo
= (b
.flo
>> ediff
) | (b
.fhi
<< (64 - ediff
));
571 b
.fhi
= (t_int64
) b
.fhi
>> ediff
; /* shift b (signed) */
573 a
.fhi
= a
.fhi
+ b
.fhi
; /* do add */
574 if (a
.sign
^ b
.sign
) { /* add or subtract? */
575 if (a
.fhi
& FP_UCRY
) { /* subtract, frac -? */
576 DUNEG (a
); /* complement result */
577 a
.sign
= 1; /* result is - */
579 else a
.sign
= 0; /* result is + */
582 if (a
.sign
) { DUNEG (a
); }; /* add, src -? comp */
583 if (a
.fhi
& FP_UCRY
) { /* check for carry */
584 a
.fhi
= a
.fhi
>> 1; /* flo won't be used */
589 fnorm (&a
, FP_URNDD
); /* normalize, round */
590 AC(ac
) = fpack (&a
, &AC(p1
), FALSE
); /* pack result */
594 /* Double precision floating multiply
595 The 62b fractions are multiplied, with cross products, to produce a
596 124b fraction with two leading and two trailing 0's. Because the
597 product has 2 leading 0's, instead of the normal 1, an extra
598 normalization step is needed. Accordingly, the exponent calculation
599 increments the result exponent, to compensate for normalization.
602 void dfmp (int32 ac
, d10
*rs
)
604 int32 p1
= ADDAC (ac
, 1);
605 t_uint64 xh
, xl
, yh
, yl
, mid
;
608 funpack (AC(ac
), AC(p1
), &a
, AFRC
); /* unpack operands */
609 funpack (rs
[0], rs
[1], &b
, AFRC
);
610 if ((a
.fhi
== 0) || (b
.fhi
== 0)) { /* either 0? result 0 */
614 a
.sign
= a
.sign
^ b
.sign
; /* result sign */
615 a
.exp
= a
.exp
+ b
.exp
- FP_BIAS
+ 1; /* result exponent */
616 xh
= a
.fhi
>> 32; /* split 62b fracs */
617 xl
= a
.fhi
& MSK32
; /* into 32b halves */
620 a
.fhi
= xh
* yh
; /* hi xproduct */
621 a
.flo
= xl
* yl
; /* low xproduct */
622 mid
= (xh
* yl
) + (yh
* xl
); /* fits in 64b */
623 a
.flo
= a
.flo
+ (mid
<< 32); /* add mid lo to lo */
624 a
.fhi
= a
.fhi
+ ((mid
>> 32) & MSK32
) + (a
.flo
< (mid
<< 32));
625 fnorm (&a
, FP_URNDD
); /* normalize, round */
626 AC(ac
) = fpack (&a
, &AC(p1
), FALSE
); /* pack result */
630 /* Double precision floating divide
631 This algorithm develops a full 62 bits of quotient, plus one rounding
632 bit, in the low order 63b of a 64b number. To do this, we must assure
633 that the initial divide step generates a 1. If it would fail, shift
634 the dividend left and decrement the result exponent accordingly.
637 void dfdv (int32 ac
, d10
*rs
)
639 int32 p1
= ADDAC (ac
, 1);
644 funpack (AC(ac
), AC(p1
), &a
, AFRC
); /* unpack operands */
645 funpack (rs
[0], rs
[1], &b
, AFRC
);
646 if (a
.fhi
>= 2 * b
.fhi
) { /* will divide work? */
647 SETF (F_AOV
| F_DCK
| F_FOV
| F_T1
);
650 if (a
.fhi
) { /* dvd = 0? quo = 0 */
651 a
.sign
= a
.sign
^ b
.sign
; /* result sign */
652 a
.exp
= a
.exp
- b
.exp
+ FP_BIAS
+ 1; /* result exponent */
653 if (a
.fhi
< b
.fhi
) { /* make sure initial */
654 a
.fhi
= a
.fhi
<< 1; /* divide step will work */
657 for (i
= 0; i
< 63; i
++) { /* 63b of quotient */
658 qu
= qu
<< 1; /* shift quotient */
659 if (a
.fhi
>= b
.fhi
) { /* will div work? */
660 a
.fhi
= a
.fhi
- b
.fhi
; /* sub, quo = 1 */
663 a
.fhi
= a
.fhi
<< 1; /* shift dividend */
667 fnorm (&a
, FP_URNDD
); /* normalize, round */
668 AC(ac
) = fpack (&a
, &AC(p1
), FALSE
); /* pack result */
672 /* Unpack floating point operand */
674 void funpack (d10 h
, d10 l
, UFP
*r
, t_bool sgn
)
678 r
->sign
= GET_FPSIGN (h
);
679 r
->exp
= GET_FPEXP (h
);
682 r
->fhi
= (fphi
<< FP_V_UFHI
) | (fplo
<< FP_V_UFLO
);
685 r
->exp
= r
->exp
^ FP_M_EXP
; /* 1s comp exp */
686 if (sgn
) { /* signed frac? */
687 if (r
->fhi
) r
->fhi
= r
->fhi
| FP_UCRY
; /* extend sign */
690 r
->fhi
= FP_UCRY
| FP_UNORM
;
693 else { /* abs frac */
694 if (r
->fhi
) r
->fhi
= UNEG (r
->fhi
) & FP_UFRAC
;
704 /* Normalize and optionally round floating point operand */
706 void fnorm (UFP
*a
, t_int64 rnd
)
709 static t_uint64 normmask
[6] = {
710 0x6000000000000000, 0x7800000000000000, 0x7F80000000000000,
711 0x7FFF800000000000, 0x7FFFFFFF80000000, 0x7FFFFFFFFFFFFFFF
713 static int32 normtab
[7] = { 1, 2, 4, 8, 16, 32, 63 };
716 if (a
->fhi
& FP_UCRY
) { /* carry set? */
717 printf ("%%PDP-10 FP: carry bit set at normalization, PC = %o\n", pager_PC
);
718 a
->flo
= (a
->flo
>> 1) | ((a
->fhi
& 1) << 63); /* try to recover */
719 a
->fhi
= a
->fhi
>> 1; /* but root cause */
720 a
->exp
= a
->exp
+ 1; /* should be fixed! */
722 if ((a
->fhi
| a
->flo
) == 0) { /* if fraction = 0 */
723 a
->sign
= a
->exp
= 0; /* result is 0 */
726 while ((a
->fhi
& FP_UNORM
) == 0) { /* normalized? */
727 for (i
= 0; i
< 6; i
++) {
728 if (a
->fhi
& normmask
[i
]) break;
730 a
->fhi
= (a
->fhi
<< normtab
[i
]) | (a
->flo
>> (64 - normtab
[i
]));
731 a
->flo
= a
->flo
<< normtab
[i
];
732 a
->exp
= a
->exp
- normtab
[i
];
734 if (rnd
) { /* rounding? */
735 a
->fhi
= a
->fhi
+ rnd
; /* add round const */
736 if (a
->fhi
& FP_UCRY
) { /* if carry out, */
737 a
->fhi
= a
->fhi
>> 1; /* renormalize */
744 /* Pack floating point result */
746 d10
fpack (UFP
*r
, d10
*lo
, t_bool fdvneg
)
750 if (r
->exp
< 0) SETF (F_AOV
| F_FOV
| F_FXU
| F_T1
);
751 else if (r
->exp
> FP_M_EXP
) SETF (F_AOV
| F_FOV
| F_T1
);
752 val
[0] = (((((d10
) r
->exp
) & FP_M_EXP
) << FP_V_EXP
) |
753 ((r
->fhi
& FP_UFHI
) >> FP_V_UFHI
)) & DMASK
;
754 if (lo
) val
[1] = ((r
->fhi
& FP_UFLO
) >> FP_V_UFLO
) & MMASK
;
756 if (r
->sign
) { /* negate? */
757 if (fdvneg
) { /* fdvr special? */
758 val
[1] = ~val
[1] & MMASK
; /* 1's comp */
759 val
[0] = ~val
[0] & DMASK
;
761 else { DMOVN (val
); } /* 2's comp */
763 if (lo
) *lo
= val
[1];