/* int128.c - 128-bit integer arithmetic for C++, by Robert Munafo (and copied to rhtf/.../int128.c.txt by proj/.../MCS.per) LICENSE This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . The license (GPL version 3) for this source code is at: mrob.com/pub/ries/GPL-3.txt INSTRUCTIONS To add 128-bit integers to an existing program: Pass -DI128_COMPAT_PLUS_EQUAL to the compiler (unless on a really old compiler) Rename the program to "foo.cxx", recompile, deal with error messages Include .../int128.c Declare variables as s128 REVISION HISTORY 20040416 Fix bug in div1: Its answer was one too low when remainder is 0 20070812 Add +=, etc. operators 20070818 Add demotion to long long, and lots of operators that take s64, u64 or s32 on the RHS 20070820 Massively faster multiply algorithm 20070822 Add conversion from double 20070825 Fix bug in operator - (const s128_o &, s32) 20090619 Fix bug in div1: Its answer was one too low when quotient is a power of 2 and remainder is 0 (apparently a different case of the bug I fixed on 20040416) 20091122 Fix bug in mult1: Sign handling trickery failed on multiplying -1 * 0. For now it's back to the simple "treat it as unsigned" method. 20130410 Add GPL license header 20170215 Add #define switch I128_COMPAT_PLUS_EQUAL to avoid errors "'s128_o::s128_o' names the constructor, not the type" 20191019 Avoid comparisons between signed and unsigned, to fix bugs that appeared when using a compiler noncompliant with C language specification 6.3.1.8. (These bugs caused incorrect results when using operator+ or operator+= with s32 or s64 on the RHS) 20191022 s128_str now properly handles minimum value -(2^127) 20201029 Add "return lhs;" in the assignment operators (+=, /=, etc.) */ #ifndef _int128_c_ #define _int128_c_ /* Compile line should supply -I/Users/munafo/shared/proj/include or equivalent */ #include "rpmtypes.h" #include "int128.h" #include #include // Assignment and Assignment-Conversion operators s128_o :: s128_o ( int x ) { this->low = x; if (x < 0) { this->high = -1; } else { this->high = 0; } } s128_o :: s128_o ( s64 x ) { this->low = x; if (x < 0) { this->high = -1; } else { this->high = 0; } } s128_o :: s128_o ( double x ) { u64 t, m, h, l; if (x < -1.7014118346046e38) { // overflow negative this->high = HIGHBIT; this->low = 0; } else if (x < -9.2233720368547e18) { // take the 54 mantissa bits and shift into position t = *((u64 *) &x); m = (t & BITS51_0) | BIT_52; t = (t & BITS62_0) >> 52; // if x is 1.5 * 2^1, t will be 1024 // if x is 1.5 * 2^52, t will be 1024+51 = 1075 t = t - 1075; if (t > 64) { l = 0; h = m << (t-64); } else { l = m << t; h = m >> (64 - t); } this->low = ~l; this->high = ~h; } else if (x < 9.2233720368547e18) { // it will fit in a u64 this->low = ((u64) x); this->high = ((x<0) ? -1 : 0); } else if (x < 1.7014118346046e38) { // take the 54 mantissa bits and shift into position t = *((u64 *) &x); m = (t & BITS51_0) | BIT_52; t = (t & BITS62_0) >> 52; // if x is 1.5 * 2^1, t will be 1024 // if x is 1.5 * 2^52, t will be 1024+51 = 1075 t = t - 1075; if (t > 64) { this->low = 0; this->high = m << (t-64); } else { this->low = m << t; this->high = m >> (64 - t); } } else { // overflow positive this->high = BITS62_0; this->low = ALLBITS; } } // Infix addition s128_o operator + ( const s128_o & lhs, const s128_o & rhs ) { s128_o res; res.high = lhs.high + rhs.high; res.low = lhs.low + rhs.low; if (res.low < rhs.low) { (res.high)++; } return res; } /* 20191019 Some compilers that do not implement C spec 6.3.1.8 regarding integers of different types in binary operators; changes noted below */ s128_o operator + ( const s128_o & lhs, s64 rhs ) { s128_o res; if (rhs < 0) { res.high = lhs.high - 1; res.low = lhs.low + rhs; if (res.low < lhs.low) { // 20191019 was "(res.low lhs.low) { // borrow (res.high)--; } return res; } s128_o operator - ( const s128_o & lhs, s32 rhs ) { s128_o res; s128_o b; b.low = rhs; b.high = (rhs < 0) ? -1 : 0; res.high = lhs.high - b.high; res.low = lhs.low - b.low; if (res.low > lhs.low) { (res.high)--; } return res; } // Unary minus (negate) s128_o operator - ( const s128_o & x ) { s128_o res; res.high = ~(x.high); res.low = ~(x.low); res.low += 1; if (res.low == 0) { res.high += 1; } return res; } // Support routines for multiply, divide and modulo s128 s128_shr(s128_o x) { s128_o rv; // printf("%.16llX %016llX >> ", x.high, x.low); rv.low = (x.low >> 1) | (x.high << 63); rv.high = x.high >> 1; // printf("%.16llX %016llX\n", rv.high, rv.low); return rv; } s128 s128_shl(s128_o x) { s128_o rv; rv.high = (x.high << 1) | (x.low >> 63); rv.low = x.low << 1; return rv; } // Multiplication // old: 28.7 new: 6.3 /* Enable both of these defines to get error-checking */ #define MULT1_OLD 0 #define MULT1_NEW 1 static int m1err = 0; s128_o mult1(s128_o xi, s128_o yi) { #if MULT1_OLD s128_o x, y, rv1; int s; #endif #if MULT1_NEW s128_o rv2; u64 acc, ac2, carry, o1, o2; u64 a, b, c, d, e, f, g, h, ah, bh, ch, dh, bg, cg, dg, cf, df, de; #endif #if MULT1_OLD s = 1; x = xi; y = yi; if (x < ((s128_o) 0)) { s = - s; x = - x; } if (y < ((s128_o) 0)) { s = - s; y = - y; } rv1 = 0; while(y != ((s128_o) 0)) { if (y.low & 1) { rv1 = rv1 + x; } x = s128_shl(x); y = s128_shr(y); } if (s < 0) { rv1 = - rv1; } #endif #if MULT1_NEW /* x a b c d y e f g h --------------- ah bh ch dh bg cg dg cf df de -------------------------- -o2-- -o1-- */ d = xi.low & LO_WORD; c = (xi.low & HI_WORD) >> 32LL; b = xi.high & LO_WORD; a = (xi.high & HI_WORD) >> 32LL; h = yi.low & LO_WORD; g = (yi.low & HI_WORD) >> 32LL; f = yi.high & LO_WORD; e = (yi.high & HI_WORD) >> 32LL; acc = d * h; o1 = acc & LO_WORD; acc >>= 32LL; carry = 0; ac2 = acc + c * h; if (ac2 < acc) { carry++; } acc = ac2 + d * g; if (acc < ac2) { carry++; } rv2.low = o1 | (acc << 32LL); ac2 = (acc >> 32LL) | (carry << 32LL); carry = 0; acc = ac2 + b * h; if (acc < ac2) { carry++; } ac2 = acc + c * g; if (ac2 < acc) { carry++; } acc = ac2 + d * f; if (acc < ac2) { carry++; } o2 = acc & LO_WORD; ac2 = (acc >> 32LL) | (carry << 32LL); acc = ac2 + a * h; ac2 = acc + b * g; acc = ac2 + c * f; ac2 = acc + d * e; rv2.high = (ac2 << 32LL) | o2; #if MULT1_OLD if ((rv1.high != rv2.high) || (rv1.low != rv2.low)) { printf("m1 err1 %16llX %016llX * %16llX %016llX\n", xi.high, xi.low, yi.high, yi.low); printf("abcd %8llX %08llX %08llX %08llX\n", a, b, c, d); printf("efgh %8llX %08llX %08llX %08llX\n", e, f, g, h); printf(" -> %16llX %016llX\n", rv2.high, rv2.low); m1err++; if (m1err > 10) { exit(-1); } } #endif #endif #if MULT1_NEW return rv2; #else return rv1; #endif } s128_o operator * ( const s128_o & lhs, const s128_o & rhs ) { s128_o rv; rv = mult1(lhs, rhs); return rv; } s128_o operator * ( const s128_o & lhs, u64 rhs ) { s128_o t; s128_o rv; t = rhs; rv = mult1(lhs, t); return rv; } s128_o operator * ( const s128_o & lhs, s64 rhs ) { s128_o t; s128_o rv; t = rhs; rv = mult1(lhs, t); return rv; } s128_o operator * ( const s128_o & lhs, long int rhs ) { s128_o t; s128_o rv; t = ((s64)rhs); rv = mult1(lhs, t); return rv; } s128_o operator * ( const s128_o & lhs, s32 rhs ) { s128_o t; s128_o rv; t = rhs; rv = mult1(lhs, t); return rv; } // Division s128_o div1(s128_o x, s128_o d, s128_o *r) { int s; s128_o d1, p2, rv; //printf("divide %.16llX %016llX / %.16llX %016llX\n", x.high, x.low, d.high, d.low); /* check for divide by zero */ if ((d.low == 0) && (d.high == 0)) { rv.low = x.low / d.low; /* This will cause runtime error */ } s = 1; if (x < ((s128_o) 0)) { // notice that MININT will be unchanged, this is used below. s = - s; x = - x; } if (d < ((s128_o) 0)) { s = - s; d = - d; } if (d == ((s128_o) 1)) { /* This includes the overflow case MININT/-1 */ rv = x; x = 0; } else if (x < ((s128_o) d)) { /* x < d, so quotient is 0 and x is remainder */ rv = 0; } else { rv = 0; /* calculate biggest power of 2 times d that's <= x */ p2 = 1; d1 = d; x = x - d1; while(x >= d1) { x = x - d1; d1 = d1 + d1; p2 = p2 + p2; } x = x + d1; while(p2 != ((s128_o) 0)) { //printf("x %.16llX %016llX d1 %.16llX %016llX\n", x.high, x.low, d1.high, d1.low); if (x >= d1) { x = x - d1; rv = rv + p2; //printf("`.. %.16llX %016llX\n", rv.high, rv.low); } p2 = s128_shr(p2); d1 = s128_shr(d1); } /* whatever is left in x is the remainder */ } /* Put sign in result */ if (s < 0) { rv = - rv; } /* return remainder if they asked for it */ if (r) { *r = x; } return rv; } s128_o operator / ( const s128_o & x, const s128_o & d ) { s128_o rv; rv = div1(x, d, 0); return rv; } s128_o operator / ( const s128_o & x, u64 d ) { s128_o t; s128_o rv; t.high = 0; t.low = d; rv = div1(x, t, 0); return rv; } // Comparison operators int operator < ( const s128_o & lhs, const s128_o & rhs ) { if (lhs.high < rhs.high) return 1; if (rhs.high < lhs.high) return 0; // high components are equal if (lhs.low < rhs.low) return 1; return 0; } int operator < ( const s128_o & lhs, s32 rhs ) { s128_o r; r = rhs; if (lhs.high < r.high) return 1; if (r.high < lhs.high) return 0; // high components are equal if (lhs.low < r.low) return 1; return 0; } int operator <= ( const s128_o & lhs, const s128_o & rhs ) { if ((lhs.high == rhs.high) && (lhs.low == rhs.low)) return 1; if (lhs.high < rhs.high) return 1; if (rhs.high < lhs.high) return 0; // high components are equal if (lhs.low < rhs.low) return 1; return 0; } int operator <= ( const s128_o & lhs, s32 rhs ) { s128_o t; t = rhs; if ((lhs.high == t.high) && (lhs.low == t.low)) return 1; if (lhs.high < t.high) return 1; if (t.high < lhs.high) return 0; // high components are equal if (lhs.low < t.low) return 1; return 0; } int operator == ( const s128_o & lhs, const s128_o & rhs ) { if (lhs.high != rhs.high) return 0; if (lhs.low != rhs.low) return 0; return 1; } int operator == ( const s128_o & lhs, s32 rhs ) { s128_o t; t = rhs; return ((lhs.high == t.high) && (lhs.low == t.low)); } int operator != ( const s128_o & lhs, const s128_o & rhs ) { if (lhs.high != rhs.high) return 1; if (lhs.low != rhs.low) return 1; return 0; } int operator != ( const s128_o & lhs, s32 rhs ) { s128_o t; t = rhs; return ( (lhs.high != t.high) || (lhs.low != t.low) ); } int operator > ( const s128_o & lhs, const s128_o & rhs ) { if (lhs.high > rhs.high) return 1; if (rhs.high > lhs.high) return 0; // high components are equal if (lhs.low > rhs.low) return 1; return 0; } int operator > ( const s128_o & lhs, s32 rhs ) { s128_o t; t = rhs; if (lhs.high > t.high) return 1; if (t.high > lhs.high) return 0; // high components are equal if (lhs.low > t.low) return 1; return 0; } int operator >= ( const s128_o & lhs, const s128_o & rhs ) { if ((lhs.high == rhs.high) && (lhs.low == rhs.low)) return 1; if (lhs.high > rhs.high) return 1; if (rhs.high > lhs.high) return 0; // high components are equal if (lhs.low > rhs.low) return 1; return 0; } // Conversion for output u64 s128_u64(s128 x) { return(x.low); } u32 s128_u32(s128 x) { u32 rv; rv = x.low; return(rv); } void s128_str(s128 x, char *s) { int d, nd, going, fl1; s128 t, p10; fl1 = 0; if (x < ((s128_o) 0)) { *s = '-'; s++; x = - x; } if (x < ((s128_o) 0)) { /* It is still negative! That means we're dealing with the "MIN_INT" value, which is -(2^127). This prevents the rest of our algorithm from working, so we change it to 1-2^127, negate again to get (2^127)-1, and set a flag telling ourselves to add 1 to the final digit (from a 7 to an 8) */ x = x + 1; x = - x; fl1 = 1; } if (x == ((s128_o) 0)) { *s = '0'; s++; *s = 0; /* Terminating null */ return; } // Count number of digits in x p10 = 1; nd = 0; going = 1; while((p10 <= x) && (going)) { p10 = p10 + p10; if (p10 < ((s128_o) 0)) { going = 0; } t = p10; p10 = p10 + p10; if (p10 < ((s128_o) 0)) { going = 0; } p10 = p10 + p10; if (p10 < ((s128_o) 0)) { going = 0; } p10 = p10 + t; if (p10 < ((s128_o) 0)) { going = 0; } nd++; } // Extract each digit while (nd > 0) { int i; nd--; p10 = 1; for(i=0; i= p10) { x = x - p10; d++; } if (fl1 && (nd == 0)) { d++; } *s=('0' + d); s++; } *s=0; /* Terminating null */ } #endif /* end of int128.c */