dc6dd0b2aa19883cccfb70973582226dc3738e3d
1 /* i7094_cpu1.c: IBM 7094 CPU complex instructions
3 Copyright (c) 2003-2006, 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.
27 #include "i7094_defs.h"
29 #define FP_HIFRAC(x) ((uint32) ((x) >> FP_N_FR) & FP_FMASK)
30 #define FP_LOFRAC(x) ((uint32) (x) & FP_FMASK)
32 #define FP_PACK38(s,e,f) (((s)? AC_S: 0) | ((t_uint64) (f)) | \
33 (((t_uint64) ((e) & FP_M_ACCH)) << FP_V_CH))
34 #define FP_PACK36(s,e,f) (((s)? SIGN: 0) | ((t_uint64) (f)) | \
35 (((t_uint64) ((e) & FP_M_CH)) << FP_V_CH))
37 extern t_uint64 AC
, MQ
, SI
, KEYS
;
39 extern uint32 SLT
, SSW
;
40 extern uint32 cpu_model
, stop_illop
;
41 extern uint32 ind_ovf
, ind_dvc
, ind_ioc
, ind_mqo
;
42 extern uint32 mode_ttrap
, mode_strap
, mode_ctrap
, mode_ftrap
;
43 extern uint32 mode_storn
, mode_multi
;
44 extern uint32 chtr_pend
, chtr_inht
, chtr_inhi
;
45 extern uint32 ch_flags
[NUM_CHAN
];
47 typedef struct { /* unpacked fp */
48 uint32 s
; /* sign: 0 +, 1 - */
49 int32 ch
; /* exponent */
50 t_uint64 fr
; /* fraction (54b) */
53 uint32
op_frnd (void);
54 t_uint64
fp_fracdiv (t_uint64 dvd
, t_uint64 dvr
, t_uint64
*rem
);
55 void fp_norm (UFP
*op
);
56 void fp_unpack (t_uint64 h
, t_uint64 l
, t_bool q_ac
, UFP
*op
);
57 uint32
fp_pack (UFP
*op
, uint32 mqs
, int32 mqch
);
59 extern t_bool
fp_trap (uint32 spill
);
60 extern t_bool
sel_trap (uint32 va
);
61 extern t_stat
ch_op_reset (uint32 ch
, t_bool ch7909
);
65 Sherman: "As the result of an addition or subtraction, if the C(AC) is
66 zero, the sign of AC is unchanged." */
68 void op_add (t_uint64 op
)
70 t_uint64 mac
= AC
& AC_MMASK
; /* get magnitudes */
71 t_uint64 mop
= op
& MMASK
;
73 AC
= AC
& AC_S
; /* isolate AC sign */
74 if ((AC
? 1: 0) ^ ((op
& SIGN
)? 1: 0)) { /* signs diff? sub */
75 if (mac
>= mop
) AC
= AC
| (mac
- mop
); /* AC >= MQ */
76 else AC
= (AC
^ AC_S
) | (mop
- mac
); /* <, sign change */
79 AC
= AC
| ((mac
+ mop
) & AC_MMASK
); /* signs same, add */
80 if ((AC
^ mac
) & AC_P
) ind_ovf
= 1; /* P change? overflow */
87 void op_mpy (t_uint64 ac
, t_uint64 sr
, uint32 sc
)
91 if (sc
== 0) return; /* sc = 0? nop */
92 sign
= ((MQ
& SIGN
)? 1: 0) ^ ((sr
& SIGN
)? 1: 0); /* result sign */
93 ac
= ac
& AC_MMASK
; /* clear AC sign */
94 sr
= sr
& MMASK
; /* mpy magnitude */
95 MQ
= MQ
& MMASK
; /* MQ magnitude */
96 if (sr
&& MQ
) { /* mpy != 0? */
97 while (sc
--) { /* for sc */
98 if (MQ
& 1) ac
= (ac
+ sr
) & AC_MMASK
; /* MQ35? AC += mpy */
99 MQ
= (MQ
>> 1) | ((ac
& 1) << 34); /* AC'MQ >> 1 */
103 else ac
= MQ
= 0; /* result = 0 */
104 if (sign
) { /* negative? */
105 ac
= ac
| AC_S
; /* insert signs */
108 AC
= ac
; /* update AC */
114 t_bool
op_div (t_uint64 sr
, uint32 sc
)
118 if (sc
== 0) return FALSE
; /* sc = 0? nop */
119 signa
= (AC
& AC_S
)? 1: 0; /* get signs */
120 signm
= (sr
& SIGN
)? 1: 0;
121 sr
= sr
& MMASK
; /* get dvr magn */
122 if ((AC
& AC_MMASK
) >= sr
) return TRUE
; /* |AC| >= |sr|? */
123 AC
= AC
& AC_MMASK
; /* AC, MQ magn */
125 while (sc
--) { /* for sc */
126 AC
= ((AC
<< 1) & AC_MMASK
) | (MQ
>> 34); /* AC'MQ << 1 */
127 MQ
= (MQ
<< 1) & MMASK
;
128 if (AC
>= sr
) { /* AC >= dvr? */
129 AC
= AC
- sr
; /* AC -= dvr */
130 MQ
= MQ
| 1; /* set quo bit */
133 if (signa
^ signm
) MQ
= MQ
| SIGN
; /* quo neg? */
134 if (signa
) AC
= AC
| AC_S
; /* rem neg? */
135 return FALSE
; /* div ok */
140 void op_als (uint32 addr
)
142 uint32 sc
= addr
& SCMASK
;
144 if ((sc
>= 35)? /* shift >= 35? */
145 ((AC
& MMASK
) != 0): /* test all bits for ovf */
146 (((AC
& MMASK
) >> (35 - sc
)) != 0)) /* test only 35-sc bits */
148 if (sc
>= 37) AC
= AC
& AC_S
; /* sc >= 37? result 0 */
149 else AC
= (AC
& AC_S
) | ((AC
<< sc
) & AC_MMASK
); /* shift, save sign */
153 void op_ars (uint32 addr
)
155 uint32 sc
= addr
& SCMASK
;
157 if (sc
>= 37) AC
= AC
& AC_S
; /* sc >= 37? result 0 */
158 else AC
= (AC
& AC_S
) | ((AC
& AC_MMASK
) >> sc
); /* shift, save sign */
162 void op_lls (uint32 addr
)
164 uint32 sc
; /* get sc */
166 AC
= AC
& AC_MMASK
; /* clear AC sign */
167 for (sc
= addr
& SCMASK
; sc
!= 0; sc
--) { /* for SC */
168 AC
= ((AC
<< 1) & AC_MMASK
) | ((MQ
>> 34) & 1); /* AC'MQ << 1 */
169 MQ
= (MQ
& SIGN
) | ((MQ
<< 1) & MMASK
); /* preserve MQ sign */
170 if (AC
& AC_P
) ind_ovf
= 1; /* if P, overflow */
172 if (MQ
& SIGN
) AC
= AC
| AC_S
; /* set ACS from MQS */
176 void op_lrs (uint32 addr
)
178 uint32 sc
= addr
& SCMASK
;
181 MQ
= MQ
& MMASK
; /* get MQ magnitude */
183 mac
= AC
& AC_MMASK
; /* get AC magnitude, */
184 AC
= AC
& AC_S
; /* sign */
185 if (sc
< 35) { /* sc [1,34]? */
186 MQ
= ((MQ
>> sc
) | (mac
<< (35 - sc
))) & MMASK
; /* MQ has AC'MQ */
187 AC
= AC
| (mac
>> sc
); /* AC has AC only */
189 else if (sc
< 37) { /* sc [35:36]? */
190 MQ
= (mac
>> (sc
- 35)) & MMASK
; /* MQ has AC only */
191 AC
= AC
| (mac
>> sc
); /* AC has <QP> */
193 else if (sc
< 72) /* sc [37:71]? */
194 MQ
= (mac
>> (sc
- 35)) & MMASK
; /* MQ has AC only */
195 else MQ
= 0; /* >72? MQ = 0 */
197 if (AC
& AC_S
) MQ
= MQ
| SIGN
; /* set MQS from ACS */
201 void op_lgl (uint32 addr
)
203 uint32 sc
; /* get sc */
205 for (sc
= addr
& SCMASK
; sc
!= 0; sc
--) { /* for SC */
206 AC
= (AC
& AC_S
) | ((AC
<< 1) & AC_MMASK
) | /* AC'MQ << 1 */
207 ((MQ
>> 35) & 1); /* preserve AC sign */
208 MQ
= (MQ
<< 1) & DMASK
;
209 if (AC
& AC_P
) ind_ovf
= 1; /* if P, overflow */
214 void op_lgr (uint32 addr
)
216 uint32 sc
= addr
& SCMASK
;
220 mac
= AC
& AC_MMASK
; /* get AC magnitude, */
221 AC
= AC
& AC_S
; /* sign */
222 if (sc
< 36) { /* sc [1,35]? */
223 MQ
= ((MQ
>> sc
) | (mac
<< (36 - sc
))) & DMASK
; /* MQ has AC'MQ */
224 AC
= AC
| (mac
>> sc
); /* AC has AC only */
226 else if (sc
== 36) { /* sc [36]? */
227 MQ
= mac
& DMASK
; /* MQ = AC<P,1:35> */
228 AC
= AC
| (mac
>> 36); /* AC = AC<Q> */
230 else if (sc
< 73) /* sc [37, 72]? */
231 MQ
= (mac
>> (sc
- 36)) & DMASK
; /* MQ has AC only */
232 else MQ
= 0; /* >72, AC,MQ = 0 */
237 /* Plus sense - undefined operations are NOPs */
239 t_stat
op_pse (uint32 addr
)
245 case 00000: /* CLM */
246 if (cpu_model
& I_9X
) AC
= AC
& AC_S
; /* 709X only */
249 case 00001: /* LBT */
250 if ((AC
& 1) != 0) PC
= (PC
+ 1) & AMASK
;
253 case 00002: /* CHS */
257 case 00003: /* SSP */
261 case 00004: /* ENK */
265 case 00005: /* IOT */
266 if (ind_ioc
) ind_ioc
= 0;
267 else PC
= (PC
+ 1) & AMASK
;
270 case 00006: /* COM */
274 case 00007: /* ETM */
275 if (cpu_model
& I_9X
) mode_ttrap
= 1; /* 709X only */
278 case 00010: /* RND */
279 if ((cpu_model
& I_9X
) && (MQ
& B1
)) /* 709X only, MQ1 set? */
280 op_add ((t_uint64
) 1); /* incr AC */
283 case 00011: /* FRN */
284 if (cpu_model
& I_9X
) { /* 709X only */
286 if (spill
) fp_trap (spill
);
290 case 00012: /* DCT */
291 if (ind_dvc
) ind_dvc
= 0;
292 else PC
= (PC
+ 1) & AMASK
;
295 case 00014: /* RCT */
296 chtr_inhi
= 1; /* 1 cycle delay */
297 chtr_inht
= 0; /* clr inhibit trap */
298 chtr_pend
= 0; /* no trap now */
301 case 00016: /* LMTM */
302 if (cpu_model
& I_94
) mode_multi
= 0; /* 709X only */
305 case 00140: /* SLF */
306 if (cpu_model
& I_9X
) SLT
= 0; /* 709X only */
309 case 00141: case 00142: case 00143: case 00144: /* SLN */
310 if (cpu_model
& I_9X
) /* 709X only */
311 SLT
= SLT
| (1u << (00144 - addr
));
314 case 00161: case 00162: case 00163: /* SWT */
315 case 00164: case 00165: case 00166:
316 if ((SSW
& (1u << (00166 - addr
))) != 0)
317 PC
= (PC
+ 1) & AMASK
;
320 case 01000: case 02000: case 03000: case 04000: /* BTT */
321 case 05000: case 06000: case 07000: case 10000:
322 if (cpu_model
& I_9X
) { /* 709X only */
323 if (sel_trap (PC
)) break; /* sel trap? */
324 ch
= GET_U_CH (addr
); /* get channel */
325 if (ch_flags
[ch
] & CHF_BOT
) /* BOT? */
326 ch_flags
[ch
] &= ~CHF_BOT
; /* clear */
327 else PC
= (PC
+ 1) & AMASK
; /* else skip */
331 case 001350: case 002350: case 003350: case 004350: /* RICx */
332 case 005350: case 006350: case 007350: case 010350:
333 ch
= GET_U_CH (addr
); /* get channel */
334 return ch_op_reset (ch
, 1);
336 case 001352: case 002352: case 003352: case 004352: /* RDCx */
337 case 005352: case 006352: case 007352: case 010352:
338 ch
= GET_U_CH (addr
); /* get channel */
339 return ch_op_reset (ch
, 0);
347 t_stat
op_mse (uint32 addr
)
353 case 00000: /* CLM */
354 if (cpu_model
& I_9X
) AC
= AC
& AC_S
; /* 709X only */
357 case 00001: /* PBT */
358 if ((AC
& AC_P
) != 0) PC
= (PC
+ 1) & AMASK
;
361 case 00002: /* EFTM */
362 if (cpu_model
& I_9X
) { /* 709X only */
364 ind_mqo
= 0; /* clears MQ ovf */
368 case 00003: /* SSM */
369 if (cpu_model
& I_9X
) AC
= AC
| AC_S
; /* 709X only */
372 case 00004: /* LFTM */
373 if (cpu_model
& I_9X
) mode_ftrap
= 0; /* 709X only */
376 case 00005: /* ESTM */
377 if (cpu_model
& I_9X
) mode_strap
= 1; /* 709X only */
380 case 00006: /* ECTM */
381 if (cpu_model
& I_9X
) mode_ctrap
= 1; /* 709X only */
384 case 00007: /* LTM */
385 if (cpu_model
& I_9X
) mode_ttrap
= 0; /* 709X only */
388 case 00010: /* LSNM */
389 if (cpu_model
& I_9X
) mode_storn
= 0; /* 709X only */
392 case 00012: /* RTT (704) */
393 if (cpu_model
& I_9X
) sel_trap (PC
); /* 709X only */
396 case 00016: /* EMTM */
400 case 00140: /* SLF */
401 if (cpu_model
& I_9X
) SLT
= 0; /* 709X only */
404 case 00141: case 00142: case 00143: case 00144: /* SLT */
405 if (cpu_model
& I_9X
) { /* 709X only */
406 t
= SLT
& (1u << (00144 - addr
));
408 if (t
!= 0) PC
= (PC
+ 1) & AMASK
;
412 case 00161: case 00162: case 00163: /* SWT */
413 case 00164: case 00165: case 00166:
414 if ((cpu_model
& I_9X
) && /* 709X only */
415 ((SSW
& (1u << (00166 - addr
))) != 0))
416 PC
= (PC
+ 1) & AMASK
;
419 case 001000: case 002000: case 003000: case 004000: /* ETT */
420 case 005000: case 006000: case 007000: case 010000:
421 if (sel_trap (PC
)) break; /* sel trap? */
422 ch
= GET_U_CH (addr
); /* get channel */
423 if (ch_flags
[ch
] & CHF_EOT
) /* EOT? */
424 ch_flags
[ch
] = ch_flags
[ch
] & ~CHF_EOT
; /* clear */
425 else PC
= (PC
+ 1) & AMASK
; /* else skip */
435 - AC<Q,P> enter into the initial exponent comparison. If either is set,
436 the numbers are always swapped. AC<P> gets OR'd into AC<S> during the
437 swap, and AC<Q,P> are cleared afterwards
438 - The early end test is actually > 077 if AC <= SR and > 100 if
439 AC > SR. However, any shift >= 54 will produce a zero fraction,
440 so the difference can be ignored */
442 uint32
op_fad (t_uint64 sr
, t_bool norm
)
447 MQ
= 0; /* clear MQ */
448 fp_unpack (AC
, 0, 1, &op1
); /* unpack AC */
449 fp_unpack (sr
, 0, 0, &op2
); /* unpack sr */
450 if (op1
.ch
> op2
.ch
) { /* AC exp > SR exp? */
451 if (AC
& AC_P
) op1
.s
= 1; /* AC P or's with S */
452 t
= op1
; /* swap operands */
455 op2
.ch
= op2
.ch
& FP_M_CH
; /* clear P,Q */
457 diff
= op2
.ch
- op1
.ch
; /* exp diff */
458 if (diff
) { /* any shift? */
459 if ((diff
< 0) || (diff
> 077)) op1
.fr
= 0; /* diff > 63? */
460 else op1
.fr
= op1
.fr
>> diff
; /* no, denormalize */
462 if (op1
.s
^ op2
.s
) { /* subtract? */
463 if (op1
.fr
>= op2
.fr
) { /* op1 > op2? */
464 op2
.fr
= op1
.fr
- op2
.fr
; /* op1 - op2 */
465 op2
.s
= op1
.s
; /* op2 sign is result */
467 else op2
.fr
= op2
.fr
- op1
.fr
; /* else op2 - op1 */
470 op2
.fr
= op2
.fr
+ op1
.fr
; /* op2 + op1 */
471 if (op2
.fr
& FP_FCRY
) { /* carry? */
472 op2
.fr
= op2
.fr
>> 1; /* renormalize */
473 op2
.ch
++; /* incr exp */
476 if (norm
) { /* normalize? */
477 if (op2
.fr
) { /* non-zero frac? */
479 mqch
= op2
.ch
- FP_N_FR
;
481 else op2
.ch
= mqch
= 0; /* else true zero */
483 else mqch
= op2
.ch
- FP_N_FR
;
484 return fp_pack (&op2
, op2
.s
, mqch
); /* pack AC, MQ */
487 /* Floating multiply */
489 uint32
op_fmp (t_uint64 sr
, t_bool norm
)
495 fp_unpack (MQ
, 0, 0, &op1
); /* unpack MQ */
496 fp_unpack (sr
, 0, 0, &op2
); /* unpack sr */
497 op1
.s
= op1
.s
^ op2
.s
; /* result sign */
498 if ((op2
.ch
== 0) && (op2
.fr
== 0)) { /* sr a normal 0? */
499 AC
= op1
.s
? AC_S
: 0; /* result is 0 */
503 f1h
= FP_HIFRAC (op1
.fr
); /* get hi fracs */
504 f2h
= FP_HIFRAC (op2
.fr
);
505 op1
.fr
= ((t_uint64
) f1h
) * ((t_uint64
) f2h
); /* f1h * f2h */
506 op1
.ch
= (op1
.ch
& FP_M_CH
) + op2
.ch
- FP_BIAS
; /* result exponent */
507 if (norm
) { /* normalize? */
508 if (!(op1
.fr
& FP_FNORM
)) { /* not normalized? */
509 op1
.fr
= op1
.fr
<< 1; /* shift frac left 1 */
510 op1
.ch
--; /* decr exp */
512 if (FP_HIFRAC (op1
.fr
)) /* hi result non-zero? */
513 mqch
= op1
.ch
- FP_N_FR
; /* set MQ exp */
514 else op1
.ch
= mqch
= 0; /* clear AC, MQ exp */
516 else mqch
= op1
.ch
- FP_N_FR
; /* set MQ exp */
517 return fp_pack (&op1
, op1
.s
, mqch
); /* pack AC, MQ */
520 /* Floating divide */
522 uint32
op_fdv (t_uint64 sr
)
529 fp_unpack (AC
, 0, 1, &op1
); /* unpack AC */
530 fp_unpack (sr
, 0, 0, &op2
); /* unpack sr */
531 quos
= op1
.s
^ op2
.s
; /* quotient sign */
532 if (op1
.fr
>= (2 * op2
.fr
)) { /* |AC| >= 2*|sr|? */
533 MQ
= quos
? SIGN
: 0; /* MQ = sign only */
534 return TRAP_F_DVC
; /* divide check */
536 if (op1
.fr
== 0) { /* |AC| == 0? */
537 MQ
= quos
? SIGN
: 0; /* MQ = sign only */
538 AC
= 0; /* AC = +0 */
541 op1
.ch
= op1
.ch
& FP_M_CH
; /* remove AC<Q,P> */
542 if (op1
.fr
>= op2
.fr
) { /* |AC| >= |sr|? */
543 op1
.fr
= op1
.fr
>> 1; /* denorm AC */
546 op1
.fr
= fp_fracdiv (op1
.fr
, op2
.fr
, &rem
); /* fraction divide */
547 op1
.fr
= op1
.fr
| (rem
<< FP_N_FR
); /* rem'quo */
548 mqch
= op1
.ch
- op2
.ch
+ FP_BIAS
; /* quotient exp */
549 op1
.ch
= op1
.ch
- FP_N_FR
; /* remainder exp */
550 spill
= fp_pack (&op1
, quos
, mqch
); /* pack up */
551 return (spill
? (spill
| TRAP_F_SGL
): 0); /* if spill, set SGL */
554 /* Double floating add
557 - AC<Q,P> enter into the initial exponent comparison. If either is set,
558 the numbers are always swapped. AC<P> gets OR'd into AC<S> during the
559 swap, and AC<Q,P> are cleared afterwards
560 - For most cases, SI ends up with the high order part of the larger number
561 - The 'early end' cases (smaller number is shifted away) must be tracked
562 exactly for SI impacts. The early end cases are:
564 (a) AC > SR, diff > 0100, and AC normalized
565 (b) AC <= SR, diff > 077, and SR normalized
567 In case (a), SI is unchanged. In case (b), SI ends up with the SR sign
568 and characteristic but the MQ (!) fraction */
570 uint32
op_dfad (t_uint64 sr
, t_uint64 sr1
, t_bool norm
)
575 fp_unpack (AC
, MQ
, 1, &op1
); /* unpack AC'MQ */
576 fp_unpack (sr
, sr1
, 0, &op2
); /* unpack sr'sr1 */
577 if (op1
.ch
> op2
.ch
) { /* AC exp > SR exp? */
578 if (((op1
.ch
- op2
.ch
) > 0100) && (AC
& B9
)) ; /* early out */
579 else SI
= FP_PACK36 (op1
.s
, op1
.ch
, FP_HIFRAC (op1
.fr
));
580 if (AC
& AC_P
) op1
.s
= 1; /* AC P or's with S */
581 t
= op1
; /* swap operands */
584 op2
.ch
= op2
.ch
& FP_M_CH
; /* clear P,Q */
586 else { /* AC <= SR */
587 if (((op2
.ch
- op1
.ch
) > 077) && (sr
& B9
)) /* early out */
588 SI
= FP_PACK36 (op2
.s
, op2
.ch
, FP_LOFRAC (MQ
));
589 else SI
= FP_PACK36 (op2
.s
, op2
.ch
, FP_HIFRAC (op2
.fr
));
591 diff
= op2
.ch
- op1
.ch
; /* exp diff */
592 if (diff
) { /* any shift? */
593 if ((diff
< 0) || (diff
> 077)) op1
.fr
= 0; /* diff > 63? */
594 else op1
.fr
= op1
.fr
>> diff
; /* no, denormalize */
596 if (op1
.s
^ op2
.s
) { /* subtract? */
597 if (op1
.fr
>= op2
.fr
) { /* op1 > op2? */
598 op2
.fr
= op1
.fr
- op2
.fr
; /* op1 - op2 */
599 op2
.s
= op1
.s
; /* op2 sign is result */
601 else op2
.fr
= op2
.fr
- op1
.fr
; /* op2 - op1 */
604 op2
.fr
= op2
.fr
+ op1
.fr
; /* op2 + op1 */
605 if (op2
.fr
& FP_FCRY
) { /* carry? */
606 op2
.fr
= op2
.fr
>> 1; /* renormalize */
607 op2
.ch
++; /* incr exp */
610 if (norm
) { /* normalize? */
611 if (op2
.fr
) { /* non-zero frac? */
613 mqch
= op2
.ch
- FP_N_FR
;
615 else op2
.ch
= mqch
= 0; /* else true zero */
617 else mqch
= op2
.ch
- FP_N_FR
;
618 return fp_pack (&op2
, op2
.s
, mqch
); /* pack AC, MQ */
621 /* Double floating multiply
623 Notes (notation is A+B' * C+D', where ' denotes 2^-27):
624 - The instruction returns 0 if A and C are both zero, because B*D is never
625 done as part of the algorithm
626 - For most cases, SI ends up with B*C, with a zero sign and exponent
627 - For the A+B' both zero 'early end' case SI ends up with A or C,
628 depending on whether the operation is normalized or not */
630 uint32
op_dfmp (t_uint64 sr
, t_uint64 sr1
, t_bool norm
)
634 uint32 f1h
, f2h
, f1l
, f2l
;
637 fp_unpack (AC
, MQ
, 1, &op1
); /* unpack AC'MQ */
638 fp_unpack (sr
, sr1
, 0, &op2
); /* unpack sr'sr1 */
639 op1
.s
= op1
.s
^ op2
.s
; /* result sign */
640 f1h
= FP_HIFRAC (op1
.fr
); /* A */
641 f1l
= FP_LOFRAC (op1
.fr
); /* B */
642 f2h
= FP_HIFRAC (op2
.fr
); /* C */
643 f2l
= FP_LOFRAC (op2
.fr
); /* D */
644 if (((op1
.ch
== 0) && (op1
.fr
== 0)) || /* AC'MQ normal 0? */
645 ((op2
.ch
== 0) && (op2
.fr
== 0)) || /* sr'sr1 normal 0? */
646 ((f1h
== 0) && (f2h
== 0))) { /* both hi frac zero? */
647 AC
= op1
.s
? AC_S
: 0; /* result is 0 */
649 SI
= sr
; /* SI has C */
652 op1
.ch
= (op1
.ch
& FP_M_CH
) + op2
.ch
- FP_BIAS
; /* result exponent */
653 if (op1
.fr
) { /* A'B != 0? */
654 op1
.fr
= ((t_uint64
) f1h
) * ((t_uint64
) f2h
); /* A * C */
655 tx
= ((t_uint64
) f1h
) * ((t_uint64
) f2l
); /* A * D */
656 op1
.fr
= op1
.fr
+ (tx
>> FP_N_FR
); /* add in hi 27b */
657 tx
= ((t_uint64
) f1l
) * ((t_uint64
) f2h
); /* B * C */
658 op1
.fr
= op1
.fr
+ (tx
>> FP_N_FR
); /* add in hi 27b */
659 SI
= tx
>> FP_N_FR
; /* SI keeps B * C */
662 if (norm
) SI
= sr
; /* early out */
663 else SI
= FP_PACK36 (op2
.s
, op2
.ch
, 0);
665 if (norm
) { /* normalize? */
666 if (!(op1
.fr
& FP_FNORM
)) { /* not normalized? */
667 op1
.fr
= op1
.fr
<< 1; /* shift frac left 1 */
668 op1
.ch
--; /* decr exp */
670 if (FP_HIFRAC (op1
.fr
)) { /* non-zero? */
671 mqch
= op1
.ch
- FP_N_FR
; /* set MQ exp */
673 else op1
.ch
= mqch
= 0; /* clear AC, MQ exp */
675 else mqch
= op1
.ch
- FP_N_FR
; /* set MQ exp */
676 return fp_pack (&op1
, op1
.s
, mqch
); /* pack AC, MQ */
679 /* Double floating divide
683 - This is a Taylor series expansion (where ' denotes >> 27):
685 (A+B') * (C+D')^-1 = (A+B') * C^-1 - (A+B') * D'* C^-2 +...
687 to two terms, which can be rewritten as terms Q1, Q2:
692 - Tracking the sign of Q2' is complicated:
694 Q1 has the sign of the quotient, s_AC ^ s_SR
695 D has the sign of the divisor, s_SR
696 R has the sign of the dividend, s_AC
697 Q1*D sign is s_AC ^ s_SR ^ s^SR = s^AC
698 Therefore, R and Q1*D have the same sign, s_AC
699 Q2' sign is s^AC ^ s_SR, which is the sign of the quotient
701 - For first divide check, SI is 0
702 - For other cases, including second divide check, SI ends up with Q1
703 - R-Q1*D is only calculated to the high 27b; using the full 54b
704 throws off the result
705 - The second divide must check for divd >= divr, otherwise an extra
706 bit of quotient would be devloped, throwing off the result
707 - A late ECO added full post-normalization; single precision divide
708 does no normalization */
710 uint32
op_dfdv (t_uint64 sr
, t_uint64 sr1
)
715 t_uint64 f1h
, f2h
, tr
, tq1
, tq1d
, trmq1d
, tq2
;
717 fp_unpack (AC
, MQ
, 1, &op1
); /* unpack AC'MQ */
718 fp_unpack (sr
, 0, 0, &op2
); /* unpack sr only */
719 ac_s
= op1
.s
; /* save AC sign */
720 op1
.s
= op1
.s
^ op2
.s
; /* sign of result */
721 f1h
= FP_HIFRAC (op1
.fr
);
722 f2h
= FP_HIFRAC (op2
.fr
);
723 if (f1h
>= (2 * f2h
)) { /* |A| >= 2*|C|? */
724 SI
= 0; /* clear SI */
725 return TRAP_F_DVC
; /* divide check */
727 if (f1h
== 0) { /* |AC| == 0? */
728 SI
= MQ
= op1
.s
? SIGN
: 0; /* MQ, SI = sign only */
729 AC
= op1
.s
? AC_S
: 0; /* AC = sign only */
732 op1
.ch
= op1
.ch
& FP_M_CH
; /* remove AC<Q,P> */
733 if (f1h
>= f2h
) { /* |A| >= |C|? */
734 op1
.fr
= op1
.fr
>> 1; /* denorm AC */
737 op1
.ch
= op1
.ch
- op2
.ch
+ FP_BIAS
; /* exp of quotient */
738 tq1
= fp_fracdiv (op1
.fr
, op2
.fr
, &tr
); /* |A+B| / |C| */
739 tr
= tr
<< FP_N_FR
; /* R << 27 */
740 tq1d
= (tq1
* ((t_uint64
) FP_LOFRAC (sr1
))) & /* Q1 * D */
741 ~((t_uint64
) FP_FMASK
); /* top 27 bits */
742 csign
= (tr
< tq1d
); /* correction sign */
743 if (csign
) trmq1d
= tq1d
- tr
; /* |R|<|Q1*D|? compl */
744 else trmq1d
= tr
- tq1d
; /* no, subtr ok */
745 SI
= FP_PACK36 (op1
.s
, op1
.ch
, tq1
); /* SI has Q1 */
746 if (trmq1d
>= (2 * op2
.fr
)) { /* |R-Q1*D| >= 2*|C|? */
747 AC
= FP_PACK38 (csign
^ ac_s
, 0, FP_HIFRAC (trmq1d
)); /* AC has R-Q1*D */
748 MQ
= (csign
^ ac_s
)? SIGN
: 0; /* MQ = sign only */
749 return TRAP_F_DVC
; /* divide check */
751 tq2
= fp_fracdiv (trmq1d
, op2
.fr
, NULL
); /* |R-Q1*D| / |C| */
752 if (trmq1d
>= op2
.fr
) tq2
&= ~((t_uint64
) 1); /* can only gen 27b quo */
753 op1
.fr
= tq1
<< FP_N_FR
; /* shift Q1 into place */
754 if (csign
) op1
.fr
= op1
.fr
- tq2
; /* sub or add Q2 */
755 else op1
.fr
= op1
.fr
+ tq2
;
756 fp_norm (&op1
); /* normalize */
757 if (op1
.fr
) mqch
= op1
.ch
- FP_N_FR
; /* non-zero? */
758 else op1
.ch
= mqch
= 0; /* clear AC, MQ exp */
759 return fp_pack (&op1
, op1
.s
, mqch
); /* pack AC, MQ */
764 uint32
op_frnd (void)
769 spill
= 0; /* no error */
770 if (MQ
& B9
) { /* MQ9 set? */
771 fp_unpack (AC
, 0, 1, &op
); /* unpack AC */
772 op
.fr
= op
.fr
+ ((t_uint64
) (1 << FP_N_FR
)); /* round up */
773 if (op
.fr
& FP_FCRY
) { /* carry out? */
774 op
.fr
= op
.fr
>> 1; /* renormalize */
775 op
.ch
++; /* incr exp */
776 if (op
.ch
== (FP_M_CH
+ 1)) /* ovf with QP = 0? */
777 spill
= TRAP_F_OVF
| TRAP_F_AC
;
779 AC
= FP_PACK38 (op
.s
, op
.ch
, FP_HIFRAC (op
.fr
)); /* pack AC */
784 /* Fraction divide - 54/27'0 yielding quotient and remainder */
786 t_uint64
fp_fracdiv (t_uint64 dvd
, t_uint64 dvr
, t_uint64
*rem
)
788 dvr
= dvr
>> FP_N_FR
;
789 if (rem
) *rem
= dvd
% dvr
;
793 /* Floating point normalize */
795 void fp_norm (UFP
*op
)
797 op
->fr
= op
->fr
& FP_DFMASK
; /* mask fraction */
798 if (op
->fr
== 0) return; /* zero? */
799 while ((op
->fr
& FP_FNORM
) == 0) { /* until norm */
800 op
->fr
= op
->fr
<< 1; /* lsh 1 */
801 op
->ch
--; /* decr exp */
806 /* Floating point unpack */
808 void fp_unpack (t_uint64 h
, t_uint64 l
, t_bool q_ac
, UFP
*op
)
810 if (q_ac
) { /* AC? */
811 op
->s
= (h
& AC_S
)? 1: 0; /* get sign */
812 op
->ch
= (uint32
) ((h
>> FP_V_CH
) & FP_M_ACCH
); /* get exp */
815 op
->s
= (h
& SIGN
)? 1: 0; /* no, mem */
816 op
->ch
= (uint32
) ((h
>> FP_V_CH
) & FP_M_CH
);
818 op
->fr
= (((t_uint64
) FP_LOFRAC (h
)) << FP_N_FR
) | /* get frac hi */
819 ((t_uint64
) FP_LOFRAC (l
)); /* get frac lo */
823 /* Floating point pack */
825 uint32
fp_pack (UFP
*op
, uint32 mqs
, int32 mqch
)
829 AC
= FP_PACK38 (op
->s
, op
->ch
, FP_HIFRAC (op
->fr
)); /* pack AC */
830 MQ
= FP_PACK36 (mqs
, mqch
, FP_LOFRAC (op
->fr
)); /* pack MQ */
831 if (op
->ch
> FP_M_CH
) spill
= TRAP_F_OVF
| TRAP_F_AC
; /* check AC exp */
832 else if (op
->ch
< 0) spill
= TRAP_F_AC
;
834 if (mqch
> FP_M_CH
) spill
|= (TRAP_F_OVF
| TRAP_F_MQ
); /* check MQ exp */
835 else if (mqch
< 0) spill
|= TRAP_F_MQ
;