git.cworth.org Git - vogl/blob - src/voglcore/vogl_bigint128.h

   1 /**************************************************************************
   2  *
   3  * Copyright 2013-2014 RAD Game Tools and Valve Software
   4  * Copyright 2010-2014 Rich Geldreich and Tenacious Software LLC
   5  * All Rights Reserved.
   6  *
   7  * Permission is hereby granted, free of charge, to any person obtaining a copy
   8  * of this software and associated documentation files (the "Software"), to deal
   9  * in the Software without restriction, including without limitation the rights
  10  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  11  * copies of the Software, and to permit persons to whom the Software is
  12  * furnished to do so, subject to the following conditions:
  13  *
  14  * The above copyright notice and this permission notice shall be included in
  15  * all copies or substantial portions of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  18  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  20  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  21  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  22  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  23  * THE SOFTWARE.
  24  *
  25  **************************************************************************/
  26
  27 // File: vogl_bigint128.h
  28 #ifndef VOGL_BIGINT128_H
  29 #define VOGL_BIGINT128_H
  30
  31 #include "vogl_core.h"
  32 #include "vogl_rand.h"
  33
  34 namespace vogl
  35 {
  36     // Straightforward, portable 128-bit signed integer class.
  37     // Mostly intended to simplifying dealing with either signed or unsigned 64-bit values in a universal way.
  38     // Not particularly fast or flexible, it can't return 256 bits from 128*128 muls yet which is annoying.
  39     class bigint128 : public utils::relative_ops<bigint128>
  40     {
  41     public:
  42         enum
  43         {
  44             cMaxBYTEs = 16,
  45             cMaxQWORDs = 2,
  46             cMaxDWORDs = 4
  47         };
  48
  49         inline bigint128()
  50         {
  51         }
  52
  53         inline bigint128(int32 x)
  54         {
  55             *this = x;
  56         }
  57
  58         inline bigint128(int64_t x)
  59         {
  60             *this = x;
  61         }
  62
  63         inline bigint128(uint32 x)
  64         {
  65             *this = x;
  66         }
  67
  68         inline bigint128(uint64_t x)
  69         {
  70             *this = x;
  71         }
  72
  73         inline bigint128(uint64_t l, uint64_t h)
  74         {
  75             set(l, h);
  76         }
  77
  78         inline bigint128 &set(uint64_t l, uint64_t h)
  79         {
  80             m_v[0] = l;
  81             m_v[1] = h;
  82             return *this;
  83         }
  84
  85         bigint128(const void *p, uint type_size_in_bytes, bool type_is_signed)
  86         {
  87             set_from_ptr(p, type_size_in_bytes, type_is_signed);
  88         }
  89
  90         inline uint64_t get_qword(uint index) const
  91         {
  92             return m_v[math::open_range_check<uint>(index, cMaxQWORDs)];
  93         }
  94         inline bigint128 &set_qword(uint index, uint64_t val)
  95         {
  96             m_v[math::open_range_check<uint>(index, cMaxQWORDs)] = val;
  97             return *this;
  98         }
  99
 100         inline uint32 get_dword(uint index) const
 101         {
 102             VOGL_ASSERT(index < cMaxDWORDs);
 103             uint64_t x = get_qword(index >> 1);
 104             if (index & 1)
 105                 x >>= 32U;
 106             return static_cast<uint32>(x);
 107         }
 108
 109         inline bigint128 &set_dword(uint index, uint32 val)
 110         {
 111             VOGL_ASSERT(index < cMaxDWORDs);
 112             uint qword_index = index >> 1U;
 113             uint shift = (index & 1) ? 32U : 0U;
 114
 115             uint64_t mask = cUINT32_MAX;
 116             mask <<= shift;
 117
 118             m_v[qword_index] &= ~mask;
 119             m_v[qword_index] |= (static_cast<uint64_t>(val) << shift);
 120
 121             return *this;
 122         }
 123
 124         inline void clear()
 125         {
 126             m_v[0] = 0;
 127             m_v[1] = 0;
 128         }
 129
 130         inline bigint128 &operator=(int64_t x)
 131         {
 132             m_v[0] = x;
 133             m_v[1] = (x >> 63);
 134             return *this;
 135         }
 136
 137         inline bigint128 &operator=(int32 x)
 138         {
 139             *this = static_cast<int64_t>(x);
 140             return *this;
 141         }
 142
 143         inline bigint128 &operator=(uint64_t x)
 144         {
 145             m_v[0] = x;
 146             m_v[1] = 0;
 147             return *this;
 148         }
 149
 150         inline bigint128 &operator=(uint32 x)
 151         {
 152             *this = static_cast<uint64_t>(x);
 153             return *this;
 154         }
 155
 156         inline operator int8() const
 157         {
 158             return static_cast<int8>(m_v[0]);
 159         }
 160
 161         inline operator int16() const
 162         {
 163             return static_cast<int16>(m_v[0]);
 164         }
 165
 166         inline operator int32() const
 167         {
 168             return static_cast<int32>(m_v[0]);
 169         }
 170
 171         inline operator int64_t() const
 172         {
 173             return static_cast<int64_t>(m_v[0]);
 174         }
 175
 176         inline operator uint8() const
 177         {
 178             return static_cast<uint8>(m_v[0]);
 179         }
 180
 181         inline operator uint16() const
 182         {
 183             return static_cast<uint16>(m_v[0]);
 184         }
 185
 186         inline operator uint32() const
 187         {
 188             return static_cast<uint32>(m_v[0]);
 189         }
 190
 191         inline operator uint64_t() const
 192         {
 193             return static_cast<uint64_t>(m_v[0]);
 194         }
 195
 196         inline bool is_zero() const
 197         {
 198             return (m_v[0] == 0) && (m_v[1] == 0);
 199         }
 200
 201         inline operator bool() const
 202         {
 203             return !is_zero();
 204         }
 205
 206         // FIXME: not endian safe, assumes little endian
 207         bigint128 &set_from_ptr(const void *p, uint type_size_in_bytes, bool type_is_signed)
 208         {
 209             if (type_size_in_bytes == 0)
 210             {
 211                 clear();
 212                 return *this;
 213             }
 214             else if (type_size_in_bytes > cMaxBYTEs)
 215             {
 216                 VOGL_ASSERT_ALWAYS;
 217                 clear();
 218                 return *this;
 219             }
 220
 221             const uint8 *pSrc = reinterpret_cast<const uint8 *>(p);
 222             uint8 *pDst = reinterpret_cast<uint8 *>(m_v);
 223             memcpy(pDst, pSrc, type_size_in_bytes);
 224
 225             if (type_size_in_bytes < cMaxBYTEs)
 226             {
 227                 uint z = 0;
 228                 if (type_is_signed)
 229                 {
 230                     bool is_signed = (pSrc[type_size_in_bytes - 1] & 0x80) != 0;
 231                     z = is_signed ? 0xFF : 0;
 232                 }
 233                 memset(pDst + type_size_in_bytes, z, cMaxBYTEs - type_size_in_bytes);
 234             }
 235
 236             return *this;
 237         }
 238
 239         static inline void get_type_range(uint type_size_in_bytes, bool type_is_signed, bigint128 &min_val, bigint128 &max_val)
 240         {
 241             if ((!type_size_in_bytes) || (type_size_in_bytes > 16))
 242             {
 243                 VOGL_ASSERT_ALWAYS;
 244                 min_val = 0U;
 245                 max_val = 0U;
 246                 return;
 247             }
 248
 249             max_val = 1U;
 250             max_val <<= (type_size_in_bytes * 8U);
 251             --max_val;
 252
 253             min_val = 0U;
 254             if (type_is_signed)
 255             {
 256                 max_val = unsigned_shift_right(max_val, 1);
 257                 min_val = -max_val - 1U;
 258             }
 259         }
 260
 261         inline bool range_check(uint type_size_in_bytes, bool type_is_signed) const
 262         {
 263             if (!type_size_in_bytes)
 264             {
 265                 VOGL_ASSERT_ALWAYS;
 266                 return false;
 267             }
 268             else if (type_size_in_bytes > 16)
 269                 return true;
 270
 271             bigint128 min_val, max_val;
 272             get_type_range(type_size_in_bytes, type_is_signed, min_val, max_val);
 273             if (type_is_signed)
 274                 return (*this >= min_val) && (*this <= max_val);
 275             else
 276                 return unsigned_less_equal(max_val);
 277         }
 278
 279         inline bool is_negative() const
 280         {
 281             return static_cast<int64_t>(m_v[1]) < 0;
 282         }
 283
 284         inline bool is_positive() const
 285         {
 286             return static_cast<int64_t>(m_v[1]) >= 0;
 287         }
 288
 289         inline bigint128 operator+() const
 290         {
 291             return *this;
 292         }
 293
 294         inline bool operator==(const bigint128 &other) const
 295         {
 296             return (m_v[0] == other.m_v[0]) && (m_v[1] == other.m_v[1]);
 297         }
 298         inline bool operator==(int32 x) const
 299         {
 300             return *this == bigint128(x);
 301         }
 302         inline bool operator==(uint32 x) const
 303         {
 304             return *this == bigint128(x);
 305         }
 306         inline bool operator==(int64_t x) const
 307         {
 308             return *this == bigint128(x);
 309         }
 310         inline bool operator==(uint64_t x) const
 311         {
 312             return *this == bigint128(x);
 313         }
 314
 315         inline bool unsigned_less(const bigint128 &other) const
 316         {
 317             if (m_v[1] < other.m_v[1])
 318                 return true;
 319             else if (m_v[1] == other.m_v[1])
 320                 return m_v[0] < other.m_v[0];
 321
 322             return false;
 323         }
 324         inline bool unsigned_greater(const bigint128 &other) const
 325         {
 326             return other.unsigned_less(*this);
 327         }
 328         inline bool unsigned_less_equal(const bigint128 &other) const
 329         {
 330             return !other.unsigned_less(*this);
 331         }
 332         inline bool unsigned_greater_equal(const bigint128 &other) const
 333         {
 334             return !unsigned_less(other);
 335         }
 336
 337         inline bool operator<(const bigint128 &other) const
 338         {
 339             if (static_cast<int64_t>(m_v[1]) < static_cast<int64_t>(other.m_v[1]))
 340                 return true;
 341             else if (m_v[1] == other.m_v[1])
 342                 return m_v[0] < other.m_v[0];
 343
 344             return false;
 345         }
 346         inline bool operator<(int32 x) const
 347         {
 348             return *this < bigint128(x);
 349         }
 350         inline bool operator<(uint32 x) const
 351         {
 352             return *this < bigint128(x);
 353         }
 354         inline bool operator<(int64_t x) const
 355         {
 356             return *this < bigint128(x);
 357         }
 358         inline bool operator<(uint64_t x) const
 359         {
 360             return *this < bigint128(x);
 361         }
 362         inline bigint128 &operator-=(const bigint128 &other)
 363         {
 364             *this = *this - other;
 365             return *this;
 366         }
 367         inline bool operator<=(int32 x) const
 368         {
 369             return *this <= bigint128(x);
 370         }
 371         inline bool operator<=(uint32 x) const
 372         {
 373             return *this <= bigint128(x);
 374         }
 375         inline bool operator<=(int64_t x) const
 376         {
 377             return *this <= bigint128(x);
 378         }
 379         inline bool operator<=(uint64_t x) const
 380         {
 381             return *this <= bigint128(x);
 382         }
 383
 384         inline bool operator>(int32 x) const
 385         {
 386             return *this > bigint128(x);
 387         }
 388         inline bool operator>(uint32 x) const
 389         {
 390             return *this > bigint128(x);
 391         }
 392         inline bool operator>(int64_t x) const
 393         {
 394             return *this > bigint128(x);
 395         }
 396         inline bool operator>(uint64_t x) const
 397         {
 398             return *this > bigint128(x);
 399         }
 400
 401         inline bool operator>=(int32 x) const
 402         {
 403             return *this >= bigint128(x);
 404         }
 405         inline bool operator>=(uint32 x) const
 406         {
 407             return *this >= bigint128(x);
 408         }
 409         inline bool operator>=(int64_t x) const
 410         {
 411             return *this >= bigint128(x);
 412         }
 413         inline bool operator>=(uint64_t x) const
 414         {
 415             return *this >= bigint128(x);
 416         }
 417
 418         inline bigint128 operator~() const
 419         {
 420             bigint128 result;
 421             result.m_v[0] = ~m_v[0];
 422             result.m_v[1] = ~m_v[1];
 423             return result;
 424         }
 425
 426         inline bigint128 &operator++()
 427         {
 428             bool carry = (m_v[0] == cUINT64_MAX);
 429             m_v[0]++;
 430             m_v[1] += carry;
 431             return *this;
 432         }
 433
 434         inline bigint128 operator++(int)
 435         {
 436             bigint128 result(*this);
 437             bool carry = (m_v[0] == cUINT64_MAX);
 438             m_v[0]++;
 439             m_v[1] += carry;
 440             return result;
 441         }
 442
 443         inline bigint128 &operator--()
 444         {
 445             bool carry = (m_v[0] == 0);
 446             m_v[0]--;
 447             m_v[1] -= carry;
 448             return *this;
 449         }
 450
 451         inline bigint128 operator--(int)
 452         {
 453             bigint128 result(*this);
 454             bool carry = (m_v[0] == 0);
 455             m_v[0]--;
 456             m_v[1] -= carry;
 457             return result;
 458         }
 459
 460         inline bigint128 &operator|=(const bigint128 &other)
 461         {
 462             *this = *this | other;
 463             return *this;
 464         }
 465
 466         inline bigint128 &operator^=(const bigint128 &other)
 467         {
 468             *this = *this ^ other;
 469             return *this;
 470         }
 471
 472         inline bigint128 &operator&=(const bigint128 &other)
 473         {
 474             *this = *this & other;
 475             return *this;
 476         }
 477
 478         inline bigint128 &operator>>=(uint n)
 479         {
 480             *this = *this >> n;
 481             return *this;
 482         }
 483
 484         inline bigint128 &operator<<=(uint n)
 485         {
 486             *this = *this << n;
 487             return *this;
 488         }
 489
 490         inline bigint128 get_abs() const
 491         {
 492             return is_negative() ? (-(*this)) : (*this);
 493         }
 494
 495         friend inline bigint128 operator|(const bigint128 &lhs, const bigint128 &rhs)
 496         {
 497             return bigint128(lhs.m_v[0] | rhs.m_v[0], lhs.m_v[1] | rhs.m_v[1]);
 498         }
 499
 500         friend inline bigint128 operator^(const bigint128 &lhs, const bigint128 &rhs)
 501         {
 502             return bigint128(lhs.m_v[0] ^ rhs.m_v[0], lhs.m_v[1] ^ rhs.m_v[1]);
 503         }
 504
 505         friend inline bigint128 operator&(const bigint128 &lhs, const bigint128 &rhs)
 506         {
 507             return bigint128(lhs.m_v[0] & rhs.m_v[0], lhs.m_v[1] & rhs.m_v[1]);
 508         }
 509
 510         friend inline bigint128 operator+(const bigint128 &lhs, const bigint128 &rhs)
 511         {
 512             bigint128 result;
 513
 514             result.m_v[0] = lhs.m_v[0] + rhs.m_v[0];
 515
 516             // Being careful here to avoid known gcc overflow bugs
 517             uint64_t max_to_add_before_overflow = cUINT64_MAX - lhs.m_v[0];
 518             bool carry = rhs.m_v[0] > max_to_add_before_overflow;
 519
 520             result.m_v[1] = lhs.m_v[1] + rhs.m_v[1] + carry;
 521             return result;
 522         }
 523
 524         inline bigint128 &operator+=(const bigint128 &other)
 525         {
 526             *this = *this + other;
 527             return *this;
 528         }
 529
 530         friend inline bigint128 operator-(const bigint128 &lhs, const bigint128 &rhs)
 531         {
 532             bigint128 result;
 533
 534             result.m_v[0] = lhs.m_v[0] - rhs.m_v[0];
 535
 536             bool carry = rhs.m_v[0] > lhs.m_v[0];
 537
 538             result.m_v[1] = lhs.m_v[1] - rhs.m_v[1] - carry;
 539             return result;
 540         }
 541
 542         inline bigint128 operator-() const
 543         {
 544             bigint128 result;
 545             result.m_v[0] = ~m_v[0];
 546             result.m_v[1] = ~m_v[1];
 547
 548             if (result.m_v[0] != cUINT64_MAX)
 549                 result.m_v[0]++;
 550             else
 551             {
 552                 result.m_v[0] = 0;
 553                 result.m_v[1]++;
 554             }
 555             return result;
 556         }
 557
 558         friend inline bigint128 operator+(const bigint128 &lhs, int32 x)
 559         {
 560             return lhs + bigint128(x);
 561         }
 562         friend inline bigint128 operator+(const bigint128 &lhs, int64_t x)
 563         {
 564             return lhs + bigint128(x);
 565         }
 566         friend inline bigint128 operator+(const bigint128 &lhs, uint32 x)
 567         {
 568             return lhs + bigint128(x);
 569         }
 570         friend inline bigint128 operator+(const bigint128 &lhs, uint64_t x)
 571         {
 572             return lhs + bigint128(x);
 573         }
 574
 575         friend inline bigint128 operator-(const bigint128 &lhs, int32 x)
 576         {
 577             return lhs - bigint128(x);
 578         }
 579         friend inline bigint128 operator-(const bigint128 &lhs, int64_t x)
 580         {
 581             return lhs - bigint128(x);
 582         }
 583         friend inline bigint128 operator-(const bigint128 &lhs, uint32 x)
 584         {
 585             return lhs - bigint128(x);
 586         }
 587         friend inline bigint128 operator-(const bigint128 &lhs, uint64_t x)
 588         {
 589             return lhs - bigint128(x);
 590         }
 591
 592         // signed shift right
 593         friend inline bigint128 operator>>(const bigint128 &lhs, uint n)
 594         {
 595             if (!n)
 596                 return lhs;
 597             else if (n >= 128)
 598                 return lhs.is_negative() ? bigint128(-1) : bigint128(UINT64_C(0));
 599
 600             uint64_t l = lhs.m_v[0];
 601             uint64_t h = lhs.m_v[1];
 602
 603             if (n >= 64)
 604             {
 605                 l = h;
 606                 h = lhs.is_negative() ? cUINT64_MAX : 0;
 607                 n -= 64;
 608             }
 609
 610             // n can be 0-63 here
 611             if (n)
 612             {
 613                 l >>= n;
 614                 l |= ((h & ((UINT64_C(1) << n) - UINT64_C(1))) << (64 - static_cast<int>(n)));
 615                 h = static_cast<uint64_t>(static_cast<int64_t>(h) >> static_cast<int>(n));
 616             }
 617
 618             return bigint128(l, h);
 619         }
 620
 621         static inline bigint128 unsigned_shift_right(const bigint128 &lhs, uint n)
 622         {
 623             if (!n)
 624                 return lhs;
 625             else if (n >= 128)
 626                 return bigint128(UINT64_C(0));
 627
 628             uint64_t l = lhs.m_v[0];
 629             uint64_t h = lhs.m_v[1];
 630
 631             if (n >= 64)
 632             {
 633                 l = h;
 634                 h = 0;
 635                 n -= 64;
 636             }
 637
 638             // n can be 0-63 here
 639             if (n)
 640             {
 641                 l >>= n;
 642                 l |= ((h & ((UINT64_C(1) << n) - UINT64_C(1))) << (64 - static_cast<int>(n)));
 643                 h >>= n;
 644             }
 645
 646             return bigint128(l, h);
 647         }
 648
 649         // shift left
 650         friend inline bigint128 operator<<(const bigint128 &lhs, uint n)
 651         {
 652             if (!n)
 653                 return lhs;
 654             else if (n >= 128)
 655                 return bigint128(UINT64_C(0));
 656
 657             uint64_t l = lhs.m_v[0];
 658             uint64_t h = lhs.m_v[1];
 659
 660             if (n >= 64)
 661             {
 662                 h = l;
 663                 l = 0;
 664                 n -= 64;
 665             }
 666
 667             // n can be 0-63 here
 668             if (n)
 669             {
 670                 h <<= n;
 671                 h |= (l >> (64 - static_cast<int>(n)));
 672                 l <<= n;
 673             }
 674
 675             return bigint128(l, h);
 676         }
 677
 678         // unsigned op
 679         uint get_num_significant_bits() const
 680         {
 681             bigint128 d(*this);
 682             uint num = 0;
 683             while (d)
 684             {
 685                 num++;
 686                 d = d.unsigned_shift_right(d, 1);
 687             }
 688             return num;
 689         }
 690
 691         // true if overflowed (returns lower 128 bits of result)
 692         static bool unsigned_multiply(bigint128 x, bigint128 y, bigint128 &result)
 693         {
 694             // Shortcut if both are 32-bit.
 695             if ((x.m_v[1] | y.m_v[1]) == 0)
 696             {
 697                 if (((x.m_v[0] | y.m_v[0]) >> 32U) == 0)
 698                 {
 699                     result = x.m_v[0] * y.m_v[0];
 700                     return false;
 701                 }
 702             }
 703
 704             // Minor optimization: ensure x <= y.
 705             if (x > y)
 706             {
 707                 std::swap(x, y);
 708             }
 709
 710             bool overflow_flag = false;
 711
 712             bigint128 value(0U);
 713
 714             // Seems most portable to use unsigned mul of 32x32=64 bits as the underlying primitive.
 715             for (uint i = 0; i < cMaxDWORDs; ++i)
 716             {
 717                 const uint64_t m = static_cast<uint64_t>(x.get_dword(i));
 718                 if (!m)
 719                     continue;
 720                 for (uint j = 0; j < cMaxDWORDs; ++j)
 721                 {
 722                     uint k = i + j;
 723
 724                     uint64_t product = m * y.get_dword(j);
 725                     while (product)
 726                     {
 727                         if (k >= cMaxDWORDs)
 728                         {
 729                             overflow_flag = true;
 730                             break;
 731                         }
 732
 733                         product += value.get_dword(k);
 734                         value.set_dword(k, static_cast<uint32>(product));
 735                         product >>= 32U;
 736                         k++;
 737                     }
 738                 }
 739             }
 740             result = value;
 741
 742             return overflow_flag;
 743         }
 744
 745         // true if overflowed
 746         static bool signed_multiply(bigint128 a, bigint128 b, bigint128 &result)
 747         {
 748             bool a_is_negative = a.is_negative();
 749             bool b_is_negative = b.is_negative();
 750
 751             bool overflow_flag = unsigned_multiply(a.get_abs(), b.get_abs(), result);
 752
 753             if (a_is_negative != b_is_negative)
 754             {
 755                 if (result.unsigned_greater(bigint128(0, UINT64_C(0x8000000000000000))))
 756                     overflow_flag = true;
 757
 758                 result = -result;
 759             }
 760             else if (result.is_negative())
 761                 overflow_flag = true;
 762
 763             return overflow_flag;
 764         }
 765
 766         // true on success
 767         static bool unsigned_divide(const bigint128 &a, const bigint128 &b, bigint128 &q, bigint128 &r)
 768         {
 769             // Shortcut if both are 64-bit.
 770             if ((!a.m_v[1]) && (!b.m_v[1]))
 771             {
 772                 q = a.m_v[0] / b.m_v[0];
 773                 r = a.m_v[0] % b.m_v[0];
 774                 return true;
 775             }
 776
 777             int num_a = a.get_num_significant_bits();
 778             int num_b = b.get_num_significant_bits();
 779             if (!num_b)
 780             {
 781                 VOGL_ASSERT_ALWAYS;
 782                 q = 0;
 783                 r = 0;
 784                 return false;
 785             }
 786
 787             int cur_pos = num_a - num_b;
 788
 789             // Good old fashioned non-restoring division 1 bit at a time. Hope you're not busy!
 790             bigint128 c(a);
 791             bigint128 f(0);
 792             bigint128 k;
 793             if (cur_pos >= 0)
 794             {
 795                 k = b << static_cast<uint>(cur_pos);
 796
 797                 do
 798                 {
 799                     f <<= 1U;
 800
 801                     if (c.unsigned_greater_equal(k))
 802                     {
 803                         f |= UINT64_C(1);
 804                         c -= k;
 805                     }
 806
 807                     k = unsigned_shift_right(k, 1U);
 808                     cur_pos--;
 809                 } while (cur_pos >= 0);
 810             }
 811
 812             q = f;
 813             r = c;
 814
 815             VOGL_ASSERT(r.unsigned_less(b));
 816             return true;
 817         }
 818
 819         static bool signed_divide(const bigint128 &a, const bigint128 &b, bigint128 &q, bigint128 &r)
 820         {
 821             if (!b)
 822             {
 823                 VOGL_ASSERT_ALWAYS;
 824                 q = 0;
 825                 r = 0;
 826                 return false;
 827             }
 828
 829             bool a_is_negative = a.is_negative();
 830             bool b_is_negative = b.is_negative();
 831
 832             // abs()'s are safe here, the min 128-bit signed value can be represented as an unsigned 128-bit value
 833             unsigned_divide(a.get_abs(), b.get_abs(), q, r);
 834
 835             if (a_is_negative != b_is_negative)
 836                 q = -q;
 837
 838             if (a_is_negative)
 839                 r = -r;
 840
 841             return true;
 842         }
 843
 844         friend inline bigint128 operator*(const bigint128 &lhs, const bigint128 &rhs)
 845         {
 846             bigint128 result;
 847             signed_multiply(lhs, rhs, result);
 848             return result;
 849         }
 850         friend inline bigint128 operator*(const bigint128 &lhs, int32 x)
 851         {
 852             return lhs * bigint128(x);
 853         }
 854         friend inline bigint128 operator*(const bigint128 &lhs, int64_t x)
 855         {
 856             return lhs * bigint128(x);
 857         }
 858         friend inline bigint128 operator*(const bigint128 &lhs, uint32 x)
 859         {
 860             return lhs * bigint128(x);
 861         }
 862         friend inline bigint128 operator*(const bigint128 &lhs, uint64_t x)
 863         {
 864             return lhs * bigint128(x);
 865         }
 866
 867         bigint128 &operator*=(const bigint128 &other)
 868         {
 869             *this = *this * other;
 870             return *this;
 871         }
 872         bigint128 &operator*=(int32 other)
 873         {
 874             *this = *this * other;
 875             return *this;
 876         }
 877         bigint128 &operator*=(uint32 other)
 878         {
 879             *this = *this * other;
 880             return *this;
 881         }
 882         bigint128 &operator*=(int64_t other)
 883         {
 884             *this = *this * other;
 885             return *this;
 886         }
 887         bigint128 &operator*=(uint64_t other)
 888         {
 889             *this = *this * other;
 890             return *this;
 891         }
 892
 893         friend inline bigint128 operator/(const bigint128 &lhs, const bigint128 &rhs)
 894         {
 895             bigint128 q, r;
 896             signed_divide(lhs, rhs, q, r);
 897             return q;
 898         }
 899         friend inline bigint128 operator/(const bigint128 &lhs, int32 x)
 900         {
 901             return lhs / bigint128(x);
 902         }
 903         friend inline bigint128 operator/(const bigint128 &lhs, int64_t x)
 904         {
 905             return lhs / bigint128(x);
 906         }
 907         friend inline bigint128 operator/(const bigint128 &lhs, uint32 x)
 908         {
 909             return lhs / bigint128(x);
 910         }
 911         friend inline bigint128 operator/(const bigint128 &lhs, uint64_t x)
 912         {
 913             return lhs / bigint128(x);
 914         }
 915
 916         bigint128 &operator/=(const bigint128 &other)
 917         {
 918             *this = *this / other;
 919             return *this;
 920         }
 921         bigint128 &operator/=(int32 other)
 922         {
 923             *this = *this / other;
 924             return *this;
 925         }
 926         bigint128 &operator/=(uint32 other)
 927         {
 928             *this = *this / other;
 929             return *this;
 930         }
 931         bigint128 &operator/=(int64_t other)
 932         {
 933             *this = *this / other;
 934             return *this;
 935         }
 936         bigint128 &operator/=(uint64_t other)
 937         {
 938             *this = *this / other;
 939             return *this;
 940         }
 941
 942         friend inline bigint128 operator%(const bigint128 &lhs, const bigint128 &rhs)
 943         {
 944             bigint128 q, r;
 945             signed_divide(lhs, rhs, q, r);
 946             return r;
 947         }
 948         friend inline bigint128 operator%(const bigint128 &lhs, int32 x)
 949         {
 950             return lhs % bigint128(x);
 951         }
 952         friend inline bigint128 operator%(const bigint128 &lhs, int64_t x)
 953         {
 954             return lhs % bigint128(x);
 955         }
 956         friend inline bigint128 operator%(const bigint128 &lhs, uint32 x)
 957         {
 958             return lhs % bigint128(x);
 959         }
 960         friend inline bigint128 operator%(const bigint128 &lhs, uint64_t x)
 961         {
 962             return lhs % bigint128(x);
 963         }
 964
 965         bigint128 &operator%=(const bigint128 &other)
 966         {
 967             *this = *this % other;
 968             return *this;
 969         }
 970         bigint128 &operator%=(int32 other)
 971         {
 972             *this = *this % other;
 973             return *this;
 974         }
 975         bigint128 &operator%=(uint32 other)
 976         {
 977             *this = *this % other;
 978             return *this;
 979         }
 980         bigint128 &operator%=(int64_t other)
 981         {
 982             *this = *this % other;
 983             return *this;
 984         }
 985         bigint128 &operator%=(uint64_t other)
 986         {
 987             *this = *this % other;
 988             return *this;
 989         }
 990
 991         static bigint128 get_random(random &r)
 992         {
 993             return bigint128(r.urand64(), r.urand64());
 994         }
 995
 996         static bigint128 get_smallest()
 997         {
 998             return bigint128(0, cINT64_MIN);
 999         }
1000
1001         static bigint128 get_largest()
1002         {
1003             return bigint128(cUINT64_MAX, cINT64_MAX);
1004         }
1005
1006         static bigint128 unsigned_get_largest()
1007         {
1008             return bigint128(cUINT64_MAX, cUINT64_MAX);
1009         }
1010
1011         static bigint128 get_all_zeros()
1012         {
1013             return bigint128(0U, 0U);
1014         }
1015
1016         static bigint128 get_all_ones()
1017         {
1018             return bigint128(cUINT64_MAX, cUINT64_MAX);
1019         }
1020
1021     private:
1022         uint64_t m_v[2];
1023     };
1024
1025     VOGL_DEFINE_BITWISE_COPYABLE(bigint128);
1026
1027     inline void bigint128_check(bool success, uint &num_failures)
1028     {
1029         if (!success)
1030         {
1031             num_failures++;
1032             VOGL_VERIFY(0);
1033         }
1034     }
1035
1036     inline bool bigint128_test(uint seed = 5)
1037     {
1038         VOGL_NOTE_UNUSED(seed);
1039         VOGL_ASSUME(sizeof(bigint128) == sizeof(uint64_t) * 2U);
1040
1041         uint num_failures = 0;
1042
1043         bigint128 zero(0U);
1044         bigint128 one(1);
1045         bigint128 neg_one(-1);
1046
1047 #define BIGINT128_CHECK(x) bigint128_check(x, num_failures)
1048         BIGINT128_CHECK(neg_one.get_qword(0) == cUINT64_MAX);
1049         BIGINT128_CHECK(neg_one.get_qword(1) == cUINT64_MAX);
1050         BIGINT128_CHECK(neg_one + 2 == 1);
1051         BIGINT128_CHECK(zero.is_zero());
1052         BIGINT128_CHECK(one.is_positive() && !one.is_zero() && !one.is_negative());
1053         BIGINT128_CHECK(neg_one.is_negative() && !neg_one.is_zero() && !neg_one.is_positive());
1054         BIGINT128_CHECK(one > zero);
1055         BIGINT128_CHECK(!(zero > one));
1056         BIGINT128_CHECK(zero < one);
1057         BIGINT128_CHECK(one > neg_one);
1058         BIGINT128_CHECK(neg_one < zero);
1059         BIGINT128_CHECK(zero > neg_one);
1060         BIGINT128_CHECK(neg_one < zero);
1061         BIGINT128_CHECK(bigint128(zero)++ == zero);
1062         BIGINT128_CHECK(++bigint128(zero) == one);
1063         BIGINT128_CHECK(bigint128(one)-- == one);
1064         BIGINT128_CHECK(--bigint128(one) == zero);
1065         BIGINT128_CHECK(--bigint128(zero) < zero);
1066         BIGINT128_CHECK(--bigint128(zero) == -1);
1067         BIGINT128_CHECK(++bigint128(zero) == 1);
1068         BIGINT128_CHECK(++bigint128(cUINT64_MAX) == bigint128(0, 1));
1069         BIGINT128_CHECK(--bigint128(0, 1) == cUINT64_MAX);
1070         BIGINT128_CHECK(++bigint128(cUINT64_MAX, cUINT64_MAX) == zero);
1071         BIGINT128_CHECK(bigint128(0, 0) - bigint128(1, 0) == bigint128(-1));
1072         BIGINT128_CHECK(bigint128(1, 0) + bigint128(cUINT64_MAX, cUINT64_MAX) == bigint128(0));
1073
1074         BIGINT128_CHECK(bigint128(0, 1) < bigint128(0, 2));
1075         BIGINT128_CHECK(bigint128(0, 2) > bigint128(0, 1));
1076         BIGINT128_CHECK(bigint128(0, cUINT64_MAX) < -1);
1077         BIGINT128_CHECK(bigint128(0, cUINT64_MAX) > bigint128::get_smallest());
1078
1079         BIGINT128_CHECK(bigint128::get_smallest() < bigint128::get_largest());
1080         BIGINT128_CHECK(bigint128::get_largest() > bigint128::get_smallest());
1081         BIGINT128_CHECK(bigint128::get_smallest() < 0);
1082         BIGINT128_CHECK(bigint128::get_largest() > 0);
1083         BIGINT128_CHECK(bigint128::get_smallest() < (bigint128::get_smallest() + 1));
1084         BIGINT128_CHECK(bigint128::get_largest() > (bigint128::get_largest() - 1));
1085
1086         BIGINT128_CHECK(bigint128::get_smallest() < bigint128(0, cINT64_MAX));
1087         BIGINT128_CHECK(bigint128::get_largest() > bigint128(0, cUINT64_MAX));
1088
1089         BIGINT128_CHECK((bigint128::get_smallest() - 1) == bigint128::get_largest());
1090         BIGINT128_CHECK((bigint128::get_largest() + 1) == bigint128::get_smallest());
1091
1092         BIGINT128_CHECK((bigint128::get_smallest() - 2) == (bigint128::get_largest() - 1));
1093         BIGINT128_CHECK((bigint128::get_largest() + 2) == (bigint128::get_smallest() + 1));
1094
1095         BIGINT128_CHECK(~bigint128(0) == bigint128(-1));
1096         BIGINT128_CHECK(~~bigint128(0) == bigint128(0));
1097         BIGINT128_CHECK((bigint128(0xDEADBEEF, 0xDEADBEEF) ^ bigint128(0xDEADBEEF, 0xDEADBEEF) ^ bigint128(0xDEADBEEF, 0xDEADBEEF) ^ bigint128(0xDEADBEEF, 0xDEADBEEF)) == 0);
1098
1099         bigint128 zz(bigint128(0xDEADBEEF, 0xDEADBEEF));
1100         BIGINT128_CHECK((zz ^= zz) == 0);
1101         BIGINT128_CHECK((bigint128(0) | bigint128(1)) == 1);
1102
1103         BIGINT128_CHECK(bigint128(0, 1) - bigint128(1) == cUINT64_MAX);
1104         BIGINT128_CHECK(bigint128(0, 50) - bigint128(0, 50) == 0);
1105         BIGINT128_CHECK(bigint128(0, 50) - bigint128(1, 50) == -1);
1106         BIGINT128_CHECK(bigint128(-1) - bigint128(-1) == 0);
1107         BIGINT128_CHECK(bigint128(-1) + bigint128(1) == 0);
1108
1109         BIGINT128_CHECK(bigint128(0).range_check(1, false));
1110         BIGINT128_CHECK(bigint128(UINT64_C(0xFF)).range_check(1, false));
1111
1112         BIGINT128_CHECK(bigint128(0).range_check(2, false));
1113         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFF)).range_check(2, false));
1114
1115         BIGINT128_CHECK(bigint128(0).range_check(3, false));
1116         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFF)).range_check(3, false));
1117
1118         BIGINT128_CHECK(bigint128(0).range_check(4, false));
1119         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFF)).range_check(4, false));
1120
1121         BIGINT128_CHECK(bigint128(0).range_check(5, false));
1122         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFFFF)).range_check(5, false));
1123
1124         BIGINT128_CHECK(bigint128(0).range_check(6, false));
1125         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFFFFFF)).range_check(6, false));
1126
1127         BIGINT128_CHECK(bigint128(0).range_check(7, false));
1128         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFFFFFFFF)).range_check(7, false));
1129
1130         BIGINT128_CHECK(bigint128(0).range_check(8, false));
1131         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFFFFFFFFFF)).range_check(8, false));
1132
1133         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFFFFFFFF80), cUINT64_MAX).range_check(1, true));
1134         BIGINT128_CHECK(!bigint128(UINT64_C(0xFFFFFFFFFFFFFF7F), cUINT64_MAX).range_check(1, true));
1135         BIGINT128_CHECK(bigint128(UINT64_C(0x000000000000007F), 0).range_check(1, true));
1136         BIGINT128_CHECK(!bigint128(UINT64_C(0x0000000000000080), 0).range_check(1, true));
1137
1138         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFFFFFF8000), cUINT64_MAX).range_check(2, true));
1139         BIGINT128_CHECK(!bigint128(UINT64_C(0xFFFFFFFFFFFF7FFF), cUINT64_MAX).range_check(2, true));
1140         BIGINT128_CHECK(bigint128(UINT64_C(0x0000000000007FFF), 0).range_check(2, true));
1141         BIGINT128_CHECK(!bigint128(UINT64_C(0x0000000000008000), 0).range_check(2, true));
1142
1143         BIGINT128_CHECK(bigint128(UINT64_C(0xFFFFFFFF80000000), cUINT64_MAX).range_check(4, true));
1144         BIGINT128_CHECK(!bigint128(UINT64_C(0xFFFFFFFF7FFFFFFF), cUINT64_MAX).range_check(4, true));
1145         BIGINT128_CHECK(bigint128(UINT64_C(0x000000007FFFFFFF), 0).range_check(4, true));
1146         BIGINT128_CHECK(!bigint128(UINT64_C(0x0000000080000000), 0).range_check(4, true));
1147
1148         BIGINT128_CHECK(bigint128(UINT64_C(0x8000000000000000), cUINT64_MAX).range_check(8, true));
1149         BIGINT128_CHECK(!bigint128(UINT64_C(0x7FFFFFFFFFFFFFFF), cUINT64_MAX).range_check(8, true));
1150         BIGINT128_CHECK(!(bigint128(UINT64_C(0x7FFFFFFFFFFFFFFF), cUINT64_MAX) - 1000).range_check(8, true));
1151         BIGINT128_CHECK(bigint128(UINT64_C(0x7FFFFFFFFFFFFFFF), 0).range_check(8, true));
1152         BIGINT128_CHECK(!bigint128(UINT64_C(0x8000000000000000), 0).range_check(8, true));
1153         BIGINT128_CHECK(!(bigint128(UINT64_C(0x8000000000000000), 0) + 1000).range_check(8, true));
1154
1155         BIGINT128_CHECK(bigint128(cUINT64_MAX, UINT64_C(0x8000000000000000)).range_check(16, true));
1156         BIGINT128_CHECK(bigint128(cUINT64_MAX, UINT64_C(0x7FFFFFFFFFFFFFFF)).range_check(16, true));
1157
1158         BIGINT128_CHECK(bigint128(0, 0).range_check(16, false));
1159         BIGINT128_CHECK(bigint128(cUINT64_MAX, 0).range_check(16, false));
1160         BIGINT128_CHECK(bigint128(cUINT64_MAX, cUINT64_MAX).range_check(16, false));
1161
1162         {
1163             bigint128 x(0U);
1164             BIGINT128_CHECK(x == 0);
1165             BIGINT128_CHECK(x.is_zero());
1166             BIGINT128_CHECK(!(x < 0));
1167             BIGINT128_CHECK(!(x > 0));
1168         }
1169
1170 #if defined(VOGL_PLATFORM_PC_X64) && defined(__GNUC__)
1171         vogl::random r;
1172         r.seed(seed);
1173         for (uint i = 0; i < 1000000000; i++)
1174         {
1175             __int128 x(0), y(0), z(0);
1176
1177             uint n0 = r.irand_inclusive(1, 16);
1178             for (uint j = 0; j < n0; j++)
1179                 ((uint8 *)&x)[j] = r.urand32() >> 8;
1180             if (n0 < 16)
1181             {
1182                 if (r.urand32() & 1)
1183                     x = -x;
1184             }
1185
1186             uint n1 = r.irand_inclusive(1, 16);
1187             for (uint j = 0; j < n1; j++)
1188                 ((uint8 *)&y)[j] = r.urand32() >> 8;
1189             if (n1 < 16)
1190             {
1191                 if (r.urand32() & 1)
1192                     y = -y;
1193             }
1194
1195             if (!r.irand(0, 100))
1196                 x = 0;
1197             if (!r.irand(0, 100))
1198                 y = 0;
1199
1200             uint op = r.irand_inclusive(0, 24);
1201
1202             bigint128 bx((uint64_t)x, (uint64_t)(x >> 64U));
1203             bigint128 by((uint64_t)y, (uint64_t)(y >> 64U));
1204             bigint128 bz(0);
1205
1206             if ((op >= 21) && (op <= 24))
1207             {
1208                 if (!r.irand(0, 30))
1209                 {
1210                     bx = bigint128::get_smallest();
1211                     memcpy(&x, &bx, sizeof(x));
1212                 }
1213                 if (!r.irand(0, 30))
1214                 {
1215                     bx = bigint128::get_largest();
1216                     memcpy(&x, &bx, sizeof(x));
1217                 }
1218                 if (!r.irand(0, 30))
1219                 {
1220                     by = bigint128::get_smallest();
1221                     memcpy(&y, &by, sizeof(y));
1222                 }
1223                 if (!r.irand(0, 30))
1224                 {
1225                     by = bigint128::get_largest();
1226                     memcpy(&y, &by, sizeof(y));
1227                 }
1228             }
1229
1230             uint shift = r.irand_inclusive(0, 127);
1231
1232             switch (op)
1233             {
1234                 case 0:
1235                 {
1236                     z = x + y;
1237                     bz = bx + by;
1238                     break;
1239                 }
1240                 case 1:
1241                 {
1242                     z = x - y;
1243                     bz = bx - by;
1244                     break;
1245                 }
1246                 case 2:
1247                 {
1248                     z = x < y;
1249                     bz = bx < by;
1250                     break;
1251                 }
1252                 case 3:
1253                 {
1254                     z = x <= y;
1255                     bz = bx <= by;
1256                     break;
1257                 }
1258                 case 4:
1259                 {
1260                     z = x > y;
1261                     bz = bx > by;
1262                     break;
1263                 }
1264                 case 5:
1265                 {
1266                     z = x >= y;
1267                     bz = bx >= by;
1268                     break;
1269                 }
1270                 case 6:
1271                 {
1272                     z = y ^ x;
1273                     bz = bx ^ by;
1274                     break;
1275                 }
1276                 case 7:
1277                 {
1278                     z = y | x;
1279                     bz = bx | by;
1280                     break;
1281                 }
1282                 case 8:
1283                 {
1284                     z = x >> shift;
1285                     bz = bx >> shift;
1286                     break;
1287                 }
1288                 case 9:
1289                 {
1290                     z = x << shift;
1291                     bz = bx << shift;
1292                     break;
1293                 }
1294                 case 10:
1295                 {
1296                     z = x & y;
1297                     bz = bx & by;
1298                     break;
1299                 }
1300                 case 11:
1301                 {
1302                     x += y;
1303                     z = x;
1304                     bx += by;
1305                     bz = bx;
1306                     break;
1307                 }
1308                 case 12:
1309                 {
1310                     x -= y;
1311                     z = x;
1312                     bx -= by;
1313                     bz = bx;
1314                     break;
1315                 }
1316                 case 13:
1317                 {
1318                     z = x--;
1319                     bz = bx--;
1320                     break;
1321                 }
1322                 case 14:
1323                 {
1324                     z = --x;
1325                     bz = --bx;
1326                     break;
1327                 }
1328                 case 15:
1329                 {
1330                     z = ++x;
1331                     bz = ++bx;
1332                     break;
1333                 }
1334                 case 16:
1335                 {
1336                     z = x++;
1337                     bz = bx++;
1338                     break;
1339                 }
1340
1341                 case 17:
1342                 {
1343                     z = (__uint128_t)x < (__uint128_t)y;
1344                     bz = bx.unsigned_less(by);
1345                     break;
1346                 }
1347                 case 18:
1348                 {
1349                     z = (__uint128_t)x <= (__uint128_t)y;
1350                     bz = bx.unsigned_less_equal(by);
1351                     break;
1352                 }
1353                 case 19:
1354                 {
1355                     z = (__uint128_t)x > (__uint128_t)y;
1356                     bz = bx.unsigned_greater(by);
1357                     break;
1358                 }
1359                 case 20:
1360                 {
1361                     z = (__uint128_t)x >= (__uint128_t)y;
1362                     bz = bx.unsigned_greater_equal(by);
1363                     break;
1364                 }
1365                 case 21:
1366                 {
1367                     if (y)
1368                     {
1369                         z = (__uint128_t)x / (__uint128_t)y;
1370                         __int128 m = (__uint128_t)x % (__uint128_t)y;
1371
1372                         bigint128 r;
1373                         VOGL_VERIFY(bigint128::unsigned_divide(bx, by, bz, r));
1374
1375                         BIGINT128_CHECK((uint64_t)m == r.get_qword(0));
1376                         BIGINT128_CHECK((uint64_t)(m >> 64) == r.get_qword(1));
1377                     }
1378                     break;
1379                 }
1380                 case 22:
1381                 {
1382                     if (y)
1383                     {
1384                         z = x / y;
1385                         __int128 m = x % y;
1386
1387                         bigint128 r;
1388                         VOGL_VERIFY(bigint128::signed_divide(bx, by, bz, r));
1389                         bigint128 nr(-r);
1390                         VOGL_NOTE_UNUSED(nr);
1391
1392                         BIGINT128_CHECK((uint64_t)m == r.get_qword(0));
1393                         BIGINT128_CHECK((uint64_t)(m >> 64) == r.get_qword(1));
1394                     }
1395                     break;
1396                 }
1397                 case 23:
1398                 {
1399                     z = (__uint128_t)x * (__uint128_t)y;
1400                     bigint128::unsigned_multiply(bx, by, bz);
1401                     break;
1402                 }
1403                 case 24:
1404                 {
1405                     z = x * y;
1406                     bigint128::signed_multiply(bx, by, bz);
1407                     break;
1408                 }
1409
1410                 // TODO: unsigned comps, signed/unsigned division
1411                 default:
1412                     VOGL_VERIFY(0);
1413                     break;
1414             }
1415
1416             BIGINT128_CHECK((uint64_t)z == bz.get_qword(0));
1417             BIGINT128_CHECK((uint64_t)(z >> 64) == bz.get_qword(1));
1418         }
1419 #endif
1420
1421 #undef BIGINT128_CHECK
1422
1423         return !num_failures;
1424     }
1425
1426 } // namespace vogl
1427
1428 #endif // VOGL_BIGINT128_H