Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- uint64_t a = xxx, b = yyy;
- x = a * b;
- if (a != 0 && x / a != b) {
- // overflow handling
- }
- uint64_t hi(uint64_t x) {
- return x >> 32;
- }
- uint64_t lo(uint64_t x) {
- return ((1L << 32) - 1) & x;
- }
- void multiply(uint64_t a, uint64_t b) {
- // actually uint32_t would do, but the casting is annoying
- uint64_t s0, s1, s2, s3;
- uint64_t x = lo(a) * lo(b);
- s0 = lo(x);
- x = hi(a) * lo(b) + hi(x);
- s1 = lo(x);
- s2 = hi(x);
- x = s1 + lo(a) * hi(b);
- s1 = lo(x);
- x = s2 + hi(a) * hi(b) + hi(x);
- s2 = lo(x);
- s3 = hi(x);
- uint64_t result = s1 << 32 | s0;
- uint64_t carry = s3 << 32 | s2;
- }
- x = s2 + hi(a) * hi(b) + hi(x)
- x <= (B - 1) + (B - 1)(B - 1) + (B - 1)
- <= B*B - 1
- < B*B
- if (b > 0 && a > 18446744073709551615 / b) {
- // overflow handling
- }; else {
- c = a * b;
- }
- 18446744073709551615 == (1<<64)-1
- // split input numbers into 32-bit digits
- uint64_t a0 = a & ((1LL<<32)-1);
- uint64_t a1 = a >> 32;
- uint64_t b0 = b & ((1LL<<32)-1);
- uint64_t b1 = b >> 32;
- // The following 3 lines of code is to calculate the carry of d1
- // (d1 - 32-bit second digit of result, and it can be calculated as d1=d11+d12),
- // but to avoid overflow.
- // Actually rewriting the following 2 lines:
- // uint64_t d1 = (a0 * b0 >> 32) + a1 * b0 + a0 * b1;
- // uint64_t c1 = d1 >> 32;
- uint64_t d11 = a1 * b0 + (a0 * b0 >> 32);
- uint64_t d12 = a0 * b1;
- uint64_t c1 = (d11 > 18446744073709551615 - d12) ? 1 : 0;
- uint64_t d2 = a1 * b1 + c1;
- uint64_t carry = d2; // needed carry stored here
- /* Multiply with overflow checking, emulating clang's builtin function
- *
- * __builtin_umull_overflow
- *
- * This code benchmarks five possible schemes for doing so.
- */
- #include <stddef.h>
- #include <stdio.h>
- #include <stdlib.h>
- #include <stdint.h>
- #include <limits.h>
- #ifndef BOOL
- #define BOOL int
- #endif
- // Option 1, check for overflow a wider type
- // - Often fastest and the least code, especially on modern compilers
- // - When long is a 64-bit int, requires compiler support for 128-bits
- // ints (requires GCC >= 3.0 or Clang)
- #if LONG_BIT > 32
- typedef __uint128_t long_overflow_t ;
- #else
- typedef uint64_t long_overflow_t;
- #endif
- BOOL
- umull_overflow1(unsigned long lhs, unsigned long rhs, unsigned long* result)
- {
- long_overflow_t prod = (long_overflow_t)lhs * (long_overflow_t)rhs;
- *result = (unsigned long) prod;
- return (prod >> LONG_BIT) != 0;
- }
- // Option 2, perform long multiplication using a smaller type
- // - Sometimes the fastest (e.g., when mulitply on longs is a library
- // call).
- // - Performs at most three multiplies, and sometimes only performs one.
- // - Highly portable code; works no matter how many bits unsigned long is
- BOOL
- umull_overflow2(unsigned long lhs, unsigned long rhs, unsigned long* result)
- {
- const unsigned long HALFSIZE_MAX = (1ul << LONG_BIT/2) - 1ul;
- unsigned long lhs_high = lhs >> LONG_BIT/2;
- unsigned long lhs_low = lhs & HALFSIZE_MAX;
- unsigned long rhs_high = rhs >> LONG_BIT/2;
- unsigned long rhs_low = rhs & HALFSIZE_MAX;
- unsigned long bot_bits = lhs_low * rhs_low;
- if (!(lhs_high || rhs_high)) {
- *result = bot_bits;
- return 0;
- }
- BOOL overflowed = lhs_high && rhs_high;
- unsigned long mid_bits1 = lhs_low * rhs_high;
- unsigned long mid_bits2 = lhs_high * rhs_low;
- *result = bot_bits + ((mid_bits1+mid_bits2) << LONG_BIT/2);
- return overflowed || *result < bot_bits
- || (mid_bits1 >> LONG_BIT/2) != 0
- || (mid_bits2 >> LONG_BIT/2) != 0;
- }
- // Option 3, perform long multiplication using a smaller type (this code is
- // very similar to option 2, but calculates overflow using a different but
- // equivalent method).
- // - Sometimes the fastest (e.g., when mulitply on longs is a library
- // call; clang likes this code).
- // - Performs at most three multiplies, and sometimes only performs one.
- // - Highly portable code; works no matter how many bits unsigned long is
- BOOL
- umull_overflow3(unsigned long lhs, unsigned long rhs, unsigned long* result)
- {
- const unsigned long HALFSIZE_MAX = (1ul << LONG_BIT/2) - 1ul;
- unsigned long lhs_high = lhs >> LONG_BIT/2;
- unsigned long lhs_low = lhs & HALFSIZE_MAX;
- unsigned long rhs_high = rhs >> LONG_BIT/2;
- unsigned long rhs_low = rhs & HALFSIZE_MAX;
- unsigned long lowbits = lhs_low * rhs_low;
- if (!(lhs_high || rhs_high)) {
- *result = lowbits;
- return 0;
- }
- BOOL overflowed = lhs_high && rhs_high;
- unsigned long midbits1 = lhs_low * rhs_high;
- unsigned long midbits2 = lhs_high * rhs_low;
- unsigned long midbits = midbits1 + midbits2;
- overflowed = overflowed || midbits < midbits1 || midbits > HALFSIZE_MAX;
- unsigned long product = lowbits + (midbits << LONG_BIT/2);
- overflowed = overflowed || product < lowbits;
- *result = product;
- return overflowed;
- }
- // Option 4, checks for overflow using division
- // - Checks for overflow using division
- // - Division is slow, especially if it is a library call
- BOOL
- umull_overflow4(unsigned long lhs, unsigned long rhs, unsigned long* result)
- {
- *result = lhs * rhs;
- return rhs > 0 && (SIZE_MAX / rhs) < lhs;
- }
- // Option 5, checks for overflow using division
- // - Checks for overflow using division
- // - Avoids division when the numbers are "small enough" to trivially
- // rule out overflow
- // - Division is slow, especially if it is a library call
- BOOL
- umull_overflow5(unsigned long lhs, unsigned long rhs, unsigned long* result)
- {
- const unsigned long MUL_NO_OVERFLOW = (1ul << LONG_BIT/2) - 1ul;
- *result = lhs * rhs;
- return (lhs >= MUL_NO_OVERFLOW || rhs >= MUL_NO_OVERFLOW) &&
- rhs > 0 && SIZE_MAX / rhs < lhs;
- }
- #ifndef umull_overflow
- #define umull_overflow2
- #endif
- /*
- * This benchmark code performs a multiply at all bit sizes,
- * essentially assuming that sizes are logarithmically distributed.
- */
- int main()
- {
- unsigned long i, j, k;
- int count = 0;
- unsigned long mult;
- unsigned long total = 0;
- for (k = 0; k < 0x40000000 / LONG_BIT / LONG_BIT; ++k)
- for (i = 0; i != LONG_MAX; i = i*2+1)
- for (j = 0; j != LONG_MAX; j = j*2+1) {
- count += umull_overflow(i+k, j+k, &mult);
- total += mult;
- }
- printf("%d overflows (total %lu)n", count, total);
- }
- +------------------+----------+----------+----------+----------+----------+
- | | Option 1 | Option 2 | Option 3 | Option 4 | Option 5 |
- | | BigInt | LngMult1 | LngMult2 | Div | OptDiv |
- +------------------+----------+----------+----------+----------+----------+
- | Clang 3.5 i386 | 1.610 | 3.217 | 3.129 | 4.405 | 4.398 |
- | GCC 4.9.0 i386 | 1.488 | 3.469 | 5.853 | 4.704 | 4.712 |
- | GCC 4.2.1 i386 | 2.842 | 4.022 | 3.629 | 4.160 | 4.696 |
- | GCC 4.2.1 PPC32 | 8.227 | 7.756 | 7.242 | 20.632 | 20.481 |
- | GCC 3.3 PPC32 | 5.684 | 9.804 | 11.525 | 21.734 | 22.517 |
- +------------------+----------+----------+----------+----------+----------+
- | Clang 3.5 x86_64 | 1.584 | 2.472 | 2.449 | 9.246 | 7.280 |
- | GCC 4.9 x86_64 | 1.414 | 2.623 | 4.327 | 9.047 | 7.538 |
- | GCC 4.2.1 x86_64 | 2.143 | 2.618 | 2.750 | 9.510 | 7.389 |
- | GCC 4.2.1 PPC64 | 13.178 | 8.994 | 8.567 | 37.504 | 29.851 |
- +------------------+----------+----------+----------+----------+----------+
- x = a * b;
- if (a != 0 && x / a != b) {
- // overflow handling
- }
- #include <stdint.h>
- uint64_t mul(uint64_t a, uint64_t b) {
- uint32_t ah = a >> 32;
- uint32_t al = a; // truncates: now a = al + 2**32 * ah
- uint32_t bh = b >> 32;
- uint32_t bl = b; // truncates: now b = bl + 2**32 * bh
- // a * b = 2**64 * ah * bh + 2**32 * (ah * bl + bh * al) + al * bl
- uint64_t partial = (uint64_t) al * (uint64_t) bl;
- uint64_t mid1 = (uint64_t) ah * (uint64_t) bl;
- uint64_t mid2 = (uint64_t) al * (uint64_t) bh;
- uint64_t carry = (uint64_t) ah * (uint64_t) bh;
- // add high parts of mid1 and mid2 to carry
- // add low parts of mid1 and mid2 to partial, carrying
- // any carry bits into carry...
- }
- struct INT32struct {INT16 high, low;};
- typedef union
- {
- struct INT32struct s;
- INT32 ll;
- } INT32union;
- INT16 mulFunction(INT16 a, INT16 b)
- {
- INT32union result.ll = a * b; //32Bits result
- if(result.s.high > 0)
- Overflow();
- return (result.s.low)
- }
- INT32 mulFunction(INT32 a, INT32 b)
- {
- INT32union s_a.ll = abs(a);
- INT32union s_b.ll = abs(b); //32Bits result
- INT32union result;
- if(s_a.s.hi > 0 && s_b.s.hi > 0)
- {
- Overflow();
- }
- else if (s_a.s.hi > 0)
- {
- INT32union res1.ll = s_a.s.hi * s_b.s.lo;
- INT32union res2.ll = s_a.s.lo * s_b.s.lo;
- if (res1.hi == 0)
- {
- result.s.lo = res1.s.lo + res2.s.hi;
- if (result.s.hi == 0)
- {
- result.s.ll = result.s.lo << 16 + res2.s.lo;
- if ((a.s.hi >> 15) ^ (b.s.hi >> 15) == 1)
- {
- result.s.ll = -result.s.ll;
- }
- return result.s.ll
- }else
- {
- Overflow();
- }
- }else
- {
- Overflow();
- }
- }else if (s_b.s.hi > 0)
- {
- //Same code changing a with b
- }else
- {
- return (s_a.lo * s_b.lo);
- }
- }
- func hex128 (_ hi: UInt64, _ lo: UInt64) -> String
- {
- var s: String = String(format: "%08X", hi >> 32)
- + String(format: "%08X", hi & 0xFFFFFFFF)
- + String(format: "%08X", lo >> 32)
- + String(format: "%08X", lo & 0xFFFFFFFF)
- return (s)
- }
- func mul64to128 (_ multiplier: UInt64, _ multiplicand : UInt64)
- -> (result_hi: UInt64, result_lo: UInt64)
- {
- let x: UInt64 = multiplier
- let x_lo: UInt64 = (x & 0xffffffff)
- let x_hi: UInt64 = x >> 32
- let y: UInt64 = multiplicand
- let y_lo: UInt64 = (y & 0xffffffff)
- let y_hi: UInt64 = y >> 32
- let mul_lo: UInt64 = (x_lo * y_lo)
- let mul_hi: UInt64 = (x_hi * y_lo) + (mul_lo >> 32)
- let mul_carry: UInt64 = (x_lo * y_hi) + (mul_hi & 0xffffffff)
- let result_hi: UInt64 = (x_hi * y_hi) + (mul_hi >> 32) + (mul_carry >> 32)
- let result_lo: UInt64 = (mul_carry << 32) + (mul_lo & 0xffffffff)
- return (result_hi, result_lo)
- }
- var c: UInt64 = 0
- var d: UInt64 = 0
- (c, d) = mul64to128(0x1234567890123456, 0x9876543210987654)
- // 0AD77D742CE3C72E45FD10D81D28D038 is the result of the above example
- print(hex128(c, d))
- (c, d) = mul64to128(0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF)
- // FFFFFFFFFFFFFFFE0000000000000001 is the result of the above example
- print(hex128(c, d))
- #include <stdlib.h>
- #include <stdio.h>
- int might_be_mul_oflow(unsigned long a, unsigned long b)
- {
- if (!a || !b)
- return 0;
- a = a | (a >> 1) | (a >> 2) | (a >> 4) | (a >> 8) | (a >> 16) | (a >> 32);
- b = b | (b >> 1) | (b >> 2) | (b >> 4) | (b >> 8) | (b >> 16) | (b >> 32);
- for (;;) {
- unsigned long na = a << 1;
- if (na <= a)
- break;
- a = na;
- }
- return (a & b) ? 1 : 0;
- }
- int main(int argc, char **argv)
- {
- unsigned long a, b;
- char *endptr;
- if (argc < 3) {
- printf("supply two unsigned long integers in C formn");
- return EXIT_FAILURE;
- }
- a = strtoul(argv[1], &endptr, 0);
- if (*endptr != 0) {
- printf("%s is garbagen", argv[1]);
- return EXIT_FAILURE;
- }
- b = strtoul(argv[2], &endptr, 0);
- if (*endptr != 0) {
- printf("%s is garbagen", argv[2]);
- return EXIT_FAILURE;
- }
- if (might_be_mul_oflow(a, b))
- printf("might be multiplication overflown");
- {
- unsigned long c = a * b;
- printf("%lu * %lu = %lun", a, b, c);
- if (a != 0 && c / a != b)
- printf("confirmed multiplication overflown");
- }
- return 0;
- }
- might_be_mul_oflow
- int might_be_mul_oflow(unsigned long a, unsigned long b)
- {
- if (!a || !b)
- return 0;
- {
- unsigned long arng = ULONG_MAX >> 1;
- unsigned long brng = 1;
- while (arng != 0) {
- if (a <= arng && b <= brng)
- return 0;
- arng >>= 1;
- brng <<= 1;
- brng |= 1;
- }
- return 1;
- }
- }
- int64_t safemult(int64_t a, int64_t b) {
- double dx;
- dx = (double)a * (double)b;
- if ( fabs(dx) < (double)9007199254740992 )
- return (int64_t)dx;
- if ( (double)INT64_MAX < fabs(dx) )
- return INT64_MAX;
- return a*b;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement