I7094/i7094_cpu1.c

   1 /* i7094_cpu1.c: IBM 7094 CPU complex instructions
   2
   3    Copyright (c) 2003-2006, 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
  27 #include "i7094_defs.h"
  28
  29 #define FP_HIFRAC(x)    ((uint32) ((x) >> FP_N_FR) & FP_FMASK)
  30 #define FP_LOFRAC(x)    ((uint32) (x) & FP_FMASK)
  31
  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))
  36
  37 extern t_uint64 AC, MQ, SI, KEYS;
  38 extern uint32 PC;
  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];
  46
  47 typedef struct {                                        /* unpacked fp */
  48     uint32              s;                              /* sign: 0 +, 1 - */
  49     int32               ch;                             /* exponent */
  50     t_uint64            fr;                             /* fraction (54b) */
  51     } UFP;
  52
  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);
  58
  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);
  62
  63 /* Integer add
  64
  65    Sherman: "As the result of an addition or subtraction, if the C(AC) is
  66    zero, the sign of AC is unchanged." */
  67
  68 void op_add (t_uint64 op)
  69 {
  70 t_uint64 mac = AC & AC_MMASK;                           /* get magnitudes */
  71 t_uint64 mop = op & MMASK;
  72
  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 */
  77     }
  78 else {
  79     AC = AC | ((mac + mop) & AC_MMASK);                 /* signs same, add */
  80     if ((AC ^ mac) & AC_P) ind_ovf = 1;                 /* P change? overflow */
  81     }
  82 return;
  83 }
  84
  85 /* Multiply */
  86
  87 void op_mpy (t_uint64 ac, t_uint64 sr, uint32 sc)
  88 {
  89 uint32 sign;
  90
  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 */
 100         ac = ac >> 1;
 101         }
 102     }
 103 else ac = MQ = 0;                                       /* result = 0 */
 104 if (sign) {                                             /* negative? */
 105     ac = ac | AC_S;                                     /* insert signs */
 106     MQ = MQ | SIGN;
 107     }
 108 AC = ac;                                                /* update AC */
 109 return;
 110 }
 111
 112 /* Divide */
 113
 114 t_bool op_div (t_uint64 sr, uint32 sc)
 115 {
 116 uint32 signa, signm;
 117
 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 */
 124 MQ = MQ & MMASK;
 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 */
 131         }
 132     }
 133 if (signa ^ signm) MQ = MQ | SIGN;                      /* quo neg? */
 134 if (signa) AC = AC | AC_S;                              /* rem neg? */
 135 return FALSE;                                           /* div ok */
 136 }
 137
 138 /* Shifts */
 139
 140 void op_als (uint32 addr)
 141 {
 142 uint32 sc = addr & SCMASK;
 143
 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 */
 147     ind_ovf = 1;
 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 */
 150 return;
 151 }
 152
 153 void op_ars (uint32 addr)
 154 {
 155 uint32 sc = addr & SCMASK;
 156
 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 */
 159 return;
 160 }
 161
 162 void op_lls (uint32 addr)
 163 {
 164 uint32 sc;                                              /* get sc */
 165
 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 */
 171     }
 172 if (MQ & SIGN) AC = AC | AC_S;                          /* set ACS from MQS */
 173 return;
 174 }
 175
 176 void op_lrs (uint32 addr)
 177 {
 178 uint32 sc = addr & SCMASK;
 179 t_uint64 mac;
 180
 181 MQ = MQ & MMASK;                                        /* get MQ magnitude */
 182 if (sc != 0) {
 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 */
 188         }
 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> */
 192         }
 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 */
 196     }
 197 if (AC & AC_S) MQ = MQ | SIGN;                          /* set MQS from ACS */
 198 return;
 199 }
 200
 201 void op_lgl (uint32 addr)
 202 {
 203 uint32 sc;                                              /* get sc */
 204
 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 */
 210     }
 211 return;
 212 }
 213
 214 void op_lgr (uint32 addr)
 215 {
 216 uint32 sc = addr & SCMASK;
 217 t_uint64 mac;
 218
 219 if (sc != 0) {
 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 */
 225         }
 226     else if (sc == 36) {                                /* sc [36]? */
 227         MQ = mac & DMASK;                               /* MQ = AC<P,1:35> */
 228         AC = AC | (mac >> 36);                          /* AC = AC<Q> */
 229         }
 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 */
 233     }
 234 return;
 235 }
 236
 237 /* Plus sense - undefined operations are NOPs */
 238
 239 t_stat op_pse (uint32 addr)
 240 {
 241 uint32 ch, spill;
 242
 243 switch (addr) {
 244
 245     case 00000:                                         /* CLM */
 246         if (cpu_model & I_9X) AC = AC & AC_S;           /* 709X only */
 247         break;
 248
 249     case 00001:                                         /* LBT */
 250         if ((AC & 1) != 0) PC = (PC + 1) & AMASK;
 251         break;
 252
 253     case 00002:                                         /* CHS */
 254         AC = AC ^ AC_S;
 255         break;
 256
 257     case 00003:                                         /* SSP */
 258         AC = AC & ~AC_S;
 259         break;
 260
 261     case 00004:                                         /* ENK */
 262         MQ = KEYS;
 263         break;
 264
 265     case 00005:                                         /* IOT */
 266         if (ind_ioc) ind_ioc = 0;
 267         else PC = (PC + 1) & AMASK;
 268         break;
 269
 270     case 00006:                                         /* COM */
 271         AC = AC ^ AC_MMASK;
 272         break;
 273
 274     case 00007:                                         /* ETM */
 275         if (cpu_model & I_9X) mode_ttrap = 1;           /* 709X only */
 276         break;
 277
 278     case 00010:                                         /* RND */
 279         if ((cpu_model & I_9X) && (MQ & B1))            /* 709X only, MQ1 set? */
 280             op_add ((t_uint64) 1);                      /* incr AC */
 281         break;
 282
 283     case 00011:                                         /* FRN */
 284         if (cpu_model & I_9X) {                         /* 709X only */
 285             spill = op_frnd ();
 286             if (spill) fp_trap (spill);
 287             }
 288         break;
 289
 290     case 00012:                                         /* DCT */
 291         if (ind_dvc) ind_dvc = 0;
 292         else PC = (PC + 1) & AMASK;
 293         break;
 294
 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 */
 299         break;
 300
 301     case 00016:                                         /* LMTM */
 302         if (cpu_model & I_94) mode_multi = 0;           /* 709X only */
 303         break;
 304
 305     case 00140:                                         /* SLF */
 306         if (cpu_model & I_9X) SLT = 0;                  /* 709X only */
 307         break;
 308
 309     case 00141: case 00142: case 00143: case 00144:     /* SLN */
 310         if (cpu_model & I_9X)                           /* 709X only */
 311             SLT = SLT | (1u << (00144 - addr));
 312         break;
 313
 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;
 318         break;
 319
 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 */
 328             }
 329         break;
 330
 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);
 335
 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);
 340         }                                               /* end case */
 341
 342 return SCPE_OK;
 343 }
 344
 345 /* Minus sense */
 346
 347 t_stat op_mse (uint32 addr)
 348 {
 349 uint32 t, ch;
 350
 351 switch (addr) {
 352
 353     case 00000:                                         /* CLM */
 354         if (cpu_model & I_9X) AC = AC & AC_S;           /* 709X only */
 355         break;
 356
 357     case 00001:                                         /* PBT */
 358         if ((AC & AC_P) != 0) PC = (PC + 1) & AMASK;
 359         break;
 360
 361     case 00002:                                         /* EFTM */
 362         if (cpu_model & I_9X) {                         /* 709X only */
 363             mode_ftrap = 1;
 364             ind_mqo = 0;                                /* clears MQ ovf */
 365             }
 366         break;
 367
 368     case 00003:                                         /* SSM */
 369         if (cpu_model & I_9X) AC = AC | AC_S;           /* 709X only */
 370         break;
 371
 372     case 00004:                                         /* LFTM */
 373         if (cpu_model & I_9X) mode_ftrap = 0;           /* 709X only */
 374         break;
 375
 376     case 00005:                                         /* ESTM */
 377         if (cpu_model & I_9X) mode_strap = 1;           /* 709X only */
 378         break;
 379
 380     case 00006:                                         /* ECTM */
 381         if (cpu_model & I_9X) mode_ctrap = 1;           /* 709X only */
 382         break;
 383
 384     case 00007:                                         /* LTM */
 385         if (cpu_model & I_9X) mode_ttrap = 0;           /* 709X only */
 386         break;
 387
 388     case 00010:                                         /* LSNM */
 389         if (cpu_model & I_9X) mode_storn = 0;           /* 709X only */
 390         break;
 391
 392     case 00012:                                         /* RTT (704) */
 393         if (cpu_model & I_9X) sel_trap (PC);            /* 709X only */
 394         break;
 395
 396     case 00016:                                         /* EMTM */
 397         mode_multi = 1;
 398         break;
 399
 400     case 00140:                                         /* SLF */
 401         if (cpu_model & I_9X) SLT = 0;                  /* 709X only */
 402         break;
 403
 404     case 00141: case 00142: case 00143: case 00144:     /* SLT */
 405         if (cpu_model & I_9X) {                         /* 709X only */
 406             t = SLT & (1u << (00144 - addr));
 407             SLT = SLT & ~t;
 408             if (t != 0) PC = (PC + 1) & AMASK;
 409             }
 410         break;
 411
 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;
 417         break;
 418
 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 */
 426         break;
 427         }
 428
 429 return SCPE_OK;
 430 }
 431
 432 /* Floating add
 433
 434    Notes:
 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 */
 441
 442 uint32 op_fad (t_uint64 sr, t_bool norm)
 443 {
 444 UFP op1, op2, t;
 445 int32 mqch, diff;
 446
 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 */
 453     op1 = op2;
 454     op2 = t;
 455     op2.ch = op2.ch & FP_M_CH;                          /* clear P,Q */
 456     }
 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 */
 461     }
 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 */
 466         }
 467     else op2.fr = op2.fr - op1.fr;                      /* else op2 - op1 */
 468     }
 469 else {
 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 */
 474         }
 475     }
 476 if (norm) {                                             /* normalize? */
 477     if (op2.fr) {                                       /* non-zero frac? */
 478         fp_norm (&op2);
 479         mqch = op2.ch - FP_N_FR;
 480         }
 481     else op2.ch = mqch = 0;                             /* else true zero */
 482     }
 483 else mqch = op2.ch - FP_N_FR;
 484 return fp_pack (&op2, op2.s, mqch);                     /* pack AC, MQ */
 485 }
 486
 487 /* Floating multiply */
 488
 489 uint32 op_fmp (t_uint64 sr, t_bool norm)
 490 {
 491 UFP op1, op2;
 492 int32 mqch;
 493 uint32 f1h, f2h;
 494
 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 */
 500     MQ = op1.s? SIGN: 0;
 501     return 0;
 502     }
 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 */
 511         }
 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 */
 515     }
 516 else mqch = op1.ch - FP_N_FR;                           /* set MQ exp */
 517 return fp_pack (&op1, op1.s, mqch);                     /* pack AC, MQ */
 518 }
 519
 520 /* Floating divide */
 521
 522 uint32 op_fdv (t_uint64 sr)
 523 {
 524 UFP op1, op2;
 525 int32 mqch;
 526 uint32 spill, quos;
 527 t_uint64 rem;
 528
 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 */
 535     }
 536 if (op1.fr == 0) {                                      /* |AC| == 0? */
 537     MQ = quos? SIGN: 0;                                 /* MQ = sign only */
 538     AC = 0;                                             /* AC = +0 */
 539     return 0;                                           /* done */
 540     }
 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 */
 544     op1.ch++;
 545     }
 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 */
 552 }
 553
 554 /* Double floating add
 555
 556    Notes:
 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:
 563
 564         (a) AC > SR, diff > 0100, and AC normalized
 565         (b) AC <= SR, diff > 077, and SR normalized
 566
 567      In case (a), SI is unchanged.  In case (b), SI ends up with the SR sign
 568      and characteristic but the MQ (!) fraction */
 569
 570 uint32 op_dfad (t_uint64 sr, t_uint64 sr1, t_bool norm)
 571 {
 572 UFP op1, op2, t;
 573 int32 mqch, diff;
 574
 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 */
 582     op1 = op2;
 583     op2 = t;
 584     op2.ch = op2.ch & FP_M_CH;                          /* clear P,Q */
 585     }
 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));
 590     }
 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 */
 595     }
 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 */
 600         }
 601     else  op2.fr = op2.fr - op1.fr;                     /* op2 - op1 */
 602     }
 603 else {
 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 */
 608         }
 609     }
 610 if (norm) {                                             /* normalize? */
 611     if (op2.fr) {                                       /* non-zero frac? */
 612         fp_norm (&op2);
 613         mqch = op2.ch - FP_N_FR;
 614         }
 615     else op2.ch = mqch = 0;                             /* else true zero */
 616     }
 617 else mqch = op2.ch - FP_N_FR;
 618 return fp_pack (&op2, op2.s, mqch);                     /* pack AC, MQ */
 619 }
 620
 621 /* Double floating multiply
 622
 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 */
 629
 630 uint32 op_dfmp (t_uint64 sr, t_uint64 sr1, t_bool norm)
 631 {
 632 UFP op1, op2;
 633 int32 mqch;
 634 uint32 f1h, f2h, f1l, f2l;
 635 t_uint64 tx;
 636
 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 */
 648     MQ = op1.s? SIGN: 0;
 649     SI = sr;                                            /* SI has C */
 650     return 0;
 651     }
 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 */
 660     }
 661 else {
 662     if (norm) SI = sr;                                  /* early out */
 663     else SI = FP_PACK36 (op2.s, op2.ch, 0);
 664     }
 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 */
 669         }
 670     if (FP_HIFRAC (op1.fr)) {                           /* non-zero? */
 671         mqch = op1.ch - FP_N_FR;                        /* set MQ exp */
 672         }
 673     else op1.ch = mqch = 0;                             /* clear AC, MQ exp */
 674     }
 675 else mqch = op1.ch - FP_N_FR;                           /* set MQ exp */
 676 return fp_pack (&op1, op1.s, mqch);                     /* pack AC, MQ */
 677 }
 678
 679 /* Double floating divide
 680
 681
 682    Notes:
 683    - This is a Taylor series expansion (where ' denotes >> 27):
 684
 685         (A+B') * (C+D')^-1 = (A+B') * C^-1 - (A+B') * D'* C^-2 +...
 686
 687      to two terms, which can be rewritten as terms Q1, Q2:
 688
 689         Q1 = (A+B')/C
 690         Q2' = (R - Q1*D)'/C
 691
 692    - Tracking the sign of Q2' is complicated:
 693
 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
 700
 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 */
 709
 710 uint32 op_dfdv (t_uint64 sr, t_uint64 sr1)
 711 {
 712 UFP op1, op2;
 713 int32 mqch;
 714 uint32 csign, ac_s;
 715 t_uint64 f1h, f2h, tr, tq1, tq1d, trmq1d, tq2;
 716
 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 */
 726     }
 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 */
 730     return 0;                                           /* done */
 731     }
 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 */
 735     op1.ch++;
 736     }
 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 */
 750     }
 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 */
 760 }
 761
 762 /* Floating round */
 763
 764 uint32 op_frnd (void)
 765 {
 766 UFP op;
 767 uint32 spill;
 768
 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;
 778         }
 779     AC = FP_PACK38 (op.s, op.ch, FP_HIFRAC (op.fr));    /* pack AC */
 780     }
 781 return spill;
 782 }
 783
 784 /* Fraction divide - 54/27'0 yielding quotient and remainder */
 785
 786 t_uint64 fp_fracdiv (t_uint64 dvd, t_uint64 dvr, t_uint64 *rem)
 787 {
 788 dvr = dvr >> FP_N_FR;
 789 if (rem) *rem = dvd % dvr;
 790 return (dvd / dvr);
 791 }
 792
 793 /* Floating point normalize */
 794
 795 void fp_norm (UFP *op)
 796 {
 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 */
 802     }
 803 return;
 804 }
 805
 806 /* Floating point unpack */
 807
 808 void fp_unpack (t_uint64 h, t_uint64 l, t_bool q_ac, UFP *op)
 809 {
 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 */
 813     }
 814 else {
 815     op->s = (h & SIGN)? 1: 0;                           /* no, mem */
 816     op->ch = (uint32) ((h >> FP_V_CH) & FP_M_CH);
 817     }
 818 op->fr = (((t_uint64) FP_LOFRAC (h)) << FP_N_FR) |      /* get frac hi */
 819     ((t_uint64) FP_LOFRAC (l));                         /* get frac lo */
 820 return;
 821 }
 822
 823 /* Floating point pack */
 824
 825 uint32 fp_pack (UFP *op, uint32 mqs, int32 mqch)
 826 {
 827 uint32 spill;
 828
 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;
 833 else spill = 0;
 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;
 836 return spill;
 837 }