root/src/dps8/dps8_math.c

/* */

DEFINITIONS

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