libphobos/libdruntime/core/int128.d

   1 /* 128 bit integer arithmetic.
   2  *
   3  * Not optimized for speed.
   4  *
   5  * Copyright: Copyright D Language Foundation 2022.
   6  * License:   $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0)
   7  * Authors:   Walter Bright
   8  * Source:    $(DRUNTIMESRC core/_int128.d)
   9  */
  10
  11 module core.int128;
  12
  13 nothrow:
  14 @safe:
  15 @nogc:
  16
  17 private alias I = long;
  18 private alias U = ulong;
  19 enum Ubits = uint(U.sizeof * 8);
  20
  21 version (DigitalMars)
  22 {
  23     /* The alignment should follow target.stackAlign(),
  24      * which is `isXmmSupported() ? 16 : (is64bit ? 8 : 4)
  25      */
  26     version (D_SIMD)
  27         private enum Cent_alignment = 16;
  28     else version (X86_64)
  29         private enum Cent_alignment = 8;
  30     else
  31         private enum Cent_alignment = 4;
  32 }
  33 else
  34 {
  35     version (X86_64) private enum Cent_alignment = 16;
  36     else             private enum Cent_alignment = (size_t.sizeof * 2);
  37 }
  38
  39 /**
  40  * 128 bit integer type.
  41  * See_also: $(REF Int128, std,int128).
  42  */
  43 align(Cent_alignment) struct Cent
  44 {
  45     version (LittleEndian)
  46     {
  47         U lo;  // low 64 bits
  48         U hi;  // high 64 bits
  49     }
  50     else
  51     {
  52         U hi;  // high 64 bits
  53         U lo;  // low 64 bits
  54     }
  55 }
  56
  57 enum Cent One = { lo:1 };
  58 enum Cent Zero = { lo:0 };
  59 enum Cent MinusOne = neg(One);
  60
  61 /*****************************
  62  * Test against 0
  63  * Params:
  64  *      c = Cent to test
  65  * Returns:
  66  *      true if != 0
  67  */
  68 pure
  69 bool tst(Cent c)
  70 {
  71     return c.hi || c.lo;
  72 }
  73
  74
  75 /*****************************
  76  * Complement
  77  * Params:
  78  *      c = Cent to complement
  79  * Returns:
  80  *      complemented value
  81  */
  82 pure
  83 Cent com(Cent c)
  84 {
  85     c.lo = ~c.lo;
  86     c.hi = ~c.hi;
  87     return c;
  88 }
  89
  90 /*****************************
  91  * Negate
  92  * Params:
  93  *      c = Cent to negate
  94  * Returns:
  95  *      negated value
  96  */
  97 pure
  98 Cent neg(Cent c)
  99 {
 100     if (c.lo == 0)
 101         c.hi = -c.hi;
 102     else
 103     {
 104         c.lo = -c.lo;
 105         c.hi = ~c.hi;
 106     }
 107     return c;
 108 }
 109
 110 /*****************************
 111  * Increment
 112  * Params:
 113  *      c = Cent to increment
 114  * Returns:
 115  *      incremented value
 116  */
 117 pure
 118 Cent inc(Cent c)
 119 {
 120     return add(c, One);
 121 }
 122
 123 /*****************************
 124  * Decrement
 125  * Params:
 126  *      c = Cent to decrement
 127  * Returns:
 128  *      incremented value
 129  */
 130 pure
 131 Cent dec(Cent c)
 132 {
 133     return sub(c, One);
 134 }
 135
 136 /*****************************
 137  * Shift left one bit
 138  * Params:
 139  *      c = Cent to shift
 140  * Returns:
 141  *      shifted value
 142  */
 143 pure
 144 Cent shl1(Cent c)
 145 {
 146     c.hi = (c.hi << 1) | (cast(I)c.lo < 0);
 147     c.lo <<= 1;
 148     return c;
 149 }
 150
 151 /*****************************
 152  * Unsigned shift right one bit
 153  * Params:
 154  *      c = Cent to shift
 155  * Returns:
 156  *      shifted value
 157  */
 158 pure
 159 Cent shr1(Cent c)
 160 {
 161     c.lo = (c.lo >> 1) | ((c.hi & 1) << (Ubits - 1));
 162     c.hi >>= 1;
 163     return c;
 164 }
 165
 166
 167 /*****************************
 168  * Arithmetic shift right one bit
 169  * Params:
 170  *      c = Cent to shift
 171  * Returns:
 172  *      shifted value
 173  */
 174 pure
 175 Cent sar1(Cent c)
 176 {
 177     c.lo = (c.lo >> 1) | ((c.hi & 1) << (Ubits - 1));
 178     c.hi = cast(I)c.hi >> 1;
 179     return c;
 180 }
 181
 182 /*****************************
 183  * Shift left n bits
 184  * Params:
 185  *      c = Cent to shift
 186  *      n = number of bits to shift
 187  * Returns:
 188  *      shifted value
 189  */
 190 pure
 191 Cent shl(Cent c, uint n)
 192 {
 193     if (n >= Ubits * 2)
 194         return Zero;
 195
 196     if (n >= Ubits)
 197     {
 198         c.hi = c.lo << (n - Ubits);
 199         c.lo = 0;
 200     }
 201     else
 202     {
 203         c.hi = ((c.hi << n) | (c.lo >> (Ubits - n - 1) >> 1));
 204         c.lo = c.lo << n;
 205     }
 206     return c;
 207 }
 208
 209 /*****************************
 210  * Unsigned shift right n bits
 211  * Params:
 212  *      c = Cent to shift
 213  *      n = number of bits to shift
 214  * Returns:
 215  *      shifted value
 216  */
 217 pure
 218 Cent shr(Cent c, uint n)
 219 {
 220     if (n >= Ubits * 2)
 221         return Zero;
 222
 223     if (n >= Ubits)
 224     {
 225         c.lo = c.hi >> (n - Ubits);
 226         c.hi = 0;
 227     }
 228     else
 229     {
 230         c.lo = ((c.lo >> n) | (c.hi << (Ubits - n - 1) << 1));
 231         c.hi = c.hi >> n;
 232     }
 233     return c;
 234 }
 235
 236 /*****************************
 237  * Arithmetic shift right n bits
 238  * Params:
 239  *      c = Cent to shift
 240  *      n = number of bits to shift
 241  * Returns:
 242  *      shifted value
 243  */
 244 pure
 245 Cent sar(Cent c, uint n)
 246 {
 247     const signmask = -(c.hi >> (Ubits - 1));
 248     const signshift = (Ubits * 2) - n;
 249     c = shr(c, n);
 250
 251     // Sign extend all bits beyond the precision of Cent.
 252     if (n >= Ubits * 2)
 253     {
 254         c.hi = signmask;
 255         c.lo = signmask;
 256     }
 257     else if (signshift >= Ubits * 2)
 258     {
 259     }
 260     else if (signshift >= Ubits)
 261     {
 262         c.hi &= ~(U.max << (signshift - Ubits));
 263         c.hi |= signmask << (signshift - Ubits);
 264     }
 265     else
 266     {
 267         c.hi = signmask;
 268         c.lo &= ~(U.max << signshift);
 269         c.lo |= signmask << signshift;
 270     }
 271     return c;
 272 }
 273
 274 /*****************************
 275  * Rotate left one bit
 276  * Params:
 277  *      c = Cent to rotate
 278  * Returns:
 279  *      rotated value
 280  */
 281 pure
 282 Cent rol1(Cent c)
 283 {
 284     int carry = cast(I)c.hi < 0;
 285
 286     c.hi = (c.hi << 1) | (cast(I)c.lo < 0);
 287     c.lo = (c.lo << 1) | carry;
 288     return c;
 289 }
 290
 291 /*****************************
 292  * Rotate right one bit
 293  * Params:
 294  *      c = Cent to rotate
 295  * Returns:
 296  *      rotated value
 297  */
 298 pure
 299 Cent ror1(Cent c)
 300 {
 301     int carry = c.lo & 1;
 302     c.lo = (c.lo >> 1) | (cast(U)(c.hi & 1) << (Ubits - 1));
 303     c.hi = (c.hi >> 1) | (cast(U)carry << (Ubits - 1));
 304     return c;
 305 }
 306
 307
 308 /*****************************
 309  * Rotate left n bits
 310  * Params:
 311  *      c = Cent to rotate
 312  *      n = number of bits to rotate
 313  * Returns:
 314  *      rotated value
 315  */
 316 pure
 317 Cent rol(Cent c, uint n)
 318 {
 319     n &= Ubits * 2 - 1;
 320     Cent l = shl(c, n);
 321     Cent r = shr(c, Ubits * 2 - n);
 322     return or(l, r);
 323 }
 324
 325 /*****************************
 326  * Rotate right n bits
 327  * Params:
 328  *      c = Cent to rotate
 329  *      n = number of bits to rotate
 330  * Returns:
 331  *      rotated value
 332  */
 333 pure
 334 Cent ror(Cent c, uint n)
 335 {
 336     n &= Ubits * 2 - 1;
 337     Cent r = shr(c, n);
 338     Cent l = shl(c, Ubits * 2 - n);
 339     return or(r, l);
 340 }
 341
 342 /****************************
 343  * And c1 & c2.
 344  * Params:
 345  *      c1 = operand 1
 346  *      c2 = operand 2
 347  * Returns:
 348  *      c1 & c2
 349  */
 350 pure
 351 Cent and(Cent c1, Cent c2)
 352 {
 353     const Cent ret = { lo:c1.lo & c2.lo, hi:c1.hi & c2.hi };
 354     return ret;
 355 }
 356
 357 /****************************
 358  * Or c1 | c2.
 359  * Params:
 360  *      c1 = operand 1
 361  *      c2 = operand 2
 362  * Returns:
 363  *      c1 | c2
 364  */
 365 pure
 366 Cent or(Cent c1, Cent c2)
 367 {
 368     const Cent ret = { lo:c1.lo | c2.lo, hi:c1.hi | c2.hi };
 369     return ret;
 370 }
 371
 372 /****************************
 373  * Xor c1 ^ c2.
 374  * Params:
 375  *      c1 = operand 1
 376  *      c2 = operand 2
 377  * Returns:
 378  *      c1 ^ c2
 379  */
 380 pure
 381 Cent xor(Cent c1, Cent c2)
 382 {
 383     const Cent ret = { lo:c1.lo ^ c2.lo, hi:c1.hi ^ c2.hi };
 384     return ret;
 385 }
 386
 387 /****************************
 388  * Add c1 to c2.
 389  * Params:
 390  *      c1 = operand 1
 391  *      c2 = operand 2
 392  * Returns:
 393  *      c1 + c2
 394  */
 395 pure
 396 Cent add(Cent c1, Cent c2)
 397 {
 398     U r = cast(U)(c1.lo + c2.lo);
 399     const Cent ret = { lo:r, hi:cast(U)(c1.hi + c2.hi + (r < c1.lo)) };
 400     return ret;
 401 }
 402
 403 /****************************
 404  * Subtract c2 from c1.
 405  * Params:
 406  *      c1 = operand 1
 407  *      c2 = operand 2
 408  * Returns:
 409  *      c1 - c2
 410  */
 411 pure
 412 Cent sub(Cent c1, Cent c2)
 413 {
 414     return add(c1, neg(c2));
 415 }
 416
 417 /****************************
 418  * Multiply c1 * c2.
 419  * Params:
 420  *      c1 = operand 1
 421  *      c2 = operand 2
 422  * Returns:
 423  *      c1 * c2
 424  */
 425 pure
 426 Cent mul(Cent c1, Cent c2)
 427 {
 428     enum mulmask = (1UL << (Ubits / 2)) - 1;
 429     enum mulshift = Ubits / 2;
 430
 431     // This algorithm splits the operands into 4 words, then computes and sums
 432     // the partial products of each part.
 433     const c2l0 = c2.lo & mulmask;
 434     const c2l1 = c2.lo >> mulshift;
 435     const c2h0 = c2.hi & mulmask;
 436     const c2h1 = c2.hi >> mulshift;
 437
 438     const c1l0 = c1.lo & mulmask;
 439     U r0 = c1l0 * c2l0;
 440     U r1 = c1l0 * c2l1 + (r0 >> mulshift);
 441     U r2 = c1l0 * c2h0 + (r1 >> mulshift);
 442     U r3 = c1l0 * c2h1 + (r2 >> mulshift);
 443
 444     const c1l1 = c1.lo >> mulshift;
 445     r1 = c1l1 * c2l0 + (r1 & mulmask);
 446     r2 = c1l1 * c2l1 + (r2 & mulmask) + (r1 >> mulshift);
 447     r3 = c1l1 * c2h0 + (r3 & mulmask) + (r2 >> mulshift);
 448
 449     const c1h0 = c1.hi & mulmask;
 450     r2 = c1h0 * c2l0 + (r2 & mulmask);
 451     r3 = c1h0 * c2l1 + (r3 & mulmask) + (r2 >> mulshift);
 452
 453     const c1h1 = c1.hi >> mulshift;
 454     r3 = c1h1 * c2l0 + (r3 & mulmask);
 455
 456     const Cent ret = { lo:(r0 & mulmask) + (r1 & mulmask) * (mulmask + 1),
 457                        hi:(r2 & mulmask) + (r3 & mulmask) * (mulmask + 1) };
 458     return ret;
 459 }
 460
 461
 462 /****************************
 463  * Unsigned divide c1 / c2.
 464  * Params:
 465  *      c1 = dividend
 466  *      c2 = divisor
 467  * Returns:
 468  *      quotient c1 / c2
 469  */
 470 pure
 471 Cent udiv(Cent c1, Cent c2)
 472 {
 473     Cent modulus;
 474     return udivmod(c1, c2, modulus);
 475 }
 476
 477 /****************************
 478  * Unsigned divide c1 / c2. The remainder after division is stored to modulus.
 479  * Params:
 480  *      c1 = dividend
 481  *      c2 = divisor
 482  *      modulus = set to c1 % c2
 483  * Returns:
 484  *      quotient c1 / c2
 485  */
 486 pure
 487 Cent udivmod(Cent c1, Cent c2, out Cent modulus)
 488 {
 489     //printf("udiv c1(%llx,%llx) c2(%llx,%llx)\n", c1.lo, c1.hi, c2.lo, c2.hi);
 490     // Based on "Unsigned Doubleword Division" in Hacker's Delight
 491     import core.bitop;
 492
 493     // Divides a 128-bit dividend by a 64-bit divisor.
 494     // The result must fit in 64 bits.
 495     static U udivmod128_64(Cent c1, U c2, out U modulus)
 496     {
 497         // We work in base 2^^32
 498         enum base = 1UL << 32;
 499         enum divmask = (1UL << (Ubits / 2)) - 1;
 500         enum divshift = Ubits / 2;
 501
 502         // Check for overflow and divide by 0
 503         if (c1.hi >= c2)
 504         {
 505             modulus = 0UL;
 506             return ~0UL;
 507         }
 508
 509         // Computes [num1 num0] / den
 510         static uint udiv96_64(U num1, uint num0, U den)
 511         {
 512             // Extract both digits of the denominator
 513             const den1 = cast(uint)(den >> divshift);
 514             const den0 = cast(uint)(den & divmask);
 515             // Estimate ret as num1 / den1, and then correct it
 516             U ret = num1 / den1;
 517             const t2 = (num1 % den1) * base + num0;
 518             const t1 = ret * den0;
 519             if (t1 > t2)
 520                 ret -= (t1 - t2 > den) ? 2 : 1;
 521             return cast(uint)ret;
 522         }
 523
 524         // Determine the normalization factor. We multiply c2 by this, so that its leading
 525         // digit is at least half base. In binary this means just shifting left by the number
 526         // of leading zeros, so that there's a 1 in the MSB.
 527         // We also shift number by the same amount. This cannot overflow because c1.hi < c2.
 528         const shift = (Ubits - 1) - bsr(c2);
 529         c2 <<= shift;
 530         U num2 = c1.hi;
 531         num2 <<= shift;
 532         num2 |= (c1.lo >> (-shift & 63)) & (-cast(I)shift >> 63);
 533         c1.lo <<= shift;
 534
 535         // Extract the low digits of the numerator (after normalizing)
 536         const num1 = cast(uint)(c1.lo >> divshift);
 537         const num0 = cast(uint)(c1.lo & divmask);
 538
 539         // Compute q1 = [num2 num1] / c2
 540         const q1 = udiv96_64(num2, num1, c2);
 541         // Compute the true (partial) remainder
 542         const rem = num2 * base + num1 - q1 * c2;
 543         // Compute q0 = [rem num0] / c2
 544         const q0 = udiv96_64(rem, num0, c2);
 545
 546         modulus = (rem * base + num0 - q0 * c2) >> shift;
 547         return (cast(U)q1 << divshift) | q0;
 548     }
 549
 550     // Special cases
 551     if (!tst(c2))
 552     {
 553         // Divide by zero
 554         modulus = Zero;
 555         return com(modulus);
 556     }
 557     if (c1.hi == 0 && c2.hi == 0)
 558     {
 559         // Single precision divide
 560         const Cent rem = { lo:c1.lo % c2.lo };
 561         modulus = rem;
 562         const Cent ret = { lo:c1.lo / c2.lo };
 563         return ret;
 564     }
 565     if (c1.hi == 0)
 566     {
 567         // Numerator is smaller than the divisor
 568         modulus = c1;
 569         return Zero;
 570     }
 571     if (c2.hi == 0)
 572     {
 573         // Divisor is a 64-bit value, so we just need one 128/64 division.
 574         // If c1 / c2 would overflow, break c1 up into two halves.
 575         const q1 = (c1.hi < c2.lo) ? 0 : (c1.hi / c2.lo);
 576         if (q1)
 577             c1.hi = c1.hi % c2.lo;
 578         Cent rem;
 579         const q0 = udivmod128_64(c1, c2.lo, rem.lo);
 580         modulus = rem;
 581         const Cent ret = { lo:q0, hi:q1 };
 582         return ret;
 583     }
 584
 585     // Full cent precision division.
 586     // Here c2 >= 2^^64
 587     // We know that c2.hi != 0, so count leading zeros is OK
 588     // We have 0 <= shift <= 63
 589     const shift = (Ubits - 1) - bsr(c2.hi);
 590
 591     // Normalize the divisor so its MSB is 1
 592     // v1 = (c2 << shift) >> 64
 593     U v1 = shl(c2, shift).hi;
 594
 595     // To ensure no overflow.
 596     Cent u1 = shr1(c1);
 597
 598     // Get quotient from divide unsigned operation.
 599     U rem_ignored;
 600     const Cent q1 = { lo:udivmod128_64(u1, v1, rem_ignored) };
 601
 602     // Undo normalization and division of c1 by 2.
 603     Cent quotient = shr(shl(q1, shift), 63);
 604
 605     // Make quotient correct or too small by 1
 606     if (tst(quotient))
 607         quotient = dec(quotient);
 608
 609     // Now quotient is correct.
 610     // Compute rem = c1 - (quotient * c2);
 611     Cent rem = sub(c1, mul(quotient, c2));
 612
 613     // Check if remainder is larger than the divisor
 614     if (uge(rem, c2))
 615     {
 616         // Increment quotient
 617         quotient = inc(quotient);
 618         // Subtract c2 from remainder
 619         rem = sub(rem, c2);
 620     }
 621     modulus = rem;
 622     //printf("quotient "); print(quotient);
 623     //printf("modulus  "); print(modulus);
 624     return quotient;
 625 }
 626
 627
 628 /****************************
 629  * Signed divide c1 / c2.
 630  * Params:
 631  *      c1 = dividend
 632  *      c2 = divisor
 633  * Returns:
 634  *      quotient c1 / c2
 635  */
 636 pure
 637 Cent div(Cent c1, Cent c2)
 638 {
 639     Cent modulus;
 640     return divmod(c1, c2, modulus);
 641 }
 642
 643 /****************************
 644  * Signed divide c1 / c2. The remainder after division is stored to modulus.
 645  * Params:
 646  *      c1 = dividend
 647  *      c2 = divisor
 648  *      modulus = set to c1 % c2
 649  * Returns:
 650  *      quotient c1 / c2
 651  */
 652 pure
 653 Cent divmod(Cent c1, Cent c2, out Cent modulus)
 654 {
 655     /* Muck about with the signs so we can use the unsigned divide
 656      */
 657     if (cast(I)c1.hi < 0)
 658     {
 659         if (cast(I)c2.hi < 0)
 660         {
 661             Cent r = udivmod(neg(c1), neg(c2), modulus);
 662             modulus = neg(modulus);
 663             return r;
 664         }
 665         Cent r = neg(udivmod(neg(c1), c2, modulus));
 666         modulus = neg(modulus);
 667         return r;
 668     }
 669     else if (cast(I)c2.hi < 0)
 670     {
 671         return neg(udivmod(c1, neg(c2), modulus));
 672     }
 673     else
 674         return udivmod(c1, c2, modulus);
 675 }
 676
 677 /****************************
 678  * If c1 > c2 unsigned
 679  * Params:
 680  *      c1 = operand 1
 681  *      c2 = operand 2
 682  * Returns:
 683  *      true if c1 > c2
 684  */
 685 pure
 686 bool ugt(Cent c1, Cent c2)
 687 {
 688     return (c1.hi == c2.hi) ? (c1.lo > c2.lo) : (c1.hi > c2.hi);
 689 }
 690
 691 /****************************
 692  * If c1 >= c2 unsigned
 693  * Params:
 694  *      c1 = operand 1
 695  *      c2 = operand 2
 696  * Returns:
 697  *      true if c1 >= c2
 698  */
 699 pure
 700 bool uge(Cent c1, Cent c2)
 701 {
 702     return !ugt(c2, c1);
 703 }
 704
 705 /****************************
 706  * If c1 < c2 unsigned
 707  * Params:
 708  *      c1 = operand 1
 709  *      c2 = operand 2
 710  * Returns:
 711  *      true if c1 < c2
 712  */
 713 pure
 714 bool ult(Cent c1, Cent c2)
 715 {
 716     return ugt(c2, c1);
 717 }
 718
 719 /****************************
 720  * If c1 <= c2 unsigned
 721  * Params:
 722  *      c1 = operand 1
 723  *      c2 = operand 2
 724  * Returns:
 725  *      true if c1 <= c2
 726  */
 727 pure
 728 bool ule(Cent c1, Cent c2)
 729 {
 730     return !ugt(c1, c2);
 731 }
 732
 733 /****************************
 734  * If c1 > c2 signed
 735  * Params:
 736  *      c1 = operand 1
 737  *      c2 = operand 2
 738  * Returns:
 739  *      true if c1 > c2
 740  */
 741 pure
 742 bool gt(Cent c1, Cent c2)
 743 {
 744     return (c1.hi == c2.hi)
 745         ? (c1.lo > c2.lo)
 746         : (cast(I)c1.hi > cast(I)c2.hi);
 747 }
 748
 749 /****************************
 750  * If c1 >= c2 signed
 751  * Params:
 752  *      c1 = operand 1
 753  *      c2 = operand 2
 754  * Returns:
 755  *      true if c1 >= c2
 756  */
 757 pure
 758 bool ge(Cent c1, Cent c2)
 759 {
 760     return !gt(c2, c1);
 761 }
 762
 763 /****************************
 764  * If c1 < c2 signed
 765  * Params:
 766  *      c1 = operand 1
 767  *      c2 = operand 2
 768  * Returns:
 769  *      true if c1 < c2
 770  */
 771 pure
 772 bool lt(Cent c1, Cent c2)
 773 {
 774     return gt(c2, c1);
 775 }
 776
 777 /****************************
 778  * If c1 <= c2 signed
 779  * Params:
 780  *      c1 = operand 1
 781  *      c2 = operand 2
 782  * Returns:
 783  *      true if c1 <= c2
 784  */
 785 pure
 786 bool le(Cent c1, Cent c2)
 787 {
 788     return !gt(c1, c2);
 789 }
 790
 791 /*******************************************************/
 792
 793 version (unittest)
 794 {
 795     version (none)
 796     {
 797         import core.stdc.stdio : printf;
 798
 799         @trusted
 800         void print(Cent c)
 801         {
 802             printf("%lld, %lld\n", cast(ulong)c.lo, cast(ulong)c.hi);
 803             printf("x%llx, x%llx\n", cast(ulong)c.lo, cast(ulong)c.hi);
 804         }
 805     }
 806 }
 807
 808 unittest
 809 {
 810     const Cent C0 = Zero;
 811     const Cent C1 = One;
 812     const Cent C2 = { lo:2 };
 813     const Cent C3 = { lo:3 };
 814     const Cent C5 = { lo:5 };
 815     const Cent C10 = { lo:10 };
 816     const Cent C20 = { lo:20 };
 817     const Cent C30 = { lo:30 };
 818     const Cent C100 = { lo:100 };
 819
 820     const Cent Cm1 =  neg(One);
 821     const Cent Cm3 =  neg(C3);
 822     const Cent Cm10 = neg(C10);
 823
 824     const Cent C3_1 = { lo:1, hi:3 };
 825     const Cent C3_2 = { lo:2, hi:3 };
 826     const Cent C4_8  = { lo:8, hi:4 };
 827     const Cent C5_0  = { lo:0, hi:5 };
 828     const Cent C7_1 = { lo:1, hi:7 };
 829     const Cent C7_9 = { lo:9, hi:7 };
 830     const Cent C9_3 = { lo:3, hi:9 };
 831     const Cent C10_0 = { lo:0, hi:10 };
 832     const Cent C10_1 = { lo:1, hi:10 };
 833     const Cent C10_3 = { lo:3, hi:10 };
 834     const Cent C11_3 = { lo:3, hi:11 };
 835     const Cent C20_0 = { lo:0, hi:20 };
 836     const Cent C90_30 = { lo:30, hi:90 };
 837
 838     const Cent Cm10_0 = inc(com(C10_0)); // Cent(lo=0,  hi=-10);
 839     const Cent Cm10_1 = inc(com(C10_1)); // Cent(lo=-1, hi=-11);
 840     const Cent Cm10_3 = inc(com(C10_3)); // Cent(lo=-3, hi=-11);
 841     const Cent Cm20_0 = inc(com(C20_0)); // Cent(lo=0,  hi=-20);
 842
 843     enum Cent Cs_3 = { lo:3, hi:I.min };
 844
 845     const Cent Cbig_1 = { lo:0xa3ccac1832952398, hi:0xc3ac542864f652f8 };
 846     const Cent Cbig_2 = { lo:0x5267b85f8a42fc20, hi:0 };
 847     const Cent Cbig_3 = { lo:0xf0000000ffffffff, hi:0 };
 848
 849     /************************/
 850
 851     assert( ugt(C1, C0) );
 852     assert( ult(C1, C2) );
 853     assert( uge(C1, C0) );
 854     assert( ule(C1, C2) );
 855
 856     assert( !ugt(C0, C1) );
 857     assert( !ult(C2, C1) );
 858     assert( !uge(C0, C1) );
 859     assert( !ule(C2, C1) );
 860
 861     assert( !ugt(C1, C1) );
 862     assert( !ult(C1, C1) );
 863     assert( uge(C1, C1) );
 864     assert( ule(C2, C2) );
 865
 866     assert( ugt(C10_3, C10_1) );
 867     assert( ugt(C11_3, C10_3) );
 868     assert( !ugt(C9_3, C10_3) );
 869     assert( !ugt(C9_3, C9_3) );
 870
 871     assert( gt(C2, C1) );
 872     assert( !gt(C1, C2) );
 873     assert( !gt(C1, C1) );
 874     assert( gt(C0, Cm1) );
 875     assert( gt(Cm1, neg(C10)));
 876     assert( !gt(Cm1, Cm1) );
 877     assert( !gt(Cm1, C0) );
 878
 879     assert( !lt(C2, C1) );
 880     assert( !le(C2, C1) );
 881     assert( ge(C2, C1) );
 882
 883     assert(neg(C10_0) == Cm10_0);
 884     assert(neg(C10_1) == Cm10_1);
 885     assert(neg(C10_3) == Cm10_3);
 886
 887     assert(add(C7_1,C3_2) == C10_3);
 888     assert(sub(C1,C2) == Cm1);
 889
 890     assert(inc(C3_1) == C3_2);
 891     assert(dec(C3_2) == C3_1);
 892
 893     assert(shl(C10,0) == C10);
 894     assert(shl(C10,Ubits) == C10_0);
 895     assert(shl(C10,1) == C20);
 896     assert(shl(C10,Ubits * 2) == C0);
 897     assert(shr(C10_0,0) == C10_0);
 898     assert(shr(C10_0,Ubits) == C10);
 899     assert(shr(C10_0,Ubits - 1) == C20);
 900     assert(shr(C10_0,Ubits + 1) == C5);
 901     assert(shr(C10_0,Ubits * 2) == C0);
 902     assert(sar(C10_0,0) == C10_0);
 903     assert(sar(C10_0,Ubits) == C10);
 904     assert(sar(C10_0,Ubits - 1) == C20);
 905     assert(sar(C10_0,Ubits + 1) == C5);
 906     assert(sar(C10_0,Ubits * 2) == C0);
 907     assert(sar(Cm1,Ubits * 2) == Cm1);
 908
 909     assert(shl1(C10) == C20);
 910     assert(shr1(C10_0) == C5_0);
 911     assert(sar1(C10_0) == C5_0);
 912     assert(sar1(Cm1) == Cm1);
 913
 914     Cent modulus;
 915
 916     assert(udiv(C10,C2) == C5);
 917     assert(udivmod(C10,C2, modulus) ==  C5);   assert(modulus == C0);
 918     assert(udivmod(C10,C3, modulus) ==  C3);   assert(modulus == C1);
 919     assert(udivmod(C10,C0, modulus) == Cm1);   assert(modulus == C0);
 920     assert(udivmod(C2,C90_30, modulus) == C0); assert(modulus == C2);
 921     assert(udiv(mul(C90_30, C2), C2) == C90_30);
 922     assert(udiv(mul(C90_30, C2), C90_30) == C2);
 923
 924     assert(div(C10,C3) == C3);
 925     assert(divmod( C10,  C3, modulus) ==  C3); assert(modulus ==  C1);
 926     assert(divmod(Cm10,  C3, modulus) == Cm3); assert(modulus == Cm1);
 927     assert(divmod( C10, Cm3, modulus) == Cm3); assert(modulus ==  C1);
 928     assert(divmod(Cm10, Cm3, modulus) ==  C3); assert(modulus == Cm1);
 929     assert(divmod(C2, C90_30, modulus) == C0); assert(modulus == C2);
 930     assert(div(mul(C90_30, C2), C2) == C90_30);
 931     assert(div(mul(C90_30, C2), C90_30) == C2);
 932
 933     const Cent Cb1divb2 = { lo:0x4496aa309d4d4a2f, hi:U.max };
 934     const Cent Cb1modb2 = { lo:0xd83203d0fdc799b8, hi:U.max };
 935     assert(divmod(Cbig_1, Cbig_2, modulus) == Cb1divb2);
 936     assert(modulus == Cb1modb2);
 937
 938     const Cent Cb1udivb2 = { lo:0x5fe0e9bace2bedad, hi:2 };
 939     const Cent Cb1umodb2 = { lo:0x2c923125a68721f8, hi:0 };
 940     assert(udivmod(Cbig_1, Cbig_2, modulus) == Cb1udivb2);
 941     assert(modulus == Cb1umodb2);
 942
 943     const Cent Cb1divb3 = { lo:0xbfa6c02b5aff8b86, hi:U.max };
 944     const Cent Cb1udivb3 = { lo:0xd0b7d13b48cb350f, hi:0 };
 945     assert(div(Cbig_1, Cbig_3) == Cb1divb3);
 946     assert(udiv(Cbig_1, Cbig_3) == Cb1udivb3);
 947
 948     assert(mul(Cm10, C1) == Cm10);
 949     assert(mul(C1, Cm10) == Cm10);
 950     assert(mul(C9_3, C10) == C90_30);
 951     assert(mul(Cs_3, C10) == C30);
 952     assert(mul(Cm10, Cm10) == C100);
 953     assert(mul(C20_0, Cm1) == Cm20_0);
 954
 955     assert( or(C4_8, C3_1) == C7_9);
 956     assert(and(C4_8, C7_9) == C4_8);
 957     assert(xor(C4_8, C7_9) == C3_1);
 958
 959     assert(rol(Cm1,  1) == Cm1);
 960     assert(ror(Cm1, 45) == Cm1);
 961     assert(rol(ror(C7_9, 5), 5) == C7_9);
 962     assert(rol(C7_9, 1) == rol1(C7_9));
 963     assert(ror(C7_9, 1) == ror1(C7_9));
 964 }