src/dps8/dps8

/* */
This source file includes following definitions.
ldexpl
exp2l
EAQToIEEElongdouble
EAQToIEEEdouble
IEEElongdoubleToFloat72
MYfrexpl
IEEElongdoubleToEAQ
isHex
ufa
fno
fno_ext
fneg
ufm
fdvX
fdv
fdi
frd
fstr
fcmp
fcmg
dufa
dufm
dfdvX
dfdv
dfdi
dvf
dfrd
dfstr
dfcmp
dfcmg
   1 /*
   2  * vim: filetype=c:tabstop=4:ai:expandtab
   3  * SPDX-License-Identifier: ICU
   4  * SPDX-License-Identifier: Multics
   5  * scspell-id: 9ca9a790-f62e-11ec-a677-80ee73e9b8e7
   6  *
   7  * ---------------------------------------------------------------------------
   8  *
   9  * Copyright (c) 2007-2013 Michael Mondy
  10  * Copyright (c) 2012-2016 Harry Reed
  11  * Copyright (c) 2013-2016 Charles Anthony
  12  * Copyright (c) 2016 Michal Tomek
  13  * Copyright (c) 2021-2024 The DPS8M Development Team
  14  *
  15  * This software is made available under the terms of the ICU License.
  16  * See the LICENSE.md file at the top-level directory of this distribution.
  17  *
  18  * ---------------------------------------------------------------------------
  19  *
  20  * This source file may contain code comments that adapt, include, and/or
  21  * incorporate Multics program code and/or documentation distributed under
  22  * the Multics License.  In the event of any discrepancy between code
  23  * comments herein and the original Multics materials, the original Multics
  24  * materials should be considered authoritative unless otherwise noted.
  25  * For more details and historical background, see the LICENSE.md file at
  26  * the top-level directory of this distribution.
  27  *
  28  * ---------------------------------------------------------------------------
  29  */
  30 
  31 //-V::536
  32 
  33 #include <stdio.h>
  34 #include <math.h>
  35 
  36 #include "dps8.h"
  37 #include "dps8_sys.h"
  38 #include "dps8_iom.h"
  39 #include "dps8_cable.h"
  40 #include "dps8_cpu.h"
  41 #include "dps8_faults.h"
  42 #include "dps8_scu.h"
  43 #include "dps8_ins.h"
  44 #include "dps8_math.h"
  45 #include "dps8_utils.h"
  46 
  47 #define DBG_CTR cpu.cycleCnt
  48 
  49 #if defined(__CYGWIN__)
  50 long double ldexpl(long double x, int n) {
     /*  */
  51        return __builtin_ldexpl(x, n);
  52 }
  53 
  54 long double exp2l (long double e) {
     /*  */
  55        return __builtin_exp2l(e);
  56 }
  57 #endif
  58 
  59 #if defined(NEED_128)
  60 # define ISZERO_128(m) iszero_128 (m)
  61 # define ISEQ_128(a,b) iseq_128 (a, b)
  62 #else
  63 # define ISZERO_128(m) ((m) == 0)
  64 # define ISEQ_128(a,b) ((a) == (b))
  65 #endif
  66 
  67 #if defined(__NVCOMPILER) || defined(__NVCOMPILER_LLVM__) || defined(__PGI) || defined(__PGLLVM__)
  68 # pragma global novector
  69 #endif
  70 
  71 /*
  72  * convert floating point quantity in C(EAQ) to a IEEE long double ...
  73  */
  74 long double EAQToIEEElongdouble(cpu_state_t * cpup)
     /*  */
  75 {
  76     // mantissa
  77     word72 Mant = convert_to_word72 (cpu.rA, cpu.rQ);
  78 
  79     if (ISZERO_128 (Mant))
  80         return (long double) 0;
  81 #if defined(NEED_128)
  82 
  83     bool S = isnonzero_128 (and_128 (Mant, SIGN72)); // sign of mantissa
  84     if (S)
  85         Mant = and_128 (negate_128 (Mant), MASK71); // 71 bits (not 72!)
  86 #else
  87     bool S = Mant & SIGN72; // sign of mantissa
  88     if (S)
  89         Mant = (-Mant) & MASK71; // 71 bits (not 72!)
  90 #endif
  91 
  92     long double m = 0;  // mantissa value;
  93     int e = SIGNEXT8_int (cpu . rE & MASK8); // make signed
  94 
  95     if (S && ISZERO_128 (Mant))// denormalized -1.0*2^e
  96         return -exp2l(e);
  97 
  98     long double v = 0.5;
  99     for(int n = 70 ; n >= 0 ; n -= 1) // this also normalizes the mantissa
 100     {
 101 #if defined(NEED_128)
 102         if (isnonzero_128 (and_128 (Mant, lshift_128 (construct_128 (0, 1), (unsigned int) n))))
 103         {
 104             m += v;
 105         }
 106 #else
 107         if (Mant & ((word72)1 << n))
 108         {
 109             m += v;
 110         }
 111 #endif
 112         v /= 2.0;
 113     }
 114 
 115     /*if (m == 0 && e == -128)    // special case - normalized 0
 116         return 0;
 117     if (m == 0)
 118         return (S ? -1 : 1) * exp2l(e); */
 119 
 120     return (S ? -1 : 1) * ldexpl(m, e);
 121 }
 122 
 123 #if defined(__MINGW32__) || defined(__MINGW64__)
 124 
 125 // MINGW doesn't have long double support, convert to IEEE double instead
 126 double EAQToIEEEdouble(cpu_state_t * cpup)
     /*  */
 127 {
 128     // mantissa
 129     word72 Mant = convert_to_word72 (cpu.rA, cpu.rQ);
 130 
 131     if (ISZERO_128 (Mant))
 132         return 0;
 133 
 134 # if defined(NEED_128)
 135     bool S = isnonzero_128 (and_128 (Mant, SIGN72)); // sign of mantissa
 136     if (S)
 137         Mant = and_128 (negate_128 (Mant), MASK71); // 71 bits (not 72!)
 138 # else
 139     bool S = Mant & SIGN72; // sign of mantissa
 140     if (S)
 141         Mant = (-Mant) & MASK71; // 71 bits (not 72!)
 142 # endif
 143 
 144     double m = 0;  // mantissa value
 145     int e = SIGNEXT8_int (cpu . rE & MASK8); // make signed
 146 
 147     if (S && ISZERO_128 (Mant))// denormalized -1.0*2^e
 148         return -exp2(e);
 149 
 150     double v = 0.5;
 151     for(int n = 70 ; n >= 0 ; n -= 1) // this also normalizes the mantissa
 152     {
 153 # if defined(NEED_128)
 154         if (isnonzero_128 (and_128 (Mant, lshift_128 (construct_128 (0, 1), (unsigned int) n))))
 155         {
 156             m += v;
 157         }
 158 # else
 159         if (Mant & ((word72)1 << n))
 160         {
 161             m += v;
 162         }
 163 # endif
 164         v /= 2.0;
 165     }
 166 
 167     return (S ? -1 : 1) * ldexp(m, e);
 168 }
 169 #endif
 170 
 171 #if !defined(QUIET_UNUSED)
 172 /*
 173  * return normalized dps8 representation of IEEE double f0 ...
 174  */
 175 float72 IEEElongdoubleToFloat72(long double f0)
     /*  */
 176 {
 177     if (f0 == 0)
 178         return (float72)((float72)0400000000000LL << 36);
 179 
 180     bool sign = f0 < 0 ? true : false;
 181     long double f = fabsl(f0);
 182 
 183     int exp;
 184     long double mant = frexpl(f, &exp);
 185     //sim_printf (sign=%d f0=%Lf mant=%Lf exp=%d\n", sign, f0, mant, exp);
 186 
 187     word72 result = 0;
 188 
 189     // now let's examine the mantissa and assign bits as necessary...
 190 
 191     if (sign && mant == 0.5)
 192     {
 193         //result = bitfieldInsert72(result, 1, 63, 1);
 194         putbits72 (& result, 71-62, 1, 1);
 195         exp -= 1;
 196         mant -= 0.5;
 197     }
 198 
 199     long double bitval = 0.5;    ///< start at 1/2 and go down .....
 200     for(int n = 62 ; n >= 0 && mant > 0; n -= 1)
 201     {
 202         if (mant >= bitval)
 203         {
 204             //result = bitfieldInsert72(result, 1, n, 1);
 205             putbits72 (& result, 71-n, 1, 1);
 206             mant -= bitval;
 207             //sim_printf ("Inserting a bit @ %d %012"PRIo64" %012"PRIo64"\n",
 208             //            n, (word36)((result >> 36) & DMASK), (word36)(result & DMASK));
 209         }
 210         bitval /= 2.0;
 211     }
 212     //sim_printf ("n=%d mant=%f\n", n, mant);
 213     //sim_printf ("result=%012"PRIo64" %012"PRIo64"\n",
 214     //            (word36)((result >> 36) & DMASK), (word36)(result & DMASK));
 215 
 216     // if f is < 0 then take 2-comp of result ...
 217     if (sign)
 218     {
 219         result = -result & (((word72)1 << 64) - 1);
 220         //sim_printf ("-result=%012"PRIo64" %012"PRIo64"\n",
 221         //            (word36)((result >> 36) & DMASK), (word36)(result & DMASK));
 222     }
 223     // insert exponent ...
 224     int e = (int)exp;
 225     //result = bitfieldInsert72(result, e & 0377, 64, 8);    ///< & 0777777777777LL;
 226     putbits72 (& result, 71-64, 8, e & 0377);
 227 
 228     // XXX TODO test for exp under/overflow ...
 229 
 230     return result;
 231 }
 232 #endif
 233 
 234 #if !defined(QUIET_UNUSED)
 235 static long double MYfrexpl(long double x, int *exp)
     /*  */
 236 {
 237     long double exponents[20], *next;
 238     int exponent, bit;
 239 
 240     /* Check for zero, nan and infinity. */
 241     if (x != x || x + x == x )
 242     {
 243         *exp = 0;
 244         return x;
 245     }
 246 
 247     if (x < 0)
 248         return -MYfrexpl(-x, exp);
 249 
 250     exponent = 0;
 251     if (x > 1.0)
 252     {
 253         for (next = exponents, exponents[0] = 2.0L, bit = 1;
 254              *next <= x + x;
 255              bit <<= 1, next[1] = next[0] * next[0], next++);
 256 
 257         for (; next >= exponents; bit >>= 1, next--)
 258             if (x + x >= *next)
 259             {
 260                 x /= *next;
 261                 exponent |= bit;
 262             }
 263 
 264     }
 265 
 266     else if (x < 0.5)
 267     {
 268         for (next = exponents, exponents[0] = 0.5L, bit = 1;
 269              *next > x;
 270              bit <<= 1, next[1] = next[0] * next[0], next++);
 271 
 272         for (; next >= exponents; bit >>= 1, next--)
 273             if (x < *next)
 274             {
 275                 x /= *next;
 276                 exponent |= bit;
 277             }
 278 
 279         exponent = -exponent;
 280     }
 281 
 282     *exp = exponent;
 283     return x;
 284 }
 285 #endif
 286 
 287 #if !defined(QUIET_UNUSED)
 288 /*
 289  * Store EAQ with normalized dps8 representation of IEEE double f0 ...
 290  */
 291 void IEEElongdoubleToEAQ(long double f0)
     /*  */
 292 {
 293     if (f0 == 0)
 294     {
 295         cpu . rA = 0;
 296 # if defined(TESTING)
 297         HDBGRegAW ("IEEEld2EAQ");
 298 # endif
 299         cpu . rQ = 0;
 300         cpu . rE = 0200U; /*-128*/
 301         return;
 302     }
 303 
 304     bool sign = f0 < 0 ? true : false;
 305     long double f = fabsl(f0);
 306     int exp;
 307     long double mant = MYfrexpl(f, &exp);
 308     word72 result = 0;
 309 
 310     // now let's examine the mantissa and assign bits as necessary...
 311     if (sign && mant == 0.5)
 312     {
 313         //result = bitfieldInsert72(result, 1, 63, 1);
 314         result = putbits72 (& result, 71-63, 1, 1);
 315         exp -= 1;
 316         mant -= 0.5;
 317     }
 318 
 319     long double bitval = 0.5;    ///< start at 1/2 and go down .....
 320     for(int n = 70 ; n >= 0 && mant > 0; n -= 1)
 321     {
 322         if (mant >= bitval)
 323         {
 324             //result = bitfieldInsert72(result, 1, n, 1);
 325             putbits72 (& result 71-n, 1, 1);
 326             mant -= bitval;
 327             //sim_printf ("Inserting a bit @ %d %012"PRIo64" %012"PRIo64"\n",
 328             //            n, (word36)((result >> 36) & DMASK), (word36)(result & DMASK));
 329         }
 330         bitval /= 2.0;
 331     }
 332 
 333     // if f is < 0 then take 2-comp of result ...
 334     if (sign)
 335         result = -result & (((word72)1 << 72) - 1);
 336 
 337     cpu . rE = exp & MASK8;
 338     cpu . rA = (result >> 36) & MASK36;
 339 # if defined(TESTING)
 340     HDBGRegAW ("IEEEld2EAQ");
 341 # endif
 342     cpu . rQ = result & MASK36;
 343 }
 344 #endif
 345 
 346 
 347 
 348 
 349 
 350 
 351 
 352 
 353 
 354 
 355 
 356 
 357 
 358 
 359 
 360 
 361 
 362 
 363 
 364 
 365 
 366 
 367 
 368 
 369 
 370 
 371 
 372 
 373 
 374 
 375 
 376 
 377 
 378 
 379 
 380 
 381 
 382 
 383 
 384 
 385 
 386 
 387 
 388 
 389 
 390 
 391 
 392 
 393 
 394 
 395 
 396 
 397 
 398 
 399 
 400 
 401 
 402 
 403 
 404 
 405 
 406 
 407 
 408 
 409 
 410 
 411 
 412 
 413 
 414 
 415 
 416 
 417 
 418 
 419 
 420 
 421 
 422 
 423 
 424 
 425 
 426 
 427 
 428 
 429 
 430 
 431 
 432 
 433 
 434 
 435 
 436 
 437 /*
 438  * single-precision arithmetic routines ...
 439  */
 440 
 441 #if defined(HEX_MODE)
 442 //#define HEX_SIGN (SIGN72 | BIT71 | BIT70 | BIT69 | BIT68)
 443 //#define HEX_MSB  (         BIT71 | BIT70 | BIT69 | BIT68)
 444 # if defined(NEED_128)
 445 #  define HEX_SIGN construct_128 (0xF0, 0)
 446 #  define HEX_MSB  construct_128 (0x70, 0)
 447 #  define HEX_NORM construct_128 (0x78, 0)
 448 # else
 449 #  define HEX_SIGN (SIGN72 | BIT71 | BIT70 | BIT69)
 450 #  define HEX_MSB  (         BIT71 | BIT70 | BIT69)
 451 #  define HEX_NORM (         BIT71 | BIT70 | BIT69 | BIT68)
 452 # endif
 453 
 454 static inline bool isHex (cpu_state_t * cpup) {
     /*  */
 455     return (! cpu.tweaks.l68_mode) && (!! cpu.options.hex_mode_installed) &&  (!! cpu.MR.hexfp) && (!! TST_I_HEX);
 456 }
 457 #endif
 458 
 459 /*!
 460  * Unnormalized floating single-precision add
 461  */
 462 void ufa (cpu_state_t * cpup, bool sub, bool normalize) {
     /*  */
 463   // C(EAQ) + C(Y) → C(EAQ)
 464   // The ufa instruction is executed as follows:
 465   //
 466   // The mantissas are aligned by shifting the mantissa of the operand having
 467   // the algebraically smaller exponent to the right the number of places
 468   // equal to the absolute value of the difference in the two exponents. Bits
 469   // shifted beyond the bit position equivalent to AQ71 are lost.
 470   //
 471   // The algebraically larger exponent replaces C(E). The sum of the
 472   // mantissas replaces C(AQ).
 473   //
 474   // If an overflow occurs during addition, then;
 475   // * C(AQ) are shifted one place to the right.
 476   // * C(AQ)0 is inverted to restore the sign.
 477   // * C(E) is increased by one.
 478 
 479   CPTUR (cptUseE);
 480 #if defined(HEX_MODE)
 481   uint shift_amt = isHex (cpup) ? 4 : 1;
 482 #else
 483   uint shift_amt = 1;
 484 #endif
 485   word72 m1 = convert_to_word72 (cpu.rA, cpu.rQ);
 486   // 28-bit mantissa (incl sign)
 487 #if defined(NEED_128)
 488   word72 m2 = lshift_128 (construct_128 (0, (uint64_t) getbits36_28 (cpu.CY, 8)), 44u); // 28-bit mantissa (incl sign)
 489 #else
 490   word72 m2 = ((word72) getbits36_28 (cpu.CY, 8)) << 44; // 28-bit mantissa (incl sign)
 491 #endif
 492 
 493   int e1 = SIGNEXT8_int (cpu.rE & MASK8);
 494   int e2 = SIGNEXT8_int (getbits36_8 (cpu.CY, 0));
 495 
 496   // RJ78: The two's complement of the subtrahend is first taken and the
 497   // smaller value is then right-shifted to equalize it (i.e. ufa).
 498 
 499   int m2zero = 0;
 500   if (sub) {
 501     // ISOLTS-735 08i asserts no carry for (-1.0 * 2^96) - (-1.0 * 2^2) but 08g
 502     // asserts carry for -0.5079365*2^78-0.0
 503     // I assume zero subtrahend is handled in a special way.
 504 
 505     if (ISZERO_128 (m2))
 506       m2zero = 1;
 507 #if defined(NEED_128)
 508     if (iseq_128 (m2, SIGN72)) {  // -1.0 -> 0.5, increase exponent, ISOLTS-735 08i,j
 509       m2 = rshift_128 (m2, shift_amt);
 510       e2 += 1;
 511     } else {
 512        m2 = and_128 (negate_128 (m2), MASK72);
 513     }
 514 #else // ! NEED_128
 515     if (m2 == SIGN72) {  // -1.0 -> 0.5, increase exponent, ISOLTS-735 08i,j
 516       m2 >>= shift_amt;
 517       e2 += 1;
 518     } else {
 519       //m2 = (-m2) & MASK72;
 520       // ubsan
 521       m2 = ((word72) (- (word72s) m2)) & MASK72;
 522     }
 523 #endif // ! NEED_128
 524   } // if sub
 525 
 526   int e3 = -1;
 527 
 528   // which exponent is smaller?
 529   L68_ (cpu.ou.cycle |= ou_GOE;)
 530   int shift_count = -1;
 531   word1 allones = 1;
 532   word1 notallzeros = 0;
 533   //word1 last = 0;
 534   (void)allones;
 535   if (e1 == e2) {
 536     shift_count = 0;
 537     e3 = e1;
 538   } else if (e1 < e2) {
 539     shift_count = abs (e2 - e1) * (int) shift_amt;
 540 #if defined(NEED_128)
 541     bool sign = isnonzero_128 (and_128 (m1, SIGN72)); // mantissa negative?
 542     for (int n = 0 ; n < shift_count ; n += 1) {
 543       //last = m1.l & 1;
 544       allones &= m1.l & 1;
 545       notallzeros |= m1.l & 1;
 546       m1 = rshift_128 (m1, 1);
 547       if (sign)
 548         m1 = or_128 (m1, SIGN72);
 549     } // for n
 550     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
 551       m1 = construct_128 (0, 0);
 552     m1 = and_128 (m1, MASK72);
 553     e3 = e2;
 554 #else // ! NEED_128
 555     bool sign = m1 & SIGN72;   // mantissa negative?
 556     for (int n = 0 ; n < shift_count ; n += 1) {
 557       //last = m1 & 1;
 558       allones &= m1 & 1;
 559       notallzeros |= m1 & 1;
 560       m1 >>= 1;
 561       if (sign)
 562         m1 |= SIGN72;
 563     } // for n
 564 
 565     if (m1 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
 566       m1 = 0;
 567     m1 &= MASK72;
 568     e3 = e2;
 569 #endif // NEED_128
 570   } else {
 571     // e2 < e1;
 572     shift_count = abs (e1 - e2) * (int) shift_amt;
 573 #if defined(NEED_128)
 574     bool sign = isnonzero_128 (and_128 (m2, SIGN72)); // mantissa negative?
 575     for (int n = 0 ; n < shift_count ; n += 1) {
 576       //last = m2.l & 1;
 577       allones &= m2.l & 1;
 578       notallzeros |= m2.l & 1;
 579       m2 = rshift_128 (m2, 1);
 580       if (sign)
 581         m2 = or_128 (m2, SIGN72);
 582     } // for n
 583     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
 584       m2 = construct_128 (0, 0);
 585     m2 = and_128 (m2, MASK72);
 586     e3 = e1;
 587 #else // ! NEED_128
 588     bool sign = m2 & SIGN72;   // mantissa negative?
 589     for (int n = 0 ; n < shift_count ; n += 1) {
 590       //last = m2 & 1;
 591       allones &= m2 & 1;
 592       notallzeros |= m2 & 1;
 593       m2 >>= 1;
 594       if (sign)
 595         m2 |= SIGN72;
 596     }
 597     if (m2 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
 598       m2 = 0;
 599     m2 &= MASK72;
 600     e3 = e1;
 601 #endif // ! NEED_128
 602   }
 603 
 604   bool ovf;
 605   word72 m3;
 606   m3 = Add72b (cpup, m1, m2, 0, I_CARRY, & cpu.cu.IR, & ovf);
 607   // ISOLTS-735 08g
 608   // if 2's complement carried, OR it in now.
 609   if (m2zero)
 610     SET_I_CARRY;
 611 
 612   if (ovf) {
 613 #if defined(NEED_128)
 614 # if defined(HEX_MODE)
 615 //    word72 signbit = m3 & sign_msk;
 616 //    m3 >>= shift_amt;
 617 //    m3 = (m3 & MASK71) | signbit;
 618 //    m3 ^= SIGN72; // C(AQ)0 is inverted to restore the sign
 619 //    e3 += 1;
 620     word72 s = and_128 (m3, SIGN72); // save the sign bit
 621     if (isHex (cpup)) {
 622       m3 = rshift_128 (m3, shift_amt); // renormalize the mantissa
 623       if (isnonzero_128 (s))
 624         // Sign is set, number should be positive; clear the sign bit and the 3 MSBs
 625         m3 = and_128 (m3, MASK68);
 626       else
 627         // Sign is clr, number should be negative; set the sign bit and the 3 MSBs
 628         m3 = or_128 (m3, HEX_SIGN);
 629     } else {
 630       word72 signbit = and_128 (m3, SIGN72);
 631       m3 = rshift_128 (m3, 1);
 632       m3 = and_128 (m3, MASK71);
 633       m3 = or_128 (m3, signbit);
 634       m3 = xor_128 (m3, SIGN72); // C(AQ)0 is inverted to restore the sign
 635     }
 636     e3 += 1;
 637 # else // ! HEX_MODE
 638     word72 signbit = and_128 (m3, SIGN72);
 639     m3 = rshift_128 (m3, 1);
 640     m3 = and_128 (m3, MASK71);
 641     m3 = or_128 (m3, signbit);
 642     m3 = xor_128 (m3, SIGN72); // C(AQ)0 is inverted to restore the sign
 643     e3 += 1;
 644 # endif // ! HEX_MODE
 645 #else // ! NEED_128
 646 # if defined(HEX_MODE)
 647 //    word72 signbit = m3 & sign_msk;
 648 //    m3 >>= shift_amt;
 649 //    m3 = (m3 & MASK71) | signbit;
 650 //    m3 ^= SIGN72; // C(AQ)0 is inverted to restore the sign
 651 //    e3 += 1;
 652     word72 s = m3 & SIGN72; // save the sign bit
 653     if (isHex (cpup)) {
 654       m3 >>= shift_amt; // renormalize the mantissa
 655       if (s)
 656         // Sign is set, number should be positive; clear the sign bit and the 3 MSBs
 657         m3 &= MASK68;
 658       else
 659         // Sign is clr, number should be negative; set the sign bit and the 3 MSBs
 660         m3 |=  HEX_SIGN;
 661     } else {
 662       word72 signbit = m3 & SIGN72;
 663       m3 >>= 1;
 664       m3 = (m3 & MASK71) | signbit;
 665       m3 ^= SIGN72; // C(AQ)0 is inverted to restore the sign
 666     }
 667     e3 += 1;
 668 # else // ! HEX_MODE
 669     word72 signbit = m3 & SIGN72;
 670     m3 >>= 1;
 671     m3 = (m3 & MASK71) | signbit;
 672     m3 ^= SIGN72; // C(AQ)0 is inverted to restore the sign
 673     e3 += 1;
 674 # endif // ! HEX_MODE
 675 #endif // ! NEED_128
 676   }
 677 
 678   convert_to_word36 (m3, & cpu.rA, & cpu.rQ);
 679 #if defined(TESTING)
 680   HDBGRegAW ("ufa");
 681 #endif
 682   cpu.rE = e3 & 0377;
 683 
 684   SC_I_NEG (cpu.rA & SIGN36); // Do this here instead of in Add72b because
 685                                 // of ovf handling above
 686   if (cpu.rA == 0 && cpu.rQ == 0) {
 687     SET_I_ZERO;
 688     cpu.rE = 0200U; /*-128*/
 689   } else {
 690     CLR_I_ZERO;
 691   }
 692 
 693   if (normalize) {
 694     fno_ext (cpup, & e3, & cpu.rE, & cpu.rA, & cpu.rQ);
 695   }
 696 
 697   // EOFL: If exponent is greater than +127, then ON
 698   if (e3 > 127) {
 699     SET_I_EOFL;
 700     if (tstOVFfault (cpup))
 701       dlyDoFault (FAULT_OFL, fst_zero, "ufa exp overflow fault");
 702   }
 703 
 704   // EUFL: If exponent is less than -128, then ON
 705   if (e3 < -128) {
 706     SET_I_EUFL;
 707     if (tstOVFfault (cpup))
 708       dlyDoFault (FAULT_OFL, fst_zero, "ufa exp underflow fault");
 709   }
 710 } // ufa
 711 
 712 /*
 713  * floating normalize ...
 714  */
 715 
 716 void fno (cpu_state_t * cpup, word8 * E, word36 * A, word36 * Q) {
     /*  */
 717   int e0 = SIGNEXT8_int ((* E) & MASK8);
 718   fno_ext (cpup, & e0, E, A, Q);
 719 }
 720 
 721 void fno_ext (cpu_state_t * cpup, int * e0, word8 * E, word36 * A, word36 * Q) {
     /*  */
 722   // The fno instruction normalizes the number in C(EAQ) if C(AQ) ≠ 0 and the
 723   // overflow indicator is OFF.
 724   //
 725   // A normalized floating number is defined as one whose mantissa lies in
 726   // the interval [0.5,1.0) such that
 727   //     0.5<= |C(AQ)| <1.0 which, in turn, requires that C(AQ)0 ≠ C(AQ)1.
 728   //
 729   // If the overflow indicator is ON, then C(AQ) is shifted one place to the
 730   // right, C(AQ)0 is inverted to reconstitute the actual sign, and the
 731   // overflow indicator is set OFF. This action makes the fno instruction
 732   // useful in correcting overflows that occur with fixed point numbers.
 733   //
 734   // Normalization is performed by shifting C(AQ)1,71 one place to the left
 735   // and reducing C(E) by 1, repeatedly, until the conditions for C(AQ)0 and
 736   // C(AQ)1 are met. Bits shifted out of AQ1 are lost.
 737   //
 738   // If C(AQ) = 0, then C(E) is set to -128 and the zero indicator is set ON.
 739 
 740   L68_ (cpu.ou.cycle |= ou_GON;)
 741 #if defined(HEX_MODE)
 742   uint shift_amt = isHex(cpup) ? 4 : 1;
 743 #endif
 744   * A &= DMASK;
 745   * Q &= DMASK;
 746   float72 m = convert_to_word72 (* A, * Q);
 747   if (TST_I_OFLOW) {
 748     CLR_I_OFLOW;
 749 #if defined(NEED_128)
 750     word72 s = and_128 (m, SIGN72); // save the sign bit
 751 # if defined(HEX_MODE)
 752     if (isHex (cpup)) {
 753       m = rshift_128 (m, shift_amt); // renormalize the mantissa
 754       if (isnonzero_128 (s)) // sign of mantissa
 755         // Sign is set, number should be positive; clear the sign bit and the 3 MSBs
 756         m = and_128 (m, MASK68);
 757       else
 758         // Sign is clr, number should be negative; set the sign bit and the 3 MSBs
 759         m = or_128 (m, HEX_SIGN);
 760     } else {
 761       m = rshift_128 (m, 1); // renormalize the mantissa
 762       m = or_128 (m, SIGN72); // set the sign bit
 763       m = xor_128 (m, s); // if the was 0, leave it 1; if it was 1, make it 0
 764     }
 765 # else // ! HEX_MODE
 766   m = rshift_128 (m, 1); // renormalize the mantissa
 767   m = or_128 (m, SIGN72); // set the sign bit
 768   m = xor_128 (m, s); // if the was 0, leave it 1; if it was 1, make it 0
 769 # endif // ! HEX_MODE
 770 
 771   // Zero: If C(AQ) = floating point 0, then ON; otherwise OFF
 772   if (iszero_128 (m)) {
 773     * E = 0200U; /* -128 */
 774     SET_I_ZERO;
 775   } else {
 776     CLR_I_ZERO;
 777     if (* E == 127) {
 778 //sim_printf ("overflow 1 E %o \r\n", * E);
 779       SET_I_EOFL;
 780       if (tstOVFfault (cpup))
 781         dlyDoFault (FAULT_OFL, fst_zero, "fno exp overflow fault");
 782     }
 783     (* E) ++;
 784     * E &= MASK8;
 785   }
 786 #else // ! NEED_128
 787   word72 s = m & SIGN72; // save the sign bit
 788 # if defined(HEX_MODE)
 789   if (isHex (cpup)) {
 790     m >>= shift_amt; // renormalize the mantissa
 791     if (s)
 792       // Sign is set, number should be positive; clear the sign bit and the 3 MSBs
 793       m &= MASK68;
 794     else
 795       // Sign is clr, number should be negative; set the sign bit and the 3 MSBs
 796       m |=  HEX_SIGN;
 797     } else { // ! isHex
 798       m >>= 1; // renormalize the mantissa
 799       m |= SIGN72; // set the sign bit
 800       m ^= s; // if the was 0, leave it 1; if it was 1, make it 0
 801     } // ! isHex
 802 # else // ! HEX_MODE
 803     m >>= 1; // renormalize the mantissa
 804     m |= SIGN72; // set the sign bit
 805     m ^= s; // if the was 0, leave it 1; if it was 1, make it 0
 806 # endif // ! HEX_MODE
 807 
 808     // Zero: If C(AQ) = floating point 0, then ON; otherwise OFF
 809     if (m == 0) {
 810       * E = 0200U; /* -128 */
 811       SET_I_ZERO;
 812     } else {
 813       CLR_I_ZERO;
 814       if (* E == 127) {
 815         SET_I_EOFL;
 816         if (tstOVFfault (cpup))
 817           dlyDoFault (FAULT_OFL, fst_zero, "fno exp overflow fault");
 818       }
 819       (* E) ++;
 820       * E &= MASK8;
 821     }
 822 #endif // ! NEED_128
 823     convert_to_word36 (m, A, Q);
 824     SC_I_NEG ((* A) & SIGN36);
 825 
 826     return;
 827   }
 828 
 829   // only normalize C(EAQ) if C(AQ) ≠ 0 and the overflow indicator is OFF
 830 #if defined(NEED_128)
 831   if (iszero_128 (m)) { /// C(AQ) == 0.
 832     //* A = (m >> 36) & MASK36;
 833     //* Q = m & MASK36;
 834     * E = 0200U; /* -128 */
 835     SET_I_ZERO;
 836     CLR_I_NEG;
 837     return;
 838   }
 839 #else // ! NEED_128
 840   if (m == 0) { // C(AQ) == 0.
 841     //* A = (m >> 36) & MASK36;
 842     //* Q = m & MASK36;
 843     * E = 0200U; /* -128 */
 844     SET_I_ZERO;
 845     CLR_I_NEG;
 846     return;
 847   }
 848 #endif // ! NEED_128
 849 
 850   //int e = SIGNEXT8_int ((* E) & MASK8);
 851   int e = * e0;
 852 #if defined(NEED_128)
 853   bool s = isnonzero_128 (and_128 (m, SIGN72));
 854 #else
 855   bool s = (m & SIGN72) != (word72)0;    ///< save sign bit
 856 #endif
 857 
 858 #if defined(HEX_MODE)
 859 // Normalized in Hex Mode: If sign is 0, bits 1-4 != 0; if sign is 1,
 860 // bits 1-4 != 017.
 861   if (isHex (cpup)) {
 862     if (s) {
 863       // Negative
 864       // Until bits 1-4 != 014
 865       // Termination guarantee: Zeros are being shifted into the right
 866       // end, so the loop will terminate when the first shifted
 867       // zero enters bits 1-4.
 868 # if defined(NEED_128)
 869       while (iseq_128 (and_128 (m, HEX_NORM), HEX_NORM)) {
 870         m = lshift_128 (m, 4);
 871         e -= 1;
 872       }
 873       m = and_128 (m, MASK71);
 874       m = or_128 (m, SIGN72);
 875 # else // ! NEED_128
 876       while ((m & HEX_NORM) == HEX_NORM) {
 877         m <<= 4;
 878         e -= 1;
 879       }
 880       m &= MASK71;
 881       m |= SIGN72;
 882 # endif // ! NEED_128
 883     } else { // ! s
 884 # if defined(NEED_128)
 885       // Positive
 886       // Until bits 1-4 != 0
 887       // Termination guarantee: m is known to be non-zero; a non-zero
 888       // bit will eventually be shifted into bits 1-4.
 889       while (iszero_128 (and_128 (m, HEX_NORM))) {
 890         m = lshift_128 (m, 4);
 891         e -= 1;
 892       }
 893       m = and_128 (m, MASK71);
 894 # else // ! NEED_128
 895       // Positive
 896       // Until bits 1-4 != 0
 897       // Termination guarantee: m is known to be non-zero; a non-zero
 898       // bit will eventually be shifted into bits 1-4.
 899       while ((m & HEX_NORM) == 0) {
 900         m <<= 4;
 901         e -= 1;
 902       }
 903       m &= MASK71;
 904 # endif // ! NEED_128
 905     } // ! s
 906   } else { // ! isHex
 907 # if defined(NEED_128)
 908     while (s == isnonzero_128 (and_128 (m, BIT71))) { // until C(AQ)0 != C(AQ)1?
 909       m = lshift_128 (m, 1);
 910       e -= 1;
 911     }
 912 
 913     m = and_128 (m, MASK71);
 914 
 915     if (s)
 916       m = or_128 (m, SIGN72);
 917 # else // ! NEED_128
 918     while (s == !! (m & BIT71)) { // until C(AQ)0 != C(AQ)1?
 919       m <<= 1;
 920       e -= 1;
 921     }
 922 
 923     m &= MASK71;
 924 
 925     if (s)
 926       m |= SIGN72;
 927 # endif // ! NEED_128
 928   } // ! isHex
 929 #else // ! HEX_MODE
 930 # if defined(NEED_128)
 931   while (s == isnonzero_128 (and_128 (m, BIT71))) { // until C(AQ)0 != C(AQ)1?
 932     m = lshift_128 (m, 1);
 933     e -= 1;
 934   }
 935 
 936   m = and_128 (m, MASK71);
 937 
 938   if (s)
 939     m = or_128 (m, SIGN72);
 940 # else // ! NEED_128
 941   while (s == !! (m & BIT71)) { // until C(AQ)0 != C(AQ)1?
 942     m <<= 1;
 943     e -= 1;
 944   }
 945 
 946   m &= MASK71;
 947 
 948   if (s)
 949     m |= SIGN72;
 950 # endif // ! NEED_128
 951 #endif // ! HEX_MODE
 952 
 953   if (e < -128) {
 954     SET_I_EUFL;
 955     if (tstOVFfault (cpup))
 956       dlyDoFault (FAULT_OFL, fst_zero, "fno exp underflow fault");
 957   }
 958 
 959   * E = (word8) e & MASK8;
 960   convert_to_word36 (m, A, Q);
 961 
 962   // EAQ is normalized, so if A is 0, so is Q, and the check can be elided
 963   if (* A == 0)    // set to normalized 0
 964     * E = 0200U; /* -128 */
 965 
 966   // Zero: If C(AQ) = floating point 0, then ON; otherwise OFF
 967   SC_I_ZERO (* A == 0 && * Q == 0);
 968 
 969   // Neg: If C(AQ)0 = 1, then ON; otherwise OFF
 970   SC_I_NEG ((* A) & SIGN36);
 971 
 972   * e0 = e;
 973 }
 974 
 975 
 976 
 977 
 978 
 979 
 980 
 981 
 982 
 983 
 984 
 985 
 986 
 987 
 988 
 989 
 990 
 991 
 992 
 993 
 994 
 995 
 996 
 997 
 998 
 999 
1000 
1001 
1002 
1003 
1004 
1005 
1006 
1007 
1008 
1009 
1010 
1011 
1012 
1013 
1014 
1015 
1016 
1017 
1018 
1019 
1020 
1021 
1022 
1023 
1024 
1025 
1026 
1027 
1028 
1029 
1030 
1031 
1032 
1033 
1034 
1035 
1036 
1037 
1038 
1039 
1040 
1041 
1042 
1043 
1044 
1045 
1046 
1047 
1048 
1049 
1050 
1051 
1052 
1053 
1054 
1055 
1056 
1057 
1058 
1059 
1060 
1061 
1062 
1063 
1064 
1065 
1066 
1067 /*
1068  * floating negate ...
1069  */
1070 void fneg (cpu_state_t * cpup) {
     /*  */
1071   // This instruction changes the number in C(EAQ) to its normalized negative
1072   // (if C(AQ) ≠ 0). The operation is executed by first forming the twos
1073   // complement of C(AQ), and then normalizing C(EAQ).
1074   //
1075   // Even if originally C(EAQ) were normalized, an exponent overflow can
1076   // still occur, namely when C(E) = +127 and C(AQ) = 100...0 which is the
1077   // twos complement approximation for the decimal value -1.0.
1078 
1079   CPTUR (cptUseE);
1080   // Form the mantissa from AQ
1081 
1082   word72 m = convert_to_word72 (cpu.rA, cpu.rQ);
1083 
1084   // If the mantissa is 4000...00 (least negative value, it is negable in
1085   // two's complement arithmetic. Divide it by 2, losing a bit of precision,
1086   // and increment the exponent.
1087   if (ISEQ_128 (m, SIGN72)) {
1088     // Negation of 400..0 / 2 is 200..0; we can get there shifting; we know
1089     // that a zero will be shifted into the sign bit because of the masking
1090     // in 'm='.
1091 #if defined(NEED_128)
1092     m = rshift_128 (m, 1);
1093 #else
1094     m >>= 1;
1095 #endif
1096     // Increment the exp, checking for overflow.
1097     if (cpu.rE == 127) {
1098       SET_I_EOFL;
1099       if (tstOVFfault (cpup))
1100         dlyDoFault (FAULT_OFL, fst_zero, "fneg exp overflow fault");
1101     }
1102     cpu.rE ++;
1103     cpu.rE &= MASK8;
1104   } else {
1105     // Do the negation
1106 #if defined(NEED_128)
1107     m = negate_128 (m);
1108 #else
1109     // m = -m;
1110     // ubsan
1111     m = (word72) (- (word72s) m);
1112 #endif
1113   }
1114   convert_to_word36 (m, & cpu.rA, & cpu.rQ);
1115   fno (cpup, & cpu.rE, & cpu.rA, & cpu.rQ);  // normalize
1116 #if defined(TESTING)
1117   HDBGRegAW ("fneg");
1118 #endif
1119 }
1120 
1121 /*
1122  * Unnormalized Floating Multiply ...
1123  */
1124 void ufm (cpu_state_t * cpup, bool normalize) {
     /*  */
1125   // The ufm instruction is executed as follows:
1126   //      C(E) + C(Y)0,7 → C(E)
1127   //      ( C(AQ) × C(Y)8,35 )0,71 → C(AQ)
1128 
1129   // Zero: If C(AQ) = 0, then ON; otherwise OFF
1130   // Neg: If C(AQ)0 = 1, then ON; otherwise OFF
1131   // Exp Ovr: If exponent is greater than +127, then ON
1132   // Exp Undr: If exponent is less than -128, then ON
1133 
1134   CPTUR (cptUseE);
1135   word72 m1 = convert_to_word72 (cpu.rA, cpu.rQ);
1136   int    e1 = SIGNEXT8_int (cpu.rE & MASK8);
1137 #if defined(NEED_128)
1138   word72 m2 = lshift_128 (construct_128 (0, (uint64_t) getbits36_28 (cpu.CY, 8)), 44u); // 28-bit mantissa (incl sign)
1139 #else
1140   word72 m2 = ((word72) getbits36_28 (cpu.CY, 8)) << 44; ///< 28-bit mantissa (incl sign)
1141 #endif
1142   int    e2 = SIGNEXT8_int (getbits36_8 (cpu.CY, 0));
1143 
1144   if (ISZERO_128 (m1) || ISZERO_128 (m2)) {
1145     SET_I_ZERO;
1146     CLR_I_NEG;
1147 
1148     cpu.rE = 0200U; /*-128*/
1149     cpu.rA = 0;
1150 #if defined(TESTING)
1151     HDBGRegAW ("ufm");
1152 #endif
1153     cpu.rQ = 0;
1154 
1155     return; // normalized 0
1156   }
1157 
1158   int e3 = e1 + e2;
1159 
1160 
1161 
1162 
1163 
1164 
1165 
1166 
1167 
1168 
1169 
1170 
1171 
1172 
1173   // RJ78: This multiplication is executed in the following way:
1174   // C(E) + C(Y)(0-7) -> C(E)
1175   // C(AQ) * C(Y)(8-35) results in a 98-bit product plus sign,
1176   // the leading 71 bits plus sign of which -> C(AQ).
1177 
1178   // shift the CY mantissa to get 98 bits precision
1179 #if defined(NEED_128)
1180   int128 t = SIGNEXT72_128(m2);
1181   uint128 ut = rshift_128 (cast_128 (t), 44);
1182   int128 m3 = multiply_s128 (SIGNEXT72_128(m1), cast_s128 (ut));
1183 //sim_debug (DBG_TRACEEXT, & cpu_dev, "m3 %016"PRIx64"%016"PRIx64"\n", m3.h, m3.l);
1184 #else
1185   int128 m3 = (SIGNEXT72_128(m1) * (SIGNEXT72_128(m2) >> 44));
1186 //sim_debug (DBG_TRACEEXT, & cpu_dev, "m3 %016"PRIx64"%016"PRIx64"\n", (uint64_t) (m3>>64), (uint64_t) m3);
1187 #endif
1188   // realign to 72bits
1189 #if defined(NEED_128)
1190   word72 m3a = and_128 (rshift_128 (cast_128 (m3), 98u - 71u), MASK72);
1191 //sim_debug (DBG_TRACEEXT, & cpu_dev, "m3a %016"PRIx64"%016"PRIx64"\n", m3a.h, m3a.l);
1192 #else
1193   word72 m3a = ((word72) (m3 >> (98-71))) & MASK72;
1194 //sim_debug (DBG_TRACEEXT, & cpu_dev, "m3a %016"PRIx64"%016"PRIx64"\n", (uint64_t) (m3a>>64), (uint64_t) m3a);
1195 #endif
1196 
1197   // A normalization is performed only in the case of both factor mantissas being 100...0
1198   // which is the twos complement approximation to the decimal value -1.0.
1199   if (ISEQ_128 (m1, SIGN72) && ISEQ_128 (m2, SIGN72)) {
1200     if (e3 == 127) {
1201       SET_I_EOFL;
1202       if (tstOVFfault (cpup))
1203         dlyDoFault (FAULT_OFL, fst_zero, "ufm exp overflow fault");
1204     }
1205 #if defined(NEED_128)
1206     m3a = rshift_128 (m3a, 1);
1207 #else
1208     m3a >>= 1;
1209 #endif
1210     e3 += 1;
1211   }
1212 
1213   convert_to_word36 (m3a, & cpu.rA, & cpu.rQ);
1214 #if defined(TESTING)
1215   HDBGRegAW ("ufm");
1216 #endif
1217   cpu.rE = (word8) e3 & MASK8;
1218 sim_debug (DBG_TRACEEXT, & cpu_dev, "fmp A %012"PRIo64" Q %012"PRIo64" E %03o\n", cpu.rA, cpu.rQ, cpu.rE);
1219   SC_I_NEG (cpu.rA & SIGN36);
1220 
1221   if (cpu.rA == 0 && cpu.rQ == 0) {
1222     SET_I_ZERO;
1223     cpu.rE = 0200U; /*-128*/
1224   } else {
1225     CLR_I_ZERO;
1226   }
1227 
1228   if (normalize) {
1229     fno_ext (cpup, & e3, & cpu.rE, & cpu.rA, & cpu.rQ);
1230   }
1231 
1232   if (e3 >  127) {
1233     SET_I_EOFL;
1234     if (tstOVFfault (cpup))
1235       dlyDoFault (FAULT_OFL, fst_zero, "ufm exp overflow fault");
1236   }
1237   if (e3 < -128) {
1238     SET_I_EUFL;
1239     if (tstOVFfault (cpup))
1240       dlyDoFault (FAULT_OFL, fst_zero, "ufm exp underflow fault");
1241   }
1242 }
1243 
1244 /*
1245  * floating divide ...
1246  */
1247 // CANFAULT
1248 static void fdvX (cpu_state_t * cpup, bool bInvert) {
     /*  */
1249   // C(EAQ) / C (Y) → C(EA)
1250   // C(Y) / C(EAQ) → C(EA) (Inverted)
1251 
1252   // 00...0 → C(Q)
1253 
1254   // The fdv instruction is executed as follows:
1255   // The dividend mantissa C(AQ) is shifted right and the dividend exponent
1256   // C(E) increased accordingly until | C(AQ)0,27 | < | C(Y)8,35 |
1257   // C(E) - C(Y)0,7 → C(E)
1258   // C(AQ) / C(Y)8,35 → C(A)
1259   // 00...0 → C(Q)
1260 
1261   CPTUR (cptUseE);
1262 #if defined(HEX_MODE)
1263   uint shift_amt = isHex(cpup) ? 4 : 1;
1264 #else
1265   uint shift_amt = 1;
1266 #endif
1267   word72 m1;
1268   int    e1;
1269 
1270   word72 m2;
1271   int    e2;
1272 
1273   bool roundovf = 0;
1274 
1275   if (! bInvert) {
1276     m1 = convert_to_word72 (cpu.rA, cpu.rQ);
1277     e1 = SIGNEXT8_int (cpu.rE & MASK8);
1278 
1279 #if defined(NEED_128)
1280     m2 = lshift_128 (construct_128 (0, (uint64_t) getbits36_28 (cpu.CY, 8)), 44u); // 28-bit mantissa (incl sign)
1281 #else
1282     m2 = ((word72) getbits36_28 (cpu.CY, 8)) << 44; ///< 28-bit mantissa (incl sign)
1283 #endif
1284     e2 = SIGNEXT8_int (getbits36_8 (cpu.CY, 0));
1285 
1286   } else { // invert
1287 
1288     m2 = convert_to_word72 (cpu.rA, cpu.rQ);
1289     e2 = SIGNEXT8_int (cpu.rE & MASK8);
1290 
1291     // round divisor per RJ78
1292     // If AQ(28-71) is not equal to 0 and A(0) = 0, then 1 is added to AQ(27).
1293     // 0 -> AQ(28-71) unconditionally. AQ(0-27) is then used as the divisor mantissa.
1294 #if defined(NEED_128)
1295     if ((iszero_128 (and_128 (m2, SIGN72))) &&
1296         (isnonzero_128 (and_128 (m2, construct_128 (0, 0377777777777777LL))))) {
1297       m2  = add_128 (m2, construct_128 (0, 0400000000000000LL));
1298       // I surmise that the divisor is taken as unsigned 28 bits in this case
1299       roundovf = 1;
1300     }
1301     m2 = and_128 (m2, lshift_128 (construct_128 (0, 0777777777400), 36));
1302 
1303     m1 = lshift_128 (construct_128 (0, getbits36_28 (cpu.CY, 8)), 44); ///< 28-bit mantissa (incl sign)
1304     e1 = SIGNEXT8_int (getbits36_8 (cpu.CY, 0));
1305 #else
1306     if (!(m2 & SIGN72) && m2 & 0377777777777777LL) {
1307         m2 += 0400000000000000LL;
1308         // I surmise that the divisor is taken as unsigned 28 bits in this case
1309       roundovf = 1;
1310     }
1311     m2 &= (word72)0777777777400 << 36;
1312 
1313     m1 = ((word72) getbits36_28 (cpu.CY, 8)) << 44; ///< 28-bit mantissa (incl sign)
1314     e1 = SIGNEXT8_int (getbits36_8 (cpu.CY, 0));
1315 #endif
1316   }
1317 
1318   if (ISZERO_128 (m1)) {
1319     SET_I_ZERO;
1320     CLR_I_NEG;
1321 
1322     cpu.rE = 0200U; /*-128*/
1323     cpu.rA = 0;
1324 #if defined(TESTING)
1325     HDBGRegAW ("fdvX");
1326 #endif
1327     cpu.rQ = 0;
1328 
1329     return; // normalized 0
1330   }
1331 
1332   // make everything positive, but save sign info for later....
1333   int sign = 1;
1334 #if defined(NEED_128)
1335   if (isnonzero_128 (and_128 (m1, SIGN72)))
1336 #else
1337   if (m1 & SIGN72)
1338 #endif
1339   {
1340     SET_I_NEG; // in case of divide fault
1341 #if defined(NEED_128)
1342     if (iseq_128 (m1, SIGN72)) {
1343       m1 = rshift_128 (m1, shift_amt);
1344       e1 += 1;
1345     } else {
1346       m1 = and_128 (negate_128 (m1), MASK72);
1347     }
1348 #else
1349     if (m1 == SIGN72) {
1350       m1 >>= shift_amt;
1351       e1 += 1;
1352     } else {
1353       m1 = (~m1 + 1) & MASK72;
1354     }
1355 #endif
1356     sign = -sign;
1357   } else {
1358     CLR_I_NEG; // in case of divide fault
1359   }
1360 
1361 #if defined(NEED_128)
1362   if ((isnonzero_128 (and_128 (m2, SIGN72))) && !roundovf) {
1363     if (iseq_128 (m2, SIGN72)) {
1364       m2 = rshift_128 (m2, shift_amt);
1365       e2 += 1;
1366     } else {
1367       m2 = and_128 (negate_128 (m2), MASK72);
1368     }
1369     sign = -sign;
1370   }
1371 #else
1372   if ((m2 & SIGN72) && !roundovf) {
1373     if (m2 == SIGN72) {
1374       m2 >>= shift_amt;
1375       e2 += 1;
1376     } else {
1377       m2 = (~m2 + 1) & MASK72;
1378     }
1379     sign = -sign;
1380   }
1381 #endif
1382 
1383   if (ISZERO_128 (m2)) {
1384     // NB: If C(Y)8,35 ==0 then the alignment loop will never exit! That's why it been moved before the alignment
1385 
1386     SET_I_ZERO;
1387     // NEG already set
1388 
1389     // FDV: If the divisor mantissa C(Y)8,35 is zero after alignment (HWR: why after?), the division does
1390     // not take place. Instead, a divide check fault occurs, C(AQ) contains the dividend magnitude, and
1391     // the negative indicator reflects the dividend sign.
1392     // FDI: If the divisor mantissa C(AQ) is zero, the division does not take place.
1393     // Instead, a divide check fault occurs and all the registers remain unchanged.
1394     if (! bInvert) {
1395       convert_to_word36 (m1, & cpu.rA, & cpu.rQ);
1396 #if defined(TESTING)
1397       HDBGRegAW ("fdvX");
1398 #endif
1399     }
1400 
1401     doFault (FAULT_DIV, fst_zero, "FDV: divide check fault");
1402   }
1403 
1404 #if defined(NEED_128)
1405   while (isge_128 (m1, m2)) {
1406     // DH02 (equivalent but perhaps clearer description):
1407     // dividend exponent C(E) increased accordingly until | C(AQ)0,71 | < | C(Y)8,35 with zero fill |
1408     // We have already taken the absolute value so just shift it
1409     m1 = rshift_128 (m1, shift_amt);
1410     e1 += 1;
1411   }
1412 #else
1413   while (m1 >= m2) {
1414     // DH02 (equivalent but perhaps clearer description):
1415     // dividend exponent C(E) increased accordingly until | C(AQ)0,71 | < | C(Y)8,35 with zero fill |
1416     // We have already taken the absolute value so just shift it
1417     m1 >>= shift_amt;
1418     e1 += 1;
1419   }
1420 #endif
1421 
1422   int e3 = e1 - e2;
1423 
1424   if (e3 > 127) {
1425     SET_I_EOFL;
1426     if (tstOVFfault (cpup))
1427       dlyDoFault (FAULT_OFL, fst_zero, "fdvX exp overflow fault");
1428   }
1429   if (e3 < -128) {
1430     SET_I_EUFL;
1431     if (tstOVFfault (cpup))
1432       dlyDoFault (FAULT_OFL, fst_zero, "fdvX exp underflow fault");
1433   }
1434 
1435   // We need 35 bits quotient + sign. Divisor is at most 28 bits.
1436   // Do a 63(28+35) by 35 fractional divide
1437   // lo 44bits are always zero
1438 #if defined(NEED_128)
1439   uint32_t divisor = rshift_128 (m2, 44).l & MASK28;
1440   word72 m3;
1441   if (divisor > MASK16)
1442     m3 = divide_128_32 (rshift_128 (m1, (44u-35u)), divisor, NULL);
1443   else
1444     m3 = divide_128_16 (rshift_128 (m1, (44u-35u)), (uint16_t) divisor, NULL);
1445 #else
1446   word72 m3 = (m1 >> (44-35)) / (m2 >> 44);
1447 #endif
1448 
1449 #if defined(NEED_128)
1450   m3 = lshift_128 (m3, 36);
1451   if (sign == -1)
1452     m3 = and_128 (negate_128 (m3), MASK72);
1453 #else
1454   m3 <<= 36; // convert back to float
1455   if (sign == -1)
1456     m3 = (~m3 + 1) & MASK72;
1457 #endif
1458 
1459   cpu.rE = (word8) e3 & MASK8;
1460 #if defined(NEED_128)
1461   cpu.rA = rshift_128 (m3, 36u).l & MASK36;
1462 #else
1463   cpu.rA = (m3 >> 36) & MASK36;
1464 #endif
1465 #if defined(TESTING)
1466   HDBGRegAW ("fdvX");
1467 #endif
1468   cpu.rQ = 0;
1469 
1470   SC_I_ZERO (cpu . rA == 0);
1471   SC_I_NEG (cpu . rA & SIGN36);
1472 
1473   if (cpu.rA == 0)    // set to normalized 0
1474     cpu.rE = 0200U; /*-128*/
1475 } // fdvX
1476 
1477 void fdv (cpu_state_t * cpup) {
     /*  */
1478   fdvX (cpup, false);    // no inversion
1479 }
1480 
1481 void fdi (cpu_state_t * cpup) {
     /*  */
1482   fdvX (cpup, true);      // invert
1483 }
1484 
1485 /*
1486  * single precision floating round ...
1487  */
1488 void frd (cpu_state_t * cpup) {
     /*  */
1489   // If C(AQ) != 0, the frd instruction performs a true round to a precision
1490   // of 28 bits and a normalization on C(EAQ).
1491   // A true round is a rounding operation such that the sum of the result of
1492   // applying the operation to two numbers of equal magnitude but opposite
1493   // sign is exactly zero.
1494 
1495   // The frd instruction is executed as follows:
1496   // C(AQ) + (11...1)29,71 -> C(AQ)
1497   // If C(AQ)0 = 0, then a carry is added at AQ71
1498   // If overflow occurs, C(AQ) is shifted one place to the right and C(E) is
1499   // increased by 1.
1500   // If overflow does not occur, C(EAQ) is normalized.
1501   // If C(AQ) = 0, C(E) is set to -128 and the zero indicator is set ON.
1502 
1503   // I believe AL39 is incorrect; bits 28-71 should be set to 0, not 29-71.
1504   // DH02-01 & Bull 9000 is correct.
1505 
1506   // test case 15.5
1507   //                 rE                     rA                                     rQ
1508   // 014174000000 00000110 000111110000000000000000000000000000 000000000000000000000000000000000000
1509   // +                                                  1111111 111111111111111111111111111111111111
1510   // =            00000110 000111110000000000000000000001111111 111111111111111111111111111111111111
1511   // If C(AQ)0 = 0, then a carry is added at AQ71
1512   // =            00000110 000111110000000000000000000010000000 000000000000000000000000000000000000
1513   // 0 → C(AQ)29,71
1514   //              00000110 000111110000000000000000000010000000 000000000000000000000000000000000000
1515   // after normalization .....
1516   // 010760000002 00000100 011111000000000000000000001000000000 000000000000000000000000000000000000
1517   // I think this is wrong
1518 
1519   // 0 -> C(AQ)28,71
1520   //              00000110 000111110000000000000000000000000000 000000000000000000000000000000000000
1521   // after normalization .....
1522   // 010760000000 00000100 011111000000000000000000000000000000 000000000000000000000000000000000000
1523   // which I think is correct
1524 
1525   //
1526   // GE CPB1004F, DH02-01 (DPS8/88) & Bull DPS9000 assembly ... have this ...
1527 
1528   // The rounding operation is performed in the following way.
1529   // -  a) A constant (all 1s) is added to bits 29-71 of the mantissa.
1530   // -  b) If the number being rounded is positive, a carry is inserted into
1531   // the least significant bit position of the adder.
1532   // -  c) If the number being rounded is negative, the carry is not inserted.
1533   // -  d) Bits 28-71 of C(AQ) are replaced by zeros.
1534   // If the mantissa overflows upon rounding, it is shifted right one place
1535   // and a corresponding correction is made to the exponent.
1536   // If the mantissa does not overflow and is nonzero upon rounding,
1537   // normalization is performed.
1538 
1539   // If the resultant mantissa is all zeros, the exponent is forced to -128
1540   // and the zero indicator is set.
1541   // If the exponent resulting from the operation is greater than +127, the
1542   // exponent Overflow indicator is set.
1543   // If the exponent resulting from the operation is less than -128, the
1544   // exponent Underflow indicator is set.
1545   // The definition of normalization is located under the description of the FNO instruction.
1546 
1547   // So, Either AL39 is wrong or the DPS8m did it wrong. (Which was fixed in
1548   // later models.) I'll assume AL39 is wrong.
1549 
1550   CPTUR (cptUseE);
1551   L68_ (cpu.ou.cycle |= ou_GOS;)
1552 
1553   word72 m = convert_to_word72 (cpu.rA, cpu.rQ);
1554   if (ISZERO_128 (m)) {
1555     cpu.rE = 0200U; /*-128*/
1556     SET_I_ZERO;
1557     CLR_I_NEG;
1558     return;
1559   }
1560 
1561   // C(AQ) + (11...1)29,71 → C(AQ)
1562   bool ovf;
1563   word18 flags1 = 0;
1564   word1 carry = 0;
1565   // If C(AQ)0 = 0, then a carry is added at AQ71
1566 #if defined(NEED_128)
1567   if (iszero_128 (and_128 (m, SIGN72)))
1568 #else
1569   if ((m & SIGN72) == 0)
1570 #endif
1571   {
1572     carry = 1;
1573   }
1574 #if defined(NEED_128)
1575   m = Add72b (cpup, m, construct_128 (0, 0177777777777777LL), carry, I_OFLOW, & flags1, & ovf);
1576 #else
1577   m = Add72b (cpup, m, 0177777777777777LL, carry, I_OFLOW, & flags1, & ovf);
1578 #endif
1579 
1580   // 0 -> C(AQ)28,71  (per. RJ78)
1581 #if defined(NEED_128)
1582   m = and_128 (m, lshift_128 (construct_128 (0, 0777777777400), 36));
1583 #else
1584   m &= ((word72)0777777777400 << 36);
1585 #endif
1586 
1587   // If overflow occurs, C(AQ) is shifted one place to the right and C(E) is
1588   // increased by 1.
1589   // If overflow does not occur, C(EAQ) is normalized.
1590   // All of this is done by fno, we just need to save the overflow flag
1591 
1592   bool savedovf = TST_I_OFLOW;
1593   SC_I_OFLOW(ovf);
1594   convert_to_word36 (m, & cpu.rA, & cpu.rQ);
1595 
1596   fno (cpup, & cpu.rE, & cpu.rA, & cpu.rQ);
1597 #if defined(TESTING)
1598   HDBGRegAW ("frd");
1599 #endif
1600   SC_I_OFLOW(savedovf);
1601 }
1602 
1603 void fstr (cpu_state_t * cpup, word36 *Y)
     /*  */
1604   {
1605     // The fstr instruction performs a true round and normalization on C(EAQ)
1606     // as it is stored.  The definition of true round is located under the
1607     // description of the frd instruction.  The definition of normalization is
1608     // located under the description of the fno instruction.  Attempted
1609     // repetition with the rpl instruction causes an illegal procedure fault.
1610 
1611     L68_ (cpu.ou.cycle |= ou_GOS;)
1612     word36 A = cpu . rA, Q = cpu . rQ;
1613     word8 E = cpu . rE;
1614     //A &= DMASK;
1615     //Q &= DMASK;
1616     //E &= (int) MASK8;
1617 
1618     float72 m = convert_to_word72 (A, Q);
1619 #if defined(NEED_128)
1620     if (iszero_128 (m))
1621 #else
1622     if (m == 0)
1623 #endif
1624       {
1625         E = (word8)-128;
1626         SET_I_ZERO;
1627         CLR_I_NEG;
1628         *Y = 0;
1629         putbits36_8 (Y, 0, (word8) E & MASK8);
1630         return;
1631       }
1632 
1633     // C(AQ) + (11...1)29,71 → C(AQ)
1634     bool ovf;
1635     word18 flags1 = 0;
1636     word1 carry = 0;
1637     // If C(AQ)0 = 0, then a carry is added at AQ71
1638 #if defined(NEED_128)
1639     if (iszero_128 (and_128 (m, SIGN72)))
1640 #else
1641     if ((m & SIGN72) == 0)
1642 #endif
1643       {
1644         carry = 1;
1645       }
1646 #if defined(NEED_128)
1647     m = Add72b (cpup, m, construct_128 (0, 0177777777777777LL), carry, I_OFLOW, & flags1, & ovf);
1648 #else
1649     m = Add72b (cpup, m, 0177777777777777LL, carry, I_OFLOW, & flags1, & ovf);
1650 #endif
1651 
1652     // 0 -> C(AQ)28,71  (per. RJ78)
1653 #if defined(NEED_128)
1654     m = and_128 (m, lshift_128 (construct_128 (0, 0777777777400), 36));
1655 #else
1656     m &= ((word72)0777777777400 << 36);
1657 #endif
1658 
1659     // If overflow occurs, C(AQ) is shifted one place to the right and C(E) is
1660     // increased by 1.
1661     // If overflow does not occur, C(EAQ) is normalized.
1662     // All of this is done by fno, we just need to save the overflow flag
1663 
1664     bool savedovf = TST_I_OFLOW;
1665     SC_I_OFLOW(ovf);
1666     convert_to_word36 (m, & A, & Q);
1667     fno (cpup, & E, & A, & Q);
1668     SC_I_OFLOW(savedovf);
1669 
1670     * Y = setbits36_8 (A >> 8, 0, (word8) E);
1671   }
1672 
1673 /*
1674  * single precision Floating Compare ...
1675  */
1676 void fcmp (cpu_state_t * cpup) {
     /*  */
1677   // C(E) :: C(Y)0,7
1678   // C(AQ)0,27 :: C(Y)8,35
1679 
1680   // Zero: If C(EAQ) = C(Y), then ON; otherwise OFF
1681   // Neg: If C(EAQ) < C(Y), then ON; otherwise OFF
1682 
1683   // Notes: The fcmp instruction is executed as follows:
1684   // The mantissas are aligned by shifting the mantissa of the operand with
1685   // the algebraically smaller exponent to the right the number of places
1686   // equal to the difference in the two exponents.
1687   // The aligned mantissas are compared and the indicators set accordingly.
1688 
1689   CPTUR (cptUseE);
1690 #if defined(NEED_128)
1691   word72 m1= lshift_128 (construct_128 (0, cpu.rA & 0777777777400), 36);
1692 #else
1693   word72 m1 = ((word72)cpu.rA & 0777777777400LL) << 36;
1694 #endif
1695   int    e1 = SIGNEXT8_int (cpu.rE & MASK8);
1696 
1697   // 28-bit mantissa (incl sign)
1698 #if defined(NEED_128)
1699   word72 m2 = lshift_128 (construct_128 (0, getbits36_28 (cpu.CY, 8)), 44);
1700 #else
1701   word72 m2 = ((word72) getbits36_28 (cpu.CY, 8)) << 44;
1702 #endif
1703   int    e2 = SIGNEXT8_int (getbits36_8 (cpu.CY, 0));
1704 
1705   // Which exponent is smaller???
1706 
1707   L68_ (cpu.ou.cycle = ou_GOE;)
1708 #if defined(HEX_MODE)
1709   uint shift_amt = isHex(cpup) ? 4 : 1;
1710 #else
1711   uint shift_amt = 1;
1712 #endif
1713   int shift_count = -1;
1714   word1 notallzeros = 0;
1715 
1716   if (e1 == e2) {
1717     shift_count = 0;
1718   } else if (e1 < e2) {
1719     L68_ (cpu.ou.cycle = ou_GOA;)
1720     shift_count = abs (e2 - e1) * (int) shift_amt;
1721     // mantissa negative?
1722 #if defined(NEED_128)
1723     bool s = isnonzero_128 (and_128 (m1, SIGN72));
1724 #else
1725     bool s = (m1 & SIGN72) != (word72)0;
1726 #endif
1727     for(int n = 0; n < shift_count; n += 1) {
1728 #if defined(NEED_128)
1729       notallzeros |= m1.l & 1;
1730       m1 = rshift_128 (m1, 1);
1731 #else
1732       notallzeros |= m1 & 1;
1733       m1 >>= 1;
1734 #endif
1735       if (s)
1736 #if defined(NEED_128)
1737         m1 = or_128 (m1, SIGN72);
1738 #else
1739         m1 |= SIGN72;
1740 #endif
1741     }
1742 #if defined(NEED_128)
1743     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1744       m1 = construct_128 (0, 0);
1745     m1 = and_128 (m1, MASK72);
1746 #else // NEED_128
1747     if (m1 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1748       m1 = 0;
1749     m1 &= MASK72;
1750 #endif // NEED_128
1751   } else { // e2 < e1;
1752     L68_ (cpu.ou.cycle = ou_GOA;)
1753     shift_count = abs (e1 - e2) * (int) shift_amt;
1754     // mantissa negative?
1755 #if defined(NEED_128)
1756     bool s = isnonzero_128 (and_128 (m2, SIGN72));
1757 #else
1758     bool s = (m2 & SIGN72) != (word72)0;
1759 #endif
1760     for (int n = 0 ; n < shift_count ; n += 1) {
1761 #if defined(NEED_128)
1762       notallzeros |= m2.l & 1;
1763       m2 = rshift_128 (m2, 1);
1764 #else
1765       notallzeros |= m2 & 1;
1766       m2 >>= 1;
1767 #endif
1768       if (s)
1769 #if defined(NEED_128)
1770         m2 = or_128 (m2, SIGN72);
1771 #else
1772         m2 |= SIGN72;
1773 #endif
1774     } // for n
1775 #if defined(NEED_128)
1776     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1777       m2 = construct_128 (0, 0);
1778     m2 = and_128 (m2, MASK72);
1779     //e3 = e1;
1780 #else
1781     if (m2 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1782       m2 = 0;
1783     m2 &= MASK72;
1784     //e3 = e1;
1785 #endif
1786   }
1787 
1788   // need to do algebraic comparisons of mantissae
1789 #if defined(NEED_128)
1790   SC_I_ZERO (iseq_128 (m1, m2));
1791   SC_I_NEG (islt_s128 (SIGNEXT72_128(m1), SIGNEXT72_128(m2)));
1792 #else
1793   SC_I_ZERO (m1 == m2);
1794   SC_I_NEG ((int128)SIGNEXT72_128(m1) < (int128)SIGNEXT72_128(m2));
1795 #endif
1796 } // fcmp
1797 
1798 /*
1799  * single precision Floating Compare magnitude ...
1800  */
1801 void fcmg (cpu_state_t * cpup) {
     /*  */
1802   // C(E) :: C(Y)0,7
1803   // | C(AQ)0,27 | :: | C(Y)8,35 |
1804   // * Zero: If | C(EAQ)| = | C(Y) |, then ON; otherwise OFF
1805   // * Neg : If | C(EAQ)| < | C(Y) |, then ON; otherwise OFF
1806 
1807   // Notes: The fcmp instruction is executed as follows:
1808   // The mantissas are aligned by shifting the mantissa of the operand with
1809   // the algebraically smaller exponent to the right the number of places
1810   // equal to the difference in the two exponents.
1811   // The aligned mantissas are compared and the indicators set accordingly.
1812 
1813   // The fcmg instruction is identical to the fcmp instruction except that the
1814   // magnitudes of the mantissas are compared instead of the algebraic values.
1815 
1816   // ISOLTS-736 01u asserts that |0.0*2^64|<|1.0| in 28bit precision
1817   // this implies that all shifts are 72 bits long
1818   // RJ78 also comments: If the number of shifts equals or exceeds 72, the
1819   // number with the lower exponent is defined as zero.
1820 
1821   CPTUR (cptUseE);
1822   L68_ (cpu.ou.cycle = ou_GOS;)
1823 #if defined(HEX_MODE)
1824   uint shift_amt = isHex(cpup) ? 4 : 1;
1825 #else
1826   uint shift_amt = 1;
1827 #endif
1828   // C(AQ)0,27
1829 #if defined(NEED_128)
1830   word72 m1 = lshift_128 (construct_128 (0, cpu.rA & 0777777777400), 36);
1831 #else
1832   word72 m1 = ((word72)cpu.rA & 0777777777400LL) << 36;
1833 #endif
1834   int    e1 = SIGNEXT8_int (cpu.rE & MASK8);
1835 
1836   // C(Y)0,7
1837   // 28-bit mantissa (incl sign)
1838 #if defined(NEED_128)
1839   word72 m2 = lshift_128 (construct_128 (0, getbits36_28 (cpu.CY, 8)), 44);
1840 #else
1841   word72 m2 = ((word72) getbits36_28 (cpu.CY, 8)) << 44;
1842 #endif
1843   int    e2 = SIGNEXT8_int (getbits36_8 (cpu.CY, 0));
1844 
1845   // Which exponent is smaller???
1846 
1847   L68_ (cpu.ou.cycle = ou_GOE;)
1848   int shift_count = -1;
1849   word1 notallzeros = 0;
1850 
1851   if (e1 == e2) {
1852     shift_count = 0;
1853   } else if (e1 < e2) {
1854     L68_ (cpu.ou.cycle = ou_GOA;)
1855     shift_count = abs (e2 - e1) * (int) shift_amt;
1856 #if defined(NEED_128)
1857     bool s = isnonzero_128 (and_128 (m1, SIGN72));
1858 #else
1859     bool s = (m1 & SIGN72) != (word72)0;    ///< save sign bit
1860 #endif
1861     for (int n = 0; n < shift_count; n += 1) {
1862 #if defined(NEED_128)
1863       notallzeros |= m1.l & 1;
1864       m1 = rshift_128 (m1, 1);
1865 #else
1866       notallzeros |= m1 & 1;
1867       m1 >>= 1;
1868 #endif
1869     if (s)
1870 #if defined(NEED_128)
1871       m1 = or_128 (m1, SIGN72);
1872 #else
1873       m1 |= SIGN72;
1874 #endif
1875     }
1876 
1877 #if defined(NEED_128)
1878     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1879       m1 = construct_128 (0, 0);
1880     m1 = and_128 (m1, MASK72);
1881 #else
1882     if (m1 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1883       m1 = 0;
1884     m1 &= MASK72;
1885 #endif
1886   } else { // e2 < e1;
1887     L68_ (cpu.ou.cycle = ou_GOA;)
1888     shift_count = abs(e1 - e2) * (int) shift_amt;
1889 #if defined(NEED_128)
1890     bool s = isnonzero_128 (and_128 (m2, SIGN72));
1891 #else
1892     bool s = (m2 & SIGN72) != (word72)0;    ///< save sign bit
1893 #endif
1894     for (int n = 0; n < shift_count; n += 1) {
1895 #if defined(NEED_128)
1896       notallzeros |= m2.l & 1;
1897       m2 = rshift_128 (m2, 1);
1898 #else
1899       notallzeros |= m2 & 1;
1900       m2 >>= 1;
1901 #endif
1902       if (s)
1903 #if defined(NEED_128)
1904         m2 = or_128 (m2, SIGN72);
1905 #else
1906         m2 |= SIGN72;
1907 #endif
1908     }
1909 #if defined(NEED_128)
1910     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1911       m2 = construct_128 (0, 0);
1912     m2 = and_128 (m2, MASK72);
1913 #else
1914     if (m2 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
1915       m2 = 0;
1916     m2 &= MASK72;
1917 #endif
1918   }
1919 
1920   SC_I_ZERO (ISEQ_128 (m1, m2));
1921 #if defined(NEED_128)
1922   int128 sm1 = SIGNEXT72_128 (m1);
1923   if (islt_s128 (sm1, construct_s128 (0, 0)))
1924     sm1 = negate_s128 (sm1);
1925   int128 sm2 = SIGNEXT72_128 (m2);
1926   if (islt_s128 (sm2, construct_s128 (0, 0)))
1927     sm2 = negate_s128 (sm2);
1928   SC_I_NEG (islt_s128 (sm1, sm2));
1929 #else
1930   int128 sm1 = SIGNEXT72_128 (m1);
1931   if (sm1 < 0)
1932     sm1 = - sm1;
1933   int128 sm2 = SIGNEXT72_128 (m2);
1934   if (sm2 < 0)
1935     sm2 = - sm2;
1936   SC_I_NEG (sm1 < sm2);
1937 #endif
1938 }
1939 
1940 /*
1941  * double-precision arithmetic routines ....
1942  */
1943 
1944 
1945 
1946 
1947 
1948 
1949 
1950 
1951 
1952 
1953 
1954 
1955 
1956 
1957 
1958 
1959 
1960 
1961 
1962 
1963 
1964 /*!
1965  * unnormalized floating double-precision add
1966  */
1967 void dufa (cpu_state_t * cpup, bool subtract, bool normalize) {
     /*  */
1968   // Except for the precision of the mantissa of the operand from main
1969   // memory, the dufa instruction is identical to the ufa instruction.
1970 
1971   // C(EAQ) + C(Y) → C(EAQ)
1972   // The ufa instruction is executed as follows:
1973   // The mantissas are aligned by shifting the mantissa of the operand having
1974   // the algebraically smaller exponent to the right the number of places
1975   // equal to the absolute value of the difference in the two exponents. Bits
1976   // shifted beyond the bit position equivalent to AQ71 are lost.
1977   // The algebraically larger exponent replaces C(E). The sum of the
1978   // mantissas replaces C(AQ).
1979   // If an overflow occurs during addition, then;
1980   // *  C(AQ) are shifted one place to the right.
1981   // *  C(AQ)0 is inverted to restore the sign.
1982   // *  C(E) is increased by one.
1983 
1984   // The dufs instruction is identical to the dufa instruction with the
1985   // exception that the twos complement of the mantissa of the operand from
1986   // main memory (op2) is used.
1987 
1988   CPTUR (cptUseE);
1989   L68_ (cpu.ou.cycle = ou_GOS;)
1990 #if defined(HEX_MODE)
1991   uint shift_amt = isHex(cpup) ? 4 : 1;
1992 #else
1993   uint shift_amt = 1;
1994 #endif
1995 
1996   word72 m1 = convert_to_word72 (cpu.rA, cpu.rQ);
1997   int e1 = SIGNEXT8_int (cpu.rE & MASK8);
1998 
1999   // 64-bit mantissa (incl sign)
2000 #if defined(NEED_128)
2001   word72 m2 = lshift_128 (construct_128 (0, (uint64_t) getbits36_28 (cpu.Ypair[0], 8)), 44u); // 28-bit mantissa (incl sign)
2002          m2 = or_128 (m2, lshift_128 (construct_128 (0, cpu.Ypair[1]), 8u));
2003 #else
2004   word72 m2 = ((word72) getbits36_28 (cpu.Ypair[0], 8)) << 44;
2005          m2 |= (word72) cpu.Ypair[1] << 8;
2006 #endif
2007 
2008   int e2 = SIGNEXT8_int (getbits36_8 (cpu.Ypair[0], 0));
2009 
2010   // see ufs
2011   int m2zero = 0;
2012   if (subtract) {
2013 #if defined(NEED_128)
2014     if (iszero_128 (m2))
2015       m2zero = 1;
2016     if (iseq_128 (m2, SIGN72)) {
2017 # if defined(HEX_MODE)
2018       m2 = rshift_128 (m2, shift_amt);
2019       e2 += 1;
2020 # else
2021       m2 = rshift_128 (m2, 1);
2022       e2 += 1;
2023 # endif
2024     } else {
2025       m2 = and_128 (negate_128 (m2), MASK72);
2026     }
2027 #else
2028     if (m2 == 0)
2029       m2zero = 1;
2030     if (m2 == SIGN72) {
2031 # if defined(HEX_MODE)
2032       m2 >>= shift_amt;
2033       e2 += 1;
2034 # else
2035       m2 >>= 1;
2036       e2 += 1;
2037 # endif
2038     } else {
2039       // m2 = (-m2) & MASK72;
2040       // ubsan
2041       m2 = ((word72) (- (word72s) m2)) & MASK72;
2042     }
2043 #endif
2044   }
2045 
2046   int e3 = -1;
2047 
2048   // Which exponent is smaller?
2049 
2050   L68_ (cpu.ou.cycle |= ou_GOE;)
2051   int shift_count = -1;
2052   word1 notallzeros = 0;
2053 
2054   if (e1 == e2) {
2055     shift_count = 0;
2056     e3 = e1;
2057   } else if (e1 < e2) {
2058     L68_ (cpu.ou.cycle |= ou_GOA;)
2059     shift_count = abs(e2 - e1) * (int) shift_amt;
2060     // mantissa negative?
2061 #if defined(NEED_128)
2062     bool s = isnonzero_128 (and_128 (m1, SIGN72));
2063 #else
2064     bool s = (m1 & SIGN72) != (word72)0;
2065 #endif
2066     for (int n = 0; n < shift_count; n += 1) {
2067 #if defined(NEED_128)
2068       notallzeros |= m1.l & 1;
2069       m1 = rshift_128 (m1, 1);
2070 #else
2071       notallzeros |= m1 & 1;
2072       m1 >>= 1;
2073 #endif
2074       if (s)
2075 #if defined(NEED_128)
2076         m1 = or_128 (m1, SIGN72);
2077 #else
2078         m1 |= SIGN72;
2079 #endif
2080     }
2081 #if defined(NEED_128)
2082     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
2083       m1 = construct_128 (0, 0);
2084     m1 = and_128 (m1, MASK72);
2085 #else
2086     if (m1 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
2087       m1 = 0;
2088     m1 &= MASK72;
2089 #endif
2090     e3 = e2;
2091   } else { // e2 < e1;
2092     L68_ (cpu.ou.cycle |= ou_GOA;)
2093     shift_count = abs(e1 - e2) * (int) shift_amt;
2094 #if defined(NEED_128)
2095     bool s = isnonzero_128 (and_128 (m2, SIGN72));
2096 #else
2097     bool s = (m2 & SIGN72) != (word72)0;    ///< save sign bit
2098 #endif
2099     for (int n = 0; n < shift_count; n += 1) {
2100 #if defined(NEED_128)
2101       notallzeros |= m2.l & 1;
2102       m2 = rshift_128 (m2, 1);
2103 #else
2104       notallzeros |= m2 & 1;
2105       m2 >>= 1;
2106 #endif
2107       if (s)
2108 #if defined(NEED_128)
2109         m2 = or_128 (m2, SIGN72);
2110 #else
2111         m2 |= SIGN72;
2112 #endif
2113     }
2114 #if defined(NEED_128)
2115     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
2116       m2 = construct_128 (0, 0);
2117     m2 = and_128 (m2, MASK72);
2118 #else
2119     if (m2 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
2120       m2 = 0;
2121     m2 &= MASK72;
2122 #endif
2123     e3 = e1;
2124   }
2125 
2126   bool ovf;
2127   word72 m3 = Add72b (cpup, m1, m2, 0, I_CARRY, & cpu.cu.IR, & ovf);
2128   if (m2zero)
2129     SET_I_CARRY;
2130 
2131   if (ovf) {
2132 #if defined(HEX_MODE)
2133 //    word72 signbit = m3 & sign_msk;
2134 //    m3 >>= shift_amt;
2135 //    m3 = (m3 & MASK71) | signbit;
2136 //    m3 ^= SIGN72; // C(AQ)0 is inverted to restore the sign
2137 //    e3 += 1;
2138     // save the sign bit
2139 # if defined(NEED_128)
2140     bool s = isnonzero_128 (and_128 (m3, SIGN72));
2141 # else
2142     bool s = (m3 & SIGN72) != (word72)0;
2143 # endif
2144     if (isHex (cpup)) {
2145 # if defined(NEED_128)
2146       m3 = rshift_128 (m3, shift_amt); // renormalize the mantissa
2147       if (s)
2148         // Sign is set, number should be positive; clear the sign bit and the 3 MSBs
2149         m3 = and_128 (m3, MASK68);
2150       else
2151         // Sign is clr, number should be negative; set the sign bit and the 3 MSBs
2152         m3 = or_128 (m3, HEX_SIGN);
2153 # else
2154       m3 >>= shift_amt; // renormalize the mantissa
2155       if (s)
2156         // Sign is set, number should be positive; clear the sign bit and the 3 MSBs
2157         m3 &= MASK68;
2158       else
2159         // Sign is clr, number should be negative; set the sign bit and the 3 MSBs
2160         m3 |=  HEX_SIGN;
2161 # endif
2162     } else {
2163 # if defined(NEED_128)
2164       word72 signbit = and_128 (m3, SIGN72);
2165       m3 = rshift_128 (m3, 1); // renormalize the mantissa
2166       m3 = or_128 (and_128 (m3, MASK71), signbit);
2167       m3 = xor_128 (m3, SIGN72); // C(AQ)0 is inverted to restore the sign
2168 # else
2169       word72 signbit = m3 & SIGN72;
2170       m3 >>= 1;
2171       m3 = (m3 & MASK71) | signbit;
2172       m3 ^= SIGN72; // C(AQ)0 is inverted to restore the sign
2173 # endif
2174     }
2175     e3 += 1;
2176 #else
2177 # if defined(NEED_128)
2178     word72 signbit = and_128 (m3, SIGN72);
2179     m3 = rshift_128 (m3, 1); // renormalize the mantissa
2180     m3 = or_128 (and_128 (m3, MASK71), signbit);
2181     m3 = xor_128 (m3, SIGN72); // C(AQ)0 is inverted to restore the sign
2182 # else
2183     word72 signbit = m3 & SIGN72;
2184     m3 >>= 1;
2185     m3 = (m3 & MASK71) | signbit;
2186     m3 ^= SIGN72; // C(AQ)0 is inverted to restore the sign
2187 # endif
2188     e3 += 1;
2189 #endif
2190   }
2191 
2192   convert_to_word36 (m3, & cpu.rA, & cpu.rQ);
2193 #if defined(TESTING)
2194   HDBGRegAW ("dufa");
2195 #endif
2196   cpu.rE = e3 & 0377;
2197 
2198   SC_I_NEG (cpu.rA & SIGN36); // Do this here instead of in Add72b because
2199                                 // of ovf handling above
2200   if (cpu.rA == 0 && cpu.rQ == 0) {
2201     SET_I_ZERO;
2202     cpu.rE = 0200U; /*-128*/
2203   } else {
2204     CLR_I_ZERO;
2205   }
2206 
2207   if (normalize) {
2208     fno_ext (cpup, & e3, & cpu.rE, & cpu.rA, & cpu.rQ);
2209   }
2210 
2211   // EOFL: If exponent is greater than +127, then ON
2212   if (e3 > 127) {
2213     SET_I_EOFL;
2214     if (tstOVFfault (cpup))
2215       dlyDoFault (FAULT_OFL, fst_zero, "dufa exp overflow fault");
2216   }
2217 
2218   // EUFL: If exponent is less than -128, then ON
2219   if (e3 < -128) {
2220     SET_I_EUFL;
2221     if (tstOVFfault (cpup))
2222         dlyDoFault (FAULT_OFL, fst_zero, "dufa exp underflow fault");
2223   }
2224 } // dufa
2225 
2226 
2227 
2228 
2229 
2230 
2231 
2232 
2233 
2234 
2235 
2236 
2237 
2238 
2239 
2240 
2241 
2242 
2243 
2244 
2245 
2246 
2247 
2248 
2249 
2250 
2251 
2252 
2253 
2254 
2255 
2256 
2257 
2258 
2259 
2260 
2261 
2262 
2263 
2264 
2265 
2266 
2267 
2268 
2269 
2270 
2271 
2272 
2273 
2274 
2275 
2276 
2277 
2278 
2279 
2280 
2281 
2282 
2283 
2284 
2285 
2286 
2287 
2288 /*
2289  * double-precision Unnormalized Floating Multiply ...
2290  */
2291 void dufm (cpu_state_t * cpup, bool normalize) {
     /*  */
2292   // Except for the precision of the mantissa of the operand from main memory,
2293   //the dufm instruction is identical to the ufm instruction.
2294 
2295   // The ufm instruction is executed as follows:
2296   //      C(E) + C(Y)0,7 → C(E)
2297   //      ( C(AQ) × C(Y)8,35 )0,71 → C(AQ)
2298 
2299   // * Zero: If C(AQ) = 0, then ON; otherwise OFF
2300   // * Neg: If C(AQ)0 = 1, then ON; otherwise OFF
2301   // * Exp Ovr: If exponent is greater than +127, then ON
2302   // * Exp Undr: If exponent is less than -128, then ON
2303 
2304   CPTUR (cptUseE);
2305   L68_ (cpu.ou.cycle |= ou_GOS;)
2306   word72 m1 = convert_to_word72 (cpu.rA, cpu.rQ);
2307   int    e1 = SIGNEXT8_int (cpu . rE & MASK8);
2308 
2309 #if !defined(NEED_128)
2310   sim_debug (DBG_TRACEEXT, & cpu_dev, "dufm e1 %d %03o m1 %012"PRIo64" %012"PRIo64"\n",
2311              e1, e1, (word36) (m1 >> 36) & MASK36, (word36) m1 & MASK36);
2312 #endif
2313    // 64-bit mantissa (incl sign)
2314 #if defined(NEED_128)
2315   word72 m2 = lshift_128 (construct_128 (0, (uint64_t) getbits36_28 (cpu.Ypair[0], 8)), 44u); // 28-bit mantissa (incl sign)
2316          m2 = or_128 (m2, lshift_128 (construct_128 (0, cpu.Ypair[1]), 8u));
2317 #else
2318   word72 m2 = ((word72) getbits36_28 (cpu.Ypair[0], 8)) << 44;
2319          m2 |= (word72) cpu.Ypair[1] << 8;
2320 #endif
2321 
2322   // 8-bit signed integer (incl sign)
2323   int    e2 = SIGNEXT8_int (getbits36_8 (cpu.Ypair[0], 0));
2324 
2325 #if !defined(NEED_128)
2326   sim_debug (DBG_TRACEEXT, & cpu_dev, "dufm e2 %d %03o m2 %012"PRIo64" %012"PRIo64"\n",
2327              e2, e2, (word36) (m2 >> 36) & MASK36, (word36) m2 & MASK36);
2328 #endif
2329 
2330   if (ISZERO_128 (m1) || ISZERO_128 (m2)) {
2331     SET_I_ZERO;
2332     CLR_I_NEG;
2333 
2334     cpu.rE = 0200U; /*-128*/
2335     cpu.rA = 0;
2336 #if defined(TESTING)
2337     HDBGRegAW ("dufm");
2338 #endif
2339     cpu.rQ = 0;
2340 
2341     return; // normalized 0
2342   }
2343 
2344   int e3 = e1 + e2;
2345 
2346 
2347 
2348 
2349 
2350 
2351 
2352 
2353 
2354 
2355 
2356 
2357 
2358 
2359   // RJ78: This multiplication is executed in the following way:
2360   // C(E) + C(Y)(0-7) -> C(E)
2361   // C(AQ) * C(Y-pair)(8-71) results in a 134-bit product plus sign. This sign plus the
2362   // leading 71 bits are loaded into the AQ.
2363 
2364   // do a 128x64 signed multiplication
2365 
2366   // fast signed multiplication algorithm without 2's complements
2367   // passes ISOLTS-745 08
2368 
2369 #if defined(NEED_128)
2370   // shift the CY mantissa
2371   int128 m2s = rshift_s128 (SIGNEXT72_128(m2), 8);
2372 
2373   // do a 128x64 signed multiplication
2374   int128 m1l = and_s128 (cast_s128 (m1), construct_128 (0, MASK64));
2375   int128 m1h = rshift_s128 (SIGNEXT72_128(m1), 64);
2376   int128 m3h = multiply_s128 (m1h, m2s); // hi partial product
2377   int128 m3l = multiply_s128 (m1l, m2s); // lo partial product
2378 
2379   // realign to 72bits
2380   m3l = rshift_s128 (m3l, 63);
2381   m3h = lshift_s128 (m3h, 1); // m3h is hi by 64, align it for addition. The result is 135 bits so this cannot overflow.
2382   word72 m3a = and_128 (add_128 (cast_128 (m3h), cast_128 (m3l)), MASK72);
2383 #else
2384   // shift the CY mantissa
2385   int128 m2s = SIGNEXT72_128(m2) >> 8;
2386 
2387   // do a 128x64 signed multiplication
2388   int128 m1l = m1 & (((uint128)1<<64)-1);
2389   int128 m1h = SIGNEXT72_128(m1) >> 64;
2390   int128 m3h = m1h * m2s; // hi partial product
2391   int128 m3l = m1l * m2s; // lo partial product
2392 
2393   // realign to 72bits
2394   m3l >>= 63;
2395   // m3h <<= 1; // m3h is hi by 64, align it for addition. The result is 135 bits so this cannot overflow.
2396   // ubsan
2397   m3h = (int128) (((uint128) m3h) << 1); // m3h is hi by 64, align it for addition. The result is 135 bits so this cannot overflow.
2398   word72 m3a = ((word72) (m3h+m3l)) & MASK72;
2399 #endif
2400 
2401   // A normalization is performed only in the case of both factor mantissas being 100...0
2402   // which is the twos complement approximation to the decimal value -1.0.
2403   if (ISEQ_128 (m1, SIGN72) && ISEQ_128 (m2, SIGN72)) {
2404     if (e3 == 127) {
2405       SET_I_EOFL;
2406       if (tstOVFfault (cpup))
2407           dlyDoFault (FAULT_OFL, fst_zero, "dufm exp overflow fault");
2408     }
2409 #if defined(NEED_128)
2410     m3a = rshift_128 (m3a, 1);
2411 #else
2412     m3a >>= 1;
2413 #endif
2414     e3 += 1;
2415   }
2416 
2417   convert_to_word36 (m3a, & cpu.rA, & cpu.rQ);
2418 #if defined(TESTING)
2419   HDBGRegAW ("dufm");
2420 #endif
2421   cpu.rE = (word8) e3 & MASK8;
2422 
2423   SC_I_NEG (cpu.rA & SIGN36);
2424 
2425   if (cpu.rA == 0 && cpu.rQ == 0) {
2426     SET_I_ZERO;
2427     cpu . rE = 0200U; /*-128*/
2428   } else {
2429     CLR_I_ZERO;
2430   }
2431 
2432   if (normalize) {
2433     fno_ext (cpup, & e3, & cpu.rE, & cpu.rA, & cpu.rQ);
2434   }
2435 
2436   if (e3 >  127) {
2437     SET_I_EOFL;
2438     if (tstOVFfault (cpup))
2439       dlyDoFault (FAULT_OFL, fst_zero, "dufm exp overflow fault");
2440   }
2441   if (e3 < -128) {
2442     SET_I_EUFL;
2443     if (tstOVFfault (cpup))
2444       dlyDoFault (FAULT_OFL, fst_zero, "dufm exp underflow fault");
2445   }
2446 }
2447 
2448 /*
2449  * floating divide ...
2450  */
2451 // CANFAULT
2452 static void dfdvX (cpu_state_t * cpup, bool bInvert) {
     /*  */
2453   // C(EAQ) / C (Y) → C(EA)
2454   // C(Y) / C(EAQ) → C(EA) (Inverted)
2455 
2456   // 00...0 → C(Q)
2457 
2458   // The dfdv instruction is executed as follows:
2459   // The dividend mantissa C(AQ) is shifted right and the dividend exponent
2460   // C(E) increased accordingly until
2461   //    | C(AQ)0,63 | < | C(Y-pair)8,71 |
2462   //    C(E) - C(Y-pair)0,7 → C(E)
2463   //    C(AQ) / C(Y-pair)8,71 → C(AQ)0,63 00...0 → C(Q)64,71
2464 
2465   CPTUR (cptUseE);
2466   L68_ (cpu.ou.cycle |= ou_GOS;)
2467 #if defined(HEX_MODE)
2468   uint shift_amt = isHex(cpup) ? 4 : 1;
2469 #else
2470   uint shift_amt = 1;
2471 #endif
2472   word72 m1;
2473   int    e1;
2474 
2475   word72 m2;
2476   int    e2;
2477 
2478   bool roundovf = 0;
2479 
2480   if (! bInvert) {
2481     m1 = convert_to_word72 (cpu.rA, cpu.rQ);
2482     e1 = SIGNEXT8_int (cpu.rE & MASK8);
2483 
2484     // 64-bit mantissa (incl sign)
2485 #if defined(NEED_128)
2486     m2 = lshift_128 (construct_128 (0, (uint64_t) getbits36_28 (cpu.Ypair[0], 8)), 44u); // 28-bit mantissa (incl sign)
2487     m2 = or_128 (m2, lshift_128 (construct_128 (0, cpu.Ypair[1]), 8u));
2488 #else
2489     m2 = ((word72) getbits36_28 (cpu.Ypair[0], 8)) << 44;
2490     m2 |= (word72) cpu.Ypair[1] << 8;
2491 #endif
2492 
2493     e2 = SIGNEXT8_int (getbits36_8 (cpu.Ypair[0], 0));
2494   } else { // invert
2495     m2 = convert_to_word72 (cpu.rA, cpu.rQ);
2496     e2 = SIGNEXT8_int (cpu.rE & MASK8);
2497 
2498     // round divisor per RJ78
2499     // If AQ(64-71) is not = 0 and A(0) = 0, a 1 is added to AQ(63). Zero is moved to
2500     // AQ(64-71), unconditionally. AQ(0-63) is then used as the divisor mantissa.
2501     // ISOLTS-745 10b
2502 #if defined(NEED_128)
2503     if ((iszero_128 (and_128 (m2, SIGN72))) && m2.l & 0377)
2504 #else
2505     if (!(m2 & SIGN72) && m2 & 0377)
2506 #endif
2507     {
2508 #if defined(NEED_128)
2509       m2 = add_128 (m2, construct_128 (0, 0400));
2510 #else
2511       m2 += 0400;
2512 #endif
2513       // ISOLTS-745 10e asserts that an overflowing addition of 400
2514       //   to 377777777777 7777777774xx does not shift the quotient (nor divisor)
2515       // I surmise that the divisor is taken as unsigned 64 bits in this case
2516       roundovf = 1;
2517     }
2518 #if defined(NEED_128)
2519     putbits72 (& m2, 64, 8, construct_128 (0, 0));
2520 #else
2521     putbits72 (& m2, 64, 8, 0);
2522 #endif
2523 
2524     // 64-bit mantissa (incl sign)
2525 #if defined(NEED_128)
2526     m1 = lshift_128 (construct_128 (0, (uint64_t) getbits36_28 (cpu.Ypair[0], 8)), 44u); // 28-bit mantissa (incl sign)
2527     m1 = or_128 (m1, lshift_128 (construct_128 (0, cpu.Ypair[1]), 8u));
2528 #else
2529     m1 = ((word72) getbits36_28 (cpu.Ypair[0], 8)) << 44;
2530     m1 |= (word72) cpu.Ypair[1] << 8;
2531 #endif
2532 
2533     e1 = SIGNEXT8_int (getbits36_8 (cpu.Ypair[0], 0));
2534   } // invert
2535 
2536   if (ISZERO_128 (m1)) {
2537     SET_I_ZERO;
2538     CLR_I_NEG;
2539 
2540     cpu.rE = 0200U; /*-128*/
2541     cpu.rA = 0;
2542 #if defined(TESTING)
2543     HDBGRegAW ("dfdvX");
2544 #endif
2545     cpu.rQ = 0;
2546 
2547     return; // normalized 0
2548   }
2549 
2550     // make everything positive, but save sign info for later....
2551   int sign = 1;
2552 #if defined(NEED_128)
2553   if (isnonzero_128 (and_128 (m1, SIGN72))) {
2554     SET_I_NEG; // in case of divide fault
2555     if (iseq_128 (m1, SIGN72)) {
2556       m1 = rshift_128 (m1, shift_amt);
2557       e1 += 1;
2558     } else {
2559       m1 = and_128 (negate_128 (m1), MASK72);
2560     }
2561     sign = -sign;
2562   } else {
2563     CLR_I_NEG; // in case of divide fault
2564   }
2565 
2566   if ((isnonzero_128 (and_128 (m2, SIGN72))) && !roundovf) {
2567     if (iseq_128 (m2, SIGN72)) {
2568       m2 = rshift_128 (m2, shift_amt);
2569       e2 += 1;
2570     } else {
2571       m2 = and_128 (negate_128 (m2), MASK72);
2572     }
2573     sign = -sign;
2574   }
2575 #else
2576   if (m1 & SIGN72) {
2577     SET_I_NEG; // in case of divide fault
2578     if (m1 == SIGN72) {
2579       m1 >>= shift_amt;
2580       e1 += 1;
2581     } else {
2582       m1 = (~m1 + 1) & MASK72;
2583     }
2584     sign = -sign;
2585   } else {
2586     CLR_I_NEG; // in case of divide fault
2587   }
2588 
2589   if ((m2 & SIGN72) && !roundovf) {
2590     if (m2 == SIGN72) {
2591       m2 >>= shift_amt;
2592       e2 += 1;
2593     } else {
2594       m2 = (~m2 + 1) & MASK72;
2595     }
2596     sign = -sign;
2597   }
2598 #endif
2599 
2600   if (ISZERO_128 (m2)) {
2601     // NB: If C(Y-pair)8,71 == 0 then the alignment loop will never exit! That's why it been moved before the alignment
2602 
2603     SET_I_ZERO;
2604     // NEG already set
2605 
2606     // FDV: If the divisor mantissa C(Y-pair)8,71 is zero after alignment (HWR: why after?), the division does
2607     // not take place. Instead, a divide check fault occurs, C(AQ) contains the dividend magnitude, and
2608     // the negative indicator reflects the dividend sign.
2609     // FDI: If the divisor mantissa C(AQ) is zero, the division does not take place.
2610     // Instead, a divide check fault occurs and all the registers remain unchanged.
2611     if (! bInvert) {
2612       convert_to_word36 (m1, & cpu.rA, & cpu.rQ);
2613 #if defined(TESTING)
2614       HDBGRegAW ("dfdvX");
2615 #endif
2616     }
2617     doFault (FAULT_DIV, fst_zero, "DFDV: divide check fault");
2618   } // m2 == 0
2619 
2620   L68_ (cpu.ou.cycle |= ou_GOA;)
2621 #if defined(NEED_128)
2622   while (isge_128 (m1, m2)) {
2623     m1 = rshift_128 (m1, shift_amt);
2624     e1 += 1;
2625   }
2626 #else
2627   while (m1 >= m2) {
2628     m1 >>= shift_amt;
2629     e1 += 1;
2630   }
2631 #endif
2632   int e3 = e1 - e2;
2633   if (e3 > 127) {
2634     SET_I_EOFL;
2635     if (tstOVFfault (cpup))
2636       dlyDoFault (FAULT_OFL, fst_zero, "dfdvX exp overflow fault");
2637    }
2638   if (e3 < -128) {
2639     SET_I_EUFL;
2640     if (tstOVFfault (cpup))
2641       dlyDoFault (FAULT_OFL, fst_zero, "dfdvX exp underflow fault");
2642   }
2643 
2644   L68_ (cpu.ou.cycle |= ou_GD1;)
2645 
2646   // We need 63 bits quotient + sign. Divisor is at most 64 bits.
2647   // Do a 127 by 64 fractional divide
2648   // lo 8bits are always zero
2649 #if defined(NEED_128)
2650   word72 m3 = divide_128 (lshift_128 (m1, 63-8), rshift_128 (m2, 8), NULL);
2651 #else
2652   word72 m3 = ((uint128)m1 << (63-8)) / ((uint128)m2 >> 8);
2653 #endif
2654   L68_ (cpu.ou.cycle |= ou_GD2;)
2655 
2656 #if defined(NEED_128)
2657   m3 = lshift_128 (m3, 8);  // convert back to float
2658   if (sign == -1)
2659     m3 = and_128 (negate_128 (m3), MASK72);
2660 #else
2661   m3 <<= 8;  // convert back to float
2662   if (sign == -1)
2663     m3 = (~m3 + 1) & MASK72;
2664 #endif
2665 
2666   convert_to_word36 (m3, & cpu.rA, & cpu.rQ);
2667 #if defined(TESTING)
2668   HDBGRegAW ("dfdvX");
2669 #endif
2670   cpu.rE = (word8) e3 & MASK8;
2671 
2672   SC_I_ZERO (cpu.rA == 0 && cpu . rQ == 0);
2673   SC_I_NEG (cpu.rA & SIGN36);
2674 
2675   if (cpu.rA == 0 && cpu.rQ == 0)    // set to normalized 0
2676       cpu.rE = 0200U; /*-128*/
2677 } // dfdvX
2678 
2679 // CANFAULT
2680 void dfdv (cpu_state_t * cpup) {
     /*  */
2681   dfdvX (cpup, false);    // no inversion
2682 }
2683 
2684 // CANFAULT
2685 void dfdi (cpu_state_t * cpup) {
     /*  */
2686   dfdvX (cpup, true); // invert
2687 }
2688 
2689 //#define DVF_HWR
2690 //#define DVF_FRACTIONAL
2691 #define DVF_CAC
2692 // CANFAULT
2693 void dvf (cpu_state_t * cpup)
     /*  */
2694 {
2695 //#if defined(DBGCAC)
2696 //sim_printf ("DVF %"PRId64" %06o:%06o\n", cpu.cycleCnt, cpu.PPR.PSR, cpu.PPR.IC);
2697 //#endif
2698     // C(AQ) / (Y)
2699     //  fractional quotient → C(A)
2700     //  fractional remainder → C(Q)
2701 
2702     // A 71-bit fractional dividend (including sign) is divided by a 36-bit
2703     // fractional divisor yielding a 36-bit fractional quotient (including
2704     // sign) and a 36-bit fractional remainder (including sign). C(AQ)71 is
2705     // ignored; bit position 35 of the remainder corresponds to bit position
2706     // 70 of the dividend. The remainder sign is equal to the dividend sign
2707     // unless the remainder is zero.
2708     //
2709     // If | dividend | >= | divisor | or if the divisor = 0, division does
2710     // not take place. Instead, a divide check fault occurs, C(AQ) contains
2711     // the dividend magnitude in absolute, and the negative indicator
2712     // reflects the dividend sign.
2713 
2714 // HWR code
2715 //HWR--#if defined(DVF_HWR)
2716 //HWR--    // m1 Dividend
2717 //HWR--    // m2 divisor
2718 //HWR--
2719 //HWR--    word72 m1 = SIGNEXT36_72((cpu . rA << 36) | (cpu . rQ & 0777777777776LLU));
2720 //HWR--    word72 m2 = SIGNEXT36_72(cpu.CY);
2721 //HWR--
2722 //HWR--//sim_debug (DBG_CAC, & cpu_dev, "[%"PRId64"]\n", cpu.cycleCnt);
2723 //HWR--//sim_debug (DBG_CAC, & cpu_dev, "m1 "); print_int128 (m1); sim_printf ("\n");
2724 //HWR--//sim_debug (DBG_CAC, & cpu_dev, "-----------------\n");
2725 //HWR--//sim_debug (DBG_CAC, & cpu_dev, "m2 "); print_int128 (m2); sim_printf ("\n");
2726 //HWR--
2727 //HWR--    if (m2 == 0)
2728 //HWR--    {
2729 //HWR--        // XXX check flags
2730 //HWR--        SET_I_ZERO;
2731 //HWR--        SC_I_NEG (cpu . rA & SIGN36);
2732 //HWR--
2733 //HWR--        cpu . rA = 0;
2734 //HWR--        cpu . rQ = 0;
2735 //HWR--
2736 //HWR--        return;
2737 //HWR--    }
2738 //HWR--
2739 //HWR--    // make everything positive, but save sign info for later....
2740 //HWR--    int sign = 1;
2741 //HWR--    int dividendSign = 1;
2742 //HWR--    if (m1 & SIGN72)
2743 //HWR--    {
2744 //HWR--        m1 = (~m1 + 1);
2745 //HWR--        sign = -sign;
2746 //HWR--        dividendSign = -1;
2747 //HWR--    }
2748 //HWR--
2749 //HWR--    if (m2 & SIGN72)
2750 //HWR--    {
2751 //HWR--        m2 = (~m2 + 1);
2752 //HWR--        sign = -sign;
2753 //HWR--    }
2754 //HWR--
2755 //HWR--    if (m1 >= m2 || m2 == 0)
2756 //HWR--    {
2757 //HWR--        //cpu . rA = m1;
2758 //HWR--        cpu . rA = (m1 >> 36) & MASK36;
2759 //HWR--        cpu . rQ = m1 & 0777777777776LLU;
2760 //HWR--
2761 //HWR--        SET_I_ZERO;
2762 //HWR--        SC_I_NEG (cpu . rA & SIGN36);
2763 //HWR--
2764 //HWR--        doFault(FAULT_DIV, fst_zero, "DVF: divide check fault");
2765 //HWR--    }
2766 //HWR--
2767 //HWR--    uint128 dividend = (uint128)m1 << 63;
2768 //HWR--    uint128 divisor = (uint128)m2;
2769 //HWR--
2770 //HWR--    //uint128 m3  = ((uint128)m1 << 63) / (uint128)m2;
2771 //HWR--    //uint128 m3r = ((uint128)m1 << 63) % (uint128)m2;
2772 //HWR--    int128 m3  = (int128)(dividend / divisor);
2773 //HWR--    int128 m3r = (int128)(dividend % divisor);
2774 //HWR--
2775 //HWR--    if (sign == -1)
2776 //HWR--        m3 = -m3;   //(~m3 + 1);
2777 //HWR--
2778 //HWR--    if (dividendSign == -1) // The remainder sign is equal to the dividend sign unless the remainder is zero.
2779 //HWR--        m3r = -m3r; //(~m3r + 1);
2780 //HWR--
2781 //HWR--    cpu . rA = (m3 >> 64) & MASK36;
2782 //HWR--    cpu . rQ = m3r & MASK36;   //01777777777LL;
2783 //HWR--#endif
2784 
2785 // canonical code
2786 #if defined(DVF_FRACTIONAL)
2787 // See: https://www.ece.ucsb.edu/~parhami/pres_folder/f31-book-arith-pres-pt4.pdf
2788 // Slide 10: Sequential Algorithm
2789 
2790     // dividend format
2791     // 0  1     70 71
2792     // s  dividend x
2793     //  C(AQ)
2794 
2795   L68_ (cpu.ou.cycle |= ou_GD1;)
2796   int sign = 1;
2797   bool dividendNegative = (getbits36_1 (cpu.rA, 0) != 0);
2798   bool divisorNegative = (getbits36_1 (cpu.CY, 0) != 0);
2799 
2800   // Get the 70 bits of the dividend (72 bits less the sign bit and the
2801   // ignored bit 71.
2802 
2803   // dividend format:   . d(0) ...  d(69)
2804 
2805   uint128 zFrac = (((uint128) (cpu.rA & MASK35)) << 35) | ((cpu.rQ >> 1) & MASK35);
2806   if (dividendNegative) {
2807     zFrac = ~zFrac + 1;
2808     sign = - sign;
2809   }
2810   zFrac &= MASK70;
2811   //char buf [128];
2812   //buf [0] = 0;
2813   //print_int128o (zFrac, buf);
2814   //sim_printf ("zFrac %s \r\n", buf);
2815 
2816   // Get the 35 bits of the divisor (36 bits less the sign bit)
2817 
2818   // divisor format: . d(0) .... d(34) 0 0 0 .... 0
2819 
2820 
2821 
2822 
2823 
2824 
2825 
2826 
2827 
2828 
2829   // divisor goes in the low half
2830   uint128 dFrac = cpu.CY & MASK35;
2831   if (divisorNegative) {
2832     dFrac = ~dFrac + 1;
2833     sign = - sign;
2834   }
2835   dFrac &= MASK35;
2836 
2837 
2838   if (dFrac == 0) {
2839 //sim_printf ("DVFa A %012"PRIo64" Q %012"PRIo64" Y %012"PRIo64"\r\n", cpu.rA, cpu.rQ, cpu.CY);
2840 // case 1: 400000000000 000000000000 000000000000 --> 400000000000 000000000000
2841 //         dFrac 000000000000 000000000000
2842 
2843     cpu.rA = (zFrac >> 35) & MASK35;
2844     cpu.rQ = (zFrac & MASK35) << 1;
2845 
2846     SC_I_ZERO (cpu.CY == 0);
2847     SC_I_NEG (cpu.rA & SIGN36);
2848 // /* if (current_running_cpu_idx) */
2849 //        sim_printf ("dvf 1 divide check fault A %012llo B %012llo Y %012llo Z %d N %d\r\n",
2850 //                    cpu.rA, cpu.rQ, cpu.CY, TST_I_ZERO, TST_I_NEG);
2851     doFault(FAULT_DIV, fst_zero, "DVF: divide check fault");
2852   }
2853 
2854   //buf [0] = 0;
2855   //print_int128o (dFrac, buf);
2856   //sim_printf ("dFrac %s \r\n", buf);
2857 
2858   L68_ (cpu.ou.cycle |= ou_GD2;)
2859   uint128 sn = zFrac;
2860   word36 quot = 0;
2861   for (uint i = 0; i < 35; i ++) {
2862     // 71 bit number
2863     uint128 s2n = sn << 1;
2864     if (s2n > dFrac) {
2865       s2n -= dFrac;
2866       quot |= (1llu << (34 - i));
2867     }
2868     sn = s2n;
2869     //buf [0] = 0;
2870     //print_int128o (sn, buf);
2871     //sim_printf ("sn %s \r\n", buf);
2872     //buf [0] = 0;
2873     //print_int128o (quot, buf);
2874     //sim_printf ("quot %s \r\n", buf);
2875   }
2876   word36 remainder = sn;
2877 
2878   if (sign == -1)
2879     quot = ~quot + 1;
2880 
2881   if (dividendNegative)
2882     remainder = ~remainder + 1;
2883 
2884   L68_ (cpu.ou.cycle |= ou_GD2;)
2885     // I am surmising that the "If | dividend | >= | divisor |" is an
2886     // overflow prediction; implement it by checking that the calculated
2887     // quotient will fit in 35 bits.
2888 
2889   if (quot & ~((uint128) MASK35)) {
2890     SC_I_ZERO (cpu.rA == 0);
2891     SC_I_NEG (cpu.rA & SIGN36);
2892 // /* if (current_running_cpu_idx) */
2893 //       sim_printf ("dvf 2 divide check fault A %012llo B %012llo Y %012llo Z %d N %d\r\n",
2894 //                   cpu.rA, cpu.rQ, cpu.CY, TST_I_ZERO, TST_I_NEG);
2895     doFault(FAULT_DIV, fst_zero, "DVF: divide check fault");
2896   }
2897   cpu.rA = quot & MASK36;
2898 # if defined(TESTING)
2899   HDBGRegAW ("dvf");
2900 # endif
2901   cpu.rQ = remainder & MASK36;
2902 
2903 #endif
2904 
2905 // MM code
2906 #if defined(DVF_CAC)
2907 
2908 
2909 
2910 
2911 
2912 
2913 
2914 
2915 
2916 
2917 
2918 
2919 
2920 
2921 
2922 // /* if (current_running_cpu_idx) */ sim_printf ("dvf 0 A %012llo B %012llo Y %012llo\r\n", cpu.rA, cpu.rQ, cpu.CY);
2923     // dividend format
2924     // 0  1     70 71
2925     // s  dividend x
2926     //  C(AQ)
2927 
2928     L68_ (cpu.ou.cycle |= ou_GD1;)
2929     int sign = 1;
2930     bool dividendNegative = (getbits36_1 (cpu . rA, 0) != 0);
2931     bool divisorNegative = (getbits36_1 (cpu.CY, 0) != 0);
2932 
2933     // Get the 70 bits of the dividend (72 bits less the sign bit and the
2934     // ignored bit 71.
2935 
2936     // dividend format:   . d(0) ...  d(69)
2937 
2938 # if defined(NEED_128)
2939     uint128 zFrac = lshift_128 (construct_128 (0, cpu.rA & MASK35), 35);
2940     zFrac = or_128 (zFrac, construct_128 (0, (cpu.rQ >> 1) & MASK35));
2941 //sim_debug (DBG_TRACEEXT, & cpu_dev, "zfrac %016"PRIx64" %016"PRIx64"\n", zFrac.h, zFrac.l);
2942 # else
2943     uint128 zFrac = ((uint128) (cpu . rA & MASK35) << 35) | ((cpu . rQ >> 1) & MASK35);
2944 //sim_debug (DBG_TRACEEXT, & cpu_dev, "zfrac %016"PRIx64" %016"PRIx64"\n", (uint64) (zFrac>>64), (uint64) zFrac);
2945 # endif
2946     //zFrac <<= 1; -- Makes Multics unbootable.
2947 
2948     if (dividendNegative)
2949       {
2950 # if defined(NEED_128)
2951         zFrac = negate_128 (zFrac);
2952 # else
2953         // zFrac = ~zFrac + 1;
2954         // ubsan
2955         zFrac = (uint128) (((int128) (~zFrac)) + 1);
2956 # endif
2957         sign = - sign;
2958       }
2959 # if defined(NEED_128)
2960     zFrac = and_128 (zFrac, MASK70);
2961 //sim_debug (DBG_TRACEEXT, & cpu_dev, "zfrac %016"PRIx64" %016"PRIx64"\n", zFrac.h, zFrac.l);
2962 # else
2963     zFrac &= MASK70;
2964 //sim_debug (DBG_TRACEEXT, & cpu_dev, "zfrac %016"PRIx64" %016"PRIx64"\n", (uint64) (zFrac>>64), (uint64) zFrac);
2965 # endif
2966 
2967   //char buf [128];
2968   //buf [0] = 0;
2969   //print_int128o (zFrac, buf);
2970   //sim_printf ("zFrac %s \r\n", buf);
2971     //char buf [128] = "";
2972     //print_int128 (zFrac, buf);
2973     //sim_debug (DBG_CAC, & cpu_dev, "zFrac %s\n", buf);
2974 
2975     // Get the 35 bits of the divisor (36 bits less the sign bit)
2976 
2977     // divisor format: . d(0) .... d(34) 0 0 0 .... 0
2978 
2979     // divisor goes in the low half
2980     uint128 dFrac = convert_to_word72 (0, cpu.CY & MASK35);
2981 # if defined(NEED_128)
2982 //sim_debug (DBG_TRACEEXT, & cpu_dev, "dfrac %016"PRIx64" %016"PRIx64"\n", dFrac.h, dFrac.l);
2983 # else
2984 //sim_debug (DBG_TRACEEXT, & cpu_dev, "dfrac %016"PRIx64" %016"PRIx64"\n", (uint64) (dFrac>>64), (uint64) dFrac);
2985 # endif
2986     if (divisorNegative)
2987       {
2988 # if defined(NEED_128)
2989         dFrac = negate_128 (dFrac);
2990 # else
2991         // dFrac = ~dFrac + 1;
2992         // ubsan
2993         dFrac = (uint128) (((int128) (~dFrac)) + 1);
2994 # endif
2995         sign = - sign;
2996       }
2997 # if defined(NEED_128)
2998     dFrac = and_128 (dFrac, construct_128 (0, MASK35));
2999 //sim_debug (DBG_TRACEEXT, & cpu_dev, "dfrac %016"PRIx64" %016"PRIx64"\n", dFrac.h, dFrac.l);
3000 # else
3001     dFrac &= MASK35;
3002 //sim_debug (DBG_TRACEEXT, & cpu_dev, "dfrac %016"PRIx64" %016"PRIx64"\n", (uint64) (dFrac>>64), (uint64) dFrac);
3003 # endif
3004 
3005     //char buf2 [128] = "";
3006     //print_int128 (dFrac, buf2);
3007     //sim_debug (DBG_CAC, & cpu_dev, "dFrac %s\n", buf2);
3008 
3009     //if (dFrac == 0 || zFrac >= dFrac)
3010     //if (dFrac == 0 || zFrac >= dFrac << 35)
3011 # if defined(NEED_128)
3012     if (iszero_128 (dFrac))
3013 # else
3014     if (dFrac == 0)
3015 # endif
3016       {
3017 // case 1: 400000000000 000000000000 000000000000 --> 400000000000 000000000000
3018 //         dFrac 000000000000 000000000000
3019 
3020         //cpu . rA = (zFrac >> 35) & MASK35;
3021         //cpu . rQ = (word36) ((zFrac & MASK35) << 1);
3022 // ISOLTS 730 expects the right to be zero and the sign
3023 // bit to be untouched.
3024         cpu.rQ = cpu.rQ & (MASK35 << 1);
3025 
3026         //SC_I_ZERO (dFrac == 0);
3027         //SC_I_NEG (cpu . rA & SIGN36);
3028         SC_I_ZERO (cpu.CY == 0);
3029         SC_I_NEG (cpu.rA & SIGN36);
3030         dlyDoFault(FAULT_DIV, fst_zero, "DVF: divide check fault");
3031         return;
3032       }
3033   //buf [0] = 0;
3034   //print_int128o (dFrac, buf);
3035   //sim_printf ("dFrac %s \r\n", buf);
3036 
3037     L68_ (cpu.ou.cycle |= ou_GD2;)
3038 # if defined(NEED_128)
3039     uint128 remainder;
3040     uint128 quot = divide_128 (zFrac, dFrac, & remainder);
3041 //sim_debug (DBG_TRACEEXT, & cpu_dev, "remainder %016"PRIx64" %016"PRIx64"\n", remainder.h, remainder.l);
3042 //sim_debug (DBG_TRACEEXT, & cpu_dev, "quot %016"PRIx64" %016"PRIx64"\n", quot.h, quot.l);
3043 # else
3044     uint128 quot = zFrac / dFrac;
3045     uint128 remainder = zFrac % dFrac;
3046 //sim_debug (DBG_TRACEEXT, & cpu_dev, "remainder %016"PRIx64" %016"PRIx64"\n", (uint64) (remainder>>64), (uint64) remainder);
3047 //sim_debug (DBG_TRACEEXT, & cpu_dev, "quot %016"PRIx64" %016"PRIx64"\n", (uint64) (quot>>64), (uint64) quot);
3048 # endif
3049 
3050     // I am surmising that the "If | dividend | >= | divisor |" is an
3051     // overflow prediction; implement it by checking that the calculated
3052     // quotient will fit in 35 bits.
3053 
3054 # if defined(NEED_128)
3055     if (isnonzero_128 (and_128 (quot, construct_128 (MASK36,  ~MASK35))))
3056 # else
3057     if (quot & ~MASK35)
3058 # endif
3059       {
3060 //
3061 // this got:
3062 //            s/b 373737373737 373737373740 200200
3063 //            was 373737373740 373737373740 000200
3064 //                          ~~
3065 
3066         bool Aneg = (cpu.rA & SIGN36) != 0; // blood type
3067         bool AQzero = cpu.rA == 0 && cpu.rQ == 0;
3068         if (cpu.rA & SIGN36)
3069           {
3070             cpu.rA = (~cpu.rA) & MASK36;
3071             cpu.rQ = (~cpu.rQ) & MASK36;
3072             cpu.rQ += 1;
3073             if (cpu.rQ & BIT37) // overflow?
3074               {
3075                 cpu.rQ &= MASK36;
3076                 cpu.rA = (cpu.rA + 1) & MASK36;
3077               }
3078           }
3079 
3080 
3081 
3082 
3083 
3084 
3085 
3086 # if defined(TESTING)
3087         HDBGRegAW ("dvf");
3088 # endif
3089         //cpu . rA = (zFrac >> 35) & MASK35;
3090         //cpu . rQ = (word36) ((zFrac & MASK35) << 1);
3091 // ISOLTS 730 expects the right to be zero and the sign
3092 // bit to be untouched.
3093         cpu.rQ = cpu.rQ & (MASK35 << 1);
3094 
3095         //SC_I_ZERO (dFrac == 0);
3096         //SC_I_NEG (cpu . rA & SIGN36);
3097         SC_I_ZERO (AQzero);
3098         SC_I_NEG (Aneg);
3099 
3100         dlyDoFault(FAULT_DIV, fst_zero, "DVF: divide check fault");
3101         return;
3102       }
3103     //char buf3 [128] = "";
3104     //print_int128 (remainder, buf3);
3105     //sim_debug (DBG_CAC, & cpu_dev, "remainder %s\n", buf3);
3106 
3107     if (sign == -1)
3108 # if defined(NEED_128)
3109       quot = negate_128 (quot);
3110 # else
3111       quot = ~quot + 1;
3112 # endif
3113 
3114     if (dividendNegative)
3115 # if defined(NEED_128)
3116       remainder = negate_128 (remainder);
3117 # else
3118       remainder = ~remainder + 1;
3119 # endif
3120 
3121 # if defined(NEED_128)
3122     cpu.rA = quot.l & MASK36;
3123     cpu.rQ = remainder.l & MASK36;
3124 # else
3125     cpu . rA = quot & MASK36;
3126     cpu . rQ = remainder & MASK36;
3127 # endif
3128 # if defined(TESTING)
3129     HDBGRegAW ("dvf");
3130 # endif
3131 #endif
3132 
3133 //sim_debug (DBG_CAC, & cpu_dev, "Quotient %"PRId64" (%"PRIo64")\n", cpu . rA, cpu . rA);
3134 //sim_debug (DBG_CAC, & cpu_dev, "Remainder %"PRId64"\n", cpu . rQ);
3135     SC_I_ZERO (cpu . rA == 0 && cpu . rQ == 0);
3136     SC_I_NEG (cpu . rA & SIGN36);
3137 // /* if (current_running_cpu_idx) */
3138 //      sim_printf ("dvf 3 A %012llo B %012llo Y %012llo Z %d N %d\r\n",
3139 //                  cpu.rA, cpu.rQ, cpu.CY, TST_I_ZERO, TST_I_NEG);
3140 }
3141 
3142 /*!
3143  * double precision floating round ...
3144  */
3145 void dfrd (cpu_state_t * cpup) {
     /*  */
3146   // The dfrd instruction is identical to the frd instruction except that the
3147   // rounding constant used is (11...1)65,71 instead of (11...1)29,71.
3148 
3149   // If C(AQ) != 0, the frd instruction performs a true round to a precision
3150   // of 64 bits and a normalization on C(EAQ).
3151   // A true round is a rounding operation such that the sum of the result of
3152   // applying the operation to two numbers of equal magnitude but opposite
3153   // sign is exactly zero.
3154 
3155   // The frd instruction is executed as follows:
3156   // C(AQ) + (11...1)65,71 -> C(AQ)
3157   // * If C(AQ)0 = 0, then a carry is added at AQ71
3158   // * If overflow occurs, C(AQ) is shifted one place to the right and C(E)
3159   // is increased by 1.
3160   // * If overflow does not occur, C(EAQ) is normalized.
3161   // * If C(AQ) = 0, C(E) is set to -128 and the zero indicator is set ON.
3162 
3163   CPTUR (cptUseE);
3164   float72 m = convert_to_word72 (cpu.rA, cpu.rQ);
3165   if (ISZERO_128 (m)) {
3166     cpu.rE = 0200U; /*-128*/
3167     SET_I_ZERO;
3168     CLR_I_NEG;
3169     return;
3170   }
3171 
3172   // C(AQ) + (11...1)65,71 -> C(AQ)
3173   bool ovf;
3174   word18 flags1 = 0;
3175   word1 carry = 0;
3176   // If C(AQ)0 = 0, then a carry is added at AQ71
3177 #if defined(NEED_128)
3178   if (iszero_128 (and_128 (m, SIGN72)))
3179 #else
3180   if ((m & SIGN72) == 0)
3181 #endif
3182   {
3183     carry = 1;
3184   }
3185 #if defined(NEED_128)
3186   m = Add72b (cpup, m, construct_128 (0, 0177), carry, I_OFLOW, & flags1, & ovf);
3187 #else
3188   m = Add72b (cpup, m, 0177, carry, I_OFLOW, & flags1, & ovf);
3189 #endif
3190 
3191   // 0 -> C(AQ)64,71
3192 #if defined(NEED_128)
3193   putbits72 (& m, 64, 8, construct_128 (0, 0));  // 64-71 => 0 per DH02
3194 #else
3195   putbits72 (& m, 64, 8, 0);  // 64-71 => 0 per DH02
3196 #endif
3197 
3198   // If overflow occurs, C(AQ) is shifted one place to the right and C(E) is
3199   // increased by 1.
3200   // If overflow does not occur, C(EAQ) is normalized.
3201   // All of this is done by fno, we just need to save the overflow flag
3202 
3203   bool savedovf = TST_I_OFLOW;
3204   SC_I_OFLOW (ovf);
3205   convert_to_word36 (m, & cpu.rA, & cpu.rQ);
3206 
3207   fno (cpup, & cpu.rE, & cpu.rA, & cpu.rQ);
3208 #if defined(TESTING)
3209   HDBGRegAW ("dfrd");
3210 #endif
3211   SC_I_OFLOW (savedovf);
3212 }
3213 
3214 void dfstr (cpu_state_t * cpup, word36 *Ypair)
     /*  */
3215 {
3216     // The dfstr instruction performs a double-precision true round and normalization on C(EAQ) as it is stored.
3217     // The definition of true round is located under the description of the frd instruction.
3218     // The definition of normalization is located under the description of the fno instruction.
3219     // Except for the precision of the stored result, the dfstr instruction is identical to the fstr instruction.
3220 
3221     // The dfrd instruction is identical to the frd instruction except that the rounding constant used is
3222     // (11...1)65,71 instead of (11...1)29,71.
3223 
3224     // If C(AQ) ≠ 0, the frd instruction performs a true round to a precision of 28 bits and a normalization on C(EAQ).
3225     // A true round is a rounding operation such that the sum of the result of applying the operation to two numbers of
3226     // equal magnitude but opposite sign is exactly zero.
3227 
3228     // The frd instruction is executed as follows:
3229     // C(AQ) + (11...1)29,71 → C(AQ)
3230     // * If C(AQ)0 = 0, then a carry is added at AQ71
3231     // * If overflow occurs, C(AQ) is shifted one place to the right and C(E) is increased by 1.
3232     // * If overflow does not occur, C(EAQ) is normalized.
3233     // * If C(AQ) = 0, C(E) is set to -128 and the zero indicator is set ON.
3234 
3235     // I believe AL39 is incorrect; bits 64-71 should be set to 0, not 65-71. DH02-01 & Bull 9000 is correct.
3236 
3237     CPTUR (cptUseE);
3238     word36 A = cpu . rA, Q = cpu . rQ;
3239     word8 E = cpu . rE;
3240     //A &= DMASK;
3241     //Q &= DMASK;
3242 
3243     float72 m = convert_to_word72 (A, Q);
3244 #if defined(NEED_128)
3245     if (iszero_128 (m))
3246 #else
3247     if (m == 0)
3248 #endif
3249     {
3250         E = (word8)-128;
3251         SET_I_ZERO;
3252         CLR_I_NEG;
3253 
3254         Ypair[0] = ((word36) E & MASK8) << 28;
3255         Ypair[1] = 0;
3256 
3257         return;
3258     }
3259 
3260     // C(AQ) + (11...1)65,71 → C(AQ)
3261     bool ovf;
3262     word18 flags1 = 0;
3263     word1 carry = 0;
3264     // If C(AQ)0 = 0, then a carry is added at AQ71
3265 #if defined(NEED_128)
3266     if (iszero_128 (and_128 (m, SIGN72)))
3267 #else
3268     if ((m & SIGN72) == 0)
3269 #endif
3270       {
3271         carry = 1;
3272       }
3273 #if defined(NEED_128)
3274     m = Add72b (cpup, m, construct_128 (0, 0177), carry, I_OFLOW, & flags1, & ovf);
3275 #else
3276     m = Add72b (cpup, m, 0177, carry, I_OFLOW, & flags1, & ovf);
3277 #endif
3278 
3279     // 0 -> C(AQ)65,71  (per. RJ78)
3280 #if defined(NEED_128)
3281     putbits72 (& m, 64, 8, construct_128 (0, 0));  // 64-71 => 0 per DH02
3282 #else
3283     putbits72 (& m, 64, 8, 0);  // 64-71 => 0 per DH02
3284 
3285 #endif
3286 
3287     // If overflow occurs, C(AQ) is shifted one place to the right and C(E) is
3288     // increased by 1.
3289     // If overflow does not occur, C(EAQ) is normalized.
3290     // All of this is done by fno, we just need to save the overflow flag
3291 
3292     bool savedovf = TST_I_OFLOW;
3293     SC_I_OFLOW(ovf);
3294     convert_to_word36 (m, & A, & Q);
3295 
3296     fno (cpup, & E, & A, & Q);
3297     SC_I_OFLOW(savedovf);
3298 
3299     Ypair[0] = (((word36)E & MASK8) << 28) | ((A & 0777777777400LL) >> 8);
3300     Ypair[1] = ((A & 0377) << 28) | ((Q & 0777777777400LL) >> 8);
3301 }
3302 
3303 /*
3304  * double precision Floating Compare ...
3305  */
3306 void dfcmp (cpu_state_t * cpup) {
     /*  */
3307   // C(E) :: C(Y)0,7
3308   // C(AQ)0,63 :: C(Y-pair)8,71
3309   // * Zero: If | C(EAQ)| = | C(Y-pair) |, then ON; otherwise OFF
3310   // * Neg : If | C(EAQ)| < | C(Y-pair) |, then ON; otherwise OFF
3311 
3312   // The dfcmp instruction is identical to the fcmp instruction except for
3313   // the precision of the mantissas actually compared.
3314 
3315   // Notes: The fcmp instruction is executed as follows:
3316   // The mantissas are aligned by shifting the mantissa of the operand with
3317   // the algebraically smaller exponent to the right the number of places
3318   // equal to the difference in the two exponents.
3319   // The aligned mantissas are compared and the indicators set accordingly.
3320 
3321   // C(AQ)0,63
3322   CPTUR (cptUseE);
3323 #if defined(HEX_MODE)
3324   uint shift_amt = isHex(cpup) ? 4 : 1;
3325 #else
3326   uint shift_amt = 1;
3327 #endif
3328   word72 m1 = convert_to_word72 (cpu.rA, cpu.rQ & 0777777777400LL);
3329   int   e1 = SIGNEXT8_int (cpu . rE & MASK8);
3330 
3331   // C(Y-pair)8,71
3332 #if defined(NEED_128)
3333   word72 m2 = lshift_128 (construct_128 (0, getbits36_28 (cpu.Ypair[0], 8)), (36 + 8));
3334   m2 = or_128 (m2, lshift_128 (construct_128 (0, cpu.Ypair[1]), 8u));
3335 #else
3336   word72 m2 = (word72) getbits36_28 (cpu.Ypair[0], 8) << (36 + 8);
3337   m2 |= cpu.Ypair[1] << 8;
3338 #endif
3339   int   e2 = SIGNEXT8_int (getbits36_8 (cpu.Ypair[0], 0));
3340 
3341   // Which exponent is smaller???
3342 
3343   int shift_count = -1;
3344   word1 notallzeros = 0;
3345 
3346 #if defined(NEED_128)
3347   if (e1 == e2) {
3348     shift_count = 0;
3349   } else if (e1 < e2) {
3350     shift_count = abs (e2 - e1) * (int) shift_amt;
3351     bool s = isnonzero_128 (and_128 (m1, SIGN72));   ///< mantissa negative?
3352     for (int n = 0; n < shift_count; n += 1) {
3353       notallzeros |= m1.l & 1;
3354       m1 = rshift_128 (m1, 1);
3355       if (s)
3356         m1 = or_128 (m1, SIGN72);
3357     }
3358 
3359 # if defined(HEX_MODE)
3360     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3361       m1 = construct_128 (0, 0);
3362 # else
3363     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count > 71)
3364       m1 = construct_128 (0, 0);
3365 # endif
3366     m1 = and_128 (m1, MASK72);
3367   } else { // e2 < e1;
3368     shift_count = abs(e1 - e2) * (int) shift_amt;
3369     bool s = isnonzero_128 (and_128 (m2, SIGN72));   ///< mantissa negative?
3370     for (int n = 0; n < shift_count; n += 1) {
3371       notallzeros |= m2.l & 1;
3372       m2 = rshift_128 (m2, 1);
3373       if (s)
3374         m2 = or_128 (m2, SIGN72);
3375     }
3376 # if defined(HEX_MODE)
3377     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3378       m2 = construct_128 (0, 0);
3379 # else
3380     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count > 71)
3381       m2 = construct_128 (0, 0);
3382 # endif
3383     m2 = and_128 (m2, MASK72);
3384   }
3385 
3386   SC_I_ZERO (iseq_128 (m1, m2));
3387   int128 sm1 = SIGNEXT72_128 (m1);
3388   int128 sm2 = SIGNEXT72_128 (m2);
3389   SC_I_NEG (islt_s128 (sm1, sm2));
3390 #else // NEED_128
3391   if (e1 == e2) {
3392     shift_count = 0;
3393   } else if (e1 < e2) {
3394     shift_count = abs(e2 - e1) * (int) shift_amt;
3395     bool s = m1 & SIGN72;   ///< mantissa negative?
3396     for (int n = 0; n < shift_count; n += 1) {
3397       notallzeros |= m1 & 1;
3398       m1 >>= 1;
3399       if (s)
3400         m1 |= SIGN72;
3401     }
3402 
3403 # if defined(HEX_MODE)
3404     if (m1 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3405       m1 = 0;
3406 # else
3407     if (m1 == MASK72 && notallzeros == 1 && shift_count > 71)
3408       m1 = 0;
3409 # endif
3410     m1 &= MASK72;
3411   } else { // e2 < e1;
3412     shift_count = abs(e1 - e2) * (int) shift_amt;
3413     bool s = m2 & SIGN72;   ///< mantissa negative?
3414     for (int n = 0; n < shift_count; n += 1) {
3415       notallzeros |= m2 & 1;
3416       m2 >>= 1;
3417       if (s)
3418         m2 |= SIGN72;
3419     }
3420 # if defined(HEX_MODE)
3421     if (m2 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3422       m2 = 0;
3423 # else
3424     if (m2 == MASK72 && notallzeros == 1 && shift_count > 71)
3425       m2 = 0;
3426 # endif
3427     m2 &= MASK72;
3428   }
3429 
3430   SC_I_ZERO (m1 == m2);
3431   int128 sm1 = SIGNEXT72_128 (m1);
3432   int128 sm2 = SIGNEXT72_128 (m2);
3433   SC_I_NEG (sm1 < sm2);
3434 #endif // NEED_128
3435 }
3436 
3437 /*
3438  * double precision Floating Compare magnitude ...
3439  */
3440 void dfcmg (cpu_state_t * cpup) {
     /*  */
3441   // C(E) :: C(Y)0,7
3442   // | C(AQ)0,27 | :: | C(Y)8,35 |
3443   // * Zero: If | C(EAQ)| = | C(Y) |, then ON; otherwise OFF
3444   // * Neg : If | C(EAQ)| < | C(Y) |, then ON; otherwise OFF
3445 
3446   // Notes: The fcmp instruction is executed as follows:
3447   // The mantissas are aligned by shifting the mantissa of the operand with
3448   // the algebraically smaller exponent to the right the number of places
3449   // equal to the difference in the two exponents.
3450   // The aligned mantissas are compared and the indicators set accordingly.
3451 
3452   // The dfcmg instruction is identical to the dfcmp instruction except that
3453   // the magnitudes of the mantissas are compared instead of the algebraic
3454   // values.
3455 
3456   CPTUR (cptUseE);
3457 #if defined(HEX_MODE)
3458   uint shift_amt = isHex (cpup) ? 4 : 1;
3459 #else
3460   uint shift_amt = 1;
3461 #endif
3462   // C(AQ)0,63
3463   word72 m1 = convert_to_word72 (cpu.rA & MASK36, cpu.rQ & 0777777777400LL);
3464   int    e1 = SIGNEXT8_int (cpu.rE & MASK8);
3465 
3466   // C(Y-pair)8,71
3467 #if defined(NEED_128)
3468   word72 m2 = lshift_128 (construct_128 (0, getbits36_28 (cpu.Ypair[0], 8)), (36 + 8));
3469   m2 = or_128 (m2, lshift_128 (construct_128 (0, cpu.Ypair[1]), 8u));
3470 #else
3471   word72 m2 = (word72) getbits36_28 (cpu.Ypair[0], 8) << (36 + 8);
3472   m2 |= cpu.Ypair[1] << 8;
3473 #endif
3474   int    e2 = SIGNEXT8_int (getbits36_8 (cpu.Ypair[0], 0));
3475 
3476   // Which exponent is smaller???
3477   L68_ (cpu.ou.cycle = ou_GOE;)
3478   int shift_count = -1;
3479   word1 notallzeros = 0;
3480 
3481   if (e1 == e2) {
3482     shift_count = 0;
3483   } else if (e1 < e2) {
3484     L68_ (cpu.ou.cycle = ou_GOA;)
3485     shift_count = abs(e2 - e1) * (int) shift_amt;
3486 #if defined(NEED_128)
3487     bool s = isnonzero_128 (and_128 (m1, SIGN72));   ///< mantissa negative?
3488     for (int n = 0; n < shift_count; n += 1) {
3489       notallzeros |= m1.l & 1;
3490       m1 = rshift_128 (m1, 1);
3491       if (s)
3492         m1 = or_128 (m1, SIGN72);
3493     }
3494 # if defined(HEX_MODE)
3495     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3496       m1 = construct_128 (0, 0);
3497 # else
3498     if (iseq_128 (m1, MASK72) && notallzeros == 1 && shift_count > 71)
3499       m1 = construct_128 (0, 0);
3500 # endif
3501     m1 = and_128 (m1, MASK72);
3502 #else
3503     bool s = m1 & SIGN72;   ///< mantissa negative?
3504     for (int n = 0; n < shift_count; n += 1) {
3505       notallzeros |= m1 & 1;
3506       m1 >>= 1;
3507       if (s)
3508         m1 |= SIGN72;
3509     }
3510 # if defined(HEX_MODE)
3511     if (m1 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3512       m1 = 0;
3513 # else
3514     if (m1 == MASK72 && notallzeros == 1 && shift_count > 71)
3515       m1 = 0;
3516 # endif
3517     m1 &= MASK72;
3518 #endif
3519   } else { // e2 < e1;
3520     shift_count = abs(e1 - e2) * (int) shift_amt;
3521 #if defined(NEED_128)
3522     bool s = isnonzero_128 (and_128 (m2, SIGN72));   ///< mantissa negative?
3523     for (int n = 0; n < shift_count; n += 1) {
3524       notallzeros |= m2.l & 1;
3525       m2 = rshift_128 (m2, 1);
3526       if (s)
3527         m2 = or_128 (m2, SIGN72);
3528     }
3529 # if defined(HEX_MODE)
3530     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3531       m2 = construct_128 (0, 0);
3532 # else
3533     if (iseq_128 (m2, MASK72) && notallzeros == 1 && shift_count > 71)
3534       m2 = construct_128 (0, 0);
3535 # endif
3536     m2 = and_128 (m2, MASK72);
3537 #else
3538     bool s = m2 & SIGN72;   ///< mantissa negative?
3539     for (int n = 0; n < shift_count; n += 1) {
3540       notallzeros |= m2 & 1;
3541       m2 >>= 1;
3542       if (s)
3543         m2 |= SIGN72;
3544     }
3545 # if defined(HEX_MODE)
3546     if (m2 == MASK72 && notallzeros == 1 && shift_count * (int) shift_amt > 71)
3547       m2 = 0;
3548 # else
3549     if (m2 == MASK72 && notallzeros == 1 && shift_count > 71)
3550       m2 = 0;
3551 # endif
3552     m2 &= MASK72;
3553 #endif
3554   }
3555 
3556 #if defined(NEED_128)
3557   SC_I_ZERO (iseq_128 (m1, m2));
3558   int128 sm1 = SIGNEXT72_128 (m1);
3559   if (sm1.h < 0)
3560     sm1 = negate_s128 (sm1);
3561   int128 sm2 = SIGNEXT72_128 (m2);
3562   if (sm2.h < 0)
3563     sm2 = negate_s128 (sm2);
3564 
3565   SC_I_NEG (islt_s128 (sm1, sm2));
3566 #else // ! NEED_128
3567   SC_I_ZERO (m1 == m2);
3568   int128 sm1 = SIGNEXT72_128 (m1);
3569   if (sm1 < 0)
3570     sm1 = - sm1;
3571   int128 sm2 = SIGNEXT72_128 (m2);
3572   if (sm2 < 0)
3573     sm2 = - sm2;
3574 
3575   SC_I_NEG (sm1 < sm2);
3576 #endif // ! NEED_128
3577 }
3578 
3579 #if defined(__NVCOMPILER) || defined(__NVCOMPILER_LLVM__) || defined(__PGI) || defined(__PGLLVM__)
3580 # pragma global vector
3581 #endif
/* */
root/src/dps8/dps8_math.c

DEFINITIONS
