developer tip

두 개의 큰 정수를 곱하는 동안 오버플로를 잡아서 계산

copycodes 2020. 12. 3. 08:03
반응형

두 개의 큰 정수를 곱하는 동안 오버플로를 잡아서 계산


비교적 큰 수를 곱하고 결과를 하나 또는 여러 정수에 저장하는 효율적인 (선택적으로 표준, 우아하고 구현하기 쉬운) 솔루션을 찾고 있습니다.

다음과 같이 선언 된 두 개의 64 비트 정수가 있다고 가정 해 보겠습니다.

uint64_t a = xxx, b = yyy; 

할 때 a * b작업으로 인해 오버플로가 발생하는지 어떻게 감지하고이 경우 캐리를 어딘가에 보관할 수 있습니까?

것을 제발 참고 내가 어떤 큰 숫자 라이브러리를 사용하지 않으려는 내가 숫자를 저장하는 방법에 제약이 있기 때문이다.


1. 오버플로 감지 :

x = a * b;
if (a != 0 && x / a != b) {
    // overflow handling
}

편집 : 고정 분할 0(마크 감사합니다!)

2. 캐리 계산 은 상당히 복잡합니다. 한 가지 방법은 두 피연산자를 반 단어로 분할 한 다음 반 단어에 긴 곱셈적용하는 것 입니다 .

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)

하자 B = 1 << 32. 그런 다음 우리는

            x <= (B - 1) + (B - 1)(B - 1) + (B - 1)
              <= B*B - 1
               < B*B

나는 이것이 작동 할 것이라고 믿습니다-적어도 그것은 Sjlver의 테스트 케이스를 처리합니다. 그 외에도 테스트되지 않았습니다 (더 이상 C ++ 컴파일러가 없기 때문에 컴파일하지 않을 수도 있습니다).


아이디어는 적분 연산에 대해 사실 인 다음 사실을 사용하는 것입니다.

a*b > c 경우에만 a > c/b

/ 여기서 적분 분할입니다.

양수의 오버플로를 검사하는 의사 코드는 다음과 같습니다.

if (a> max_int64 / b) then "overflow"else "ok" .

0과 음수를 처리하려면 더 많은 검사를 추가해야합니다.

음수가 아닌 C 코드는 다음 ab같습니다.

if (b > 0 && a > 18446744073709551615 / b) {
     // overflow handling
}; else {
    c = a * b;
}

노트 :

18446744073709551615 == (1<<64)-1

캐리를 계산하기 위해 우리는 숫자를 두 개의 32 자리로 나누고 종이에 이렇게 곱하는 방법을 사용할 수 있습니다. 오버플로를 피하기 위해 숫자를 분할해야합니다.

코드는 다음과 같습니다.

// 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

이 질문에 대한 다른 답변이 여러 개 있었지만 그중 몇 가지는 완전히 테스트되지 않은 코드가 있으며 지금까지 가능한 다른 옵션을 적절히 비교 한 사람은 없습니다.

이러한 이유로 저는 몇 가지 가능한 구현을 작성하고 테스트했습니다 (마지막 구현 은 OpenBSD 의이 코드기반으로하며 여기 에서 Reddit 에서 논의 됨 ). 코드는 다음과 같습니다.

/* 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);
}

다음은 내가 보유한 다양한 컴파일러 및 시스템으로 테스트 한 결과입니다 (이 경우 모든 테스트는 OS X에서 수행되었지만 결과는 BSD 또는 Linux 시스템에서 유사해야합니다).

+------------------+----------+----------+----------+----------+----------+
|                  | 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 |
+------------------+----------+----------+----------+----------+----------+

이러한 결과를 바탕으로 몇 가지 결론을 도출 할 수 있습니다.

  • 분명히, 분할 기반 접근 방식은 간단하고 이식 가능하지만 느립니다.
  • 어떤 기술도 모든 경우에 확실한 승자는 없습니다.
  • 최신 컴파일러에서는 사용할 수 있다면 더 큰 int 사용 방법이 가장 좋습니다.
  • 오래된 컴파일러에서는 긴 곱셈 방식이 가장 좋습니다.
  • 놀랍게도 GCC 4.9.0은 GCC 4.2.1에 비해 성능 회귀가 있고 GCC 4.2.1은 GCC 3.3에 비해 성능 회귀가 있습니다.

a == 0 일 때도 작동하는 버전 :

    x = a * b;
    if (a != 0 && x / a != b) {
        // overflow handling
    }

오버플로를 감지 할뿐만 아니라 캐리를 캡처해야하는 경우 숫자를 32 비트 부분으로 나누는 것이 가장 좋습니다. 코드는 악몽입니다. 다음은 단지 스케치입니다.

#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...
}

문제는 부분 곱만이 아니라 합계가 오버플로 될 수 있다는 사실입니다.

이 작업을 실제로 수행해야한다면 로컬 어셈블리 언어로 확장 곱하기 루틴을 작성합니다. 즉, 예를 들어 두 개의 64 비트 정수를 곱하여 두 개의 64 비트 레지스터에 저장된 128 비트 결과를 얻습니다. 합리적인 모든 하드웨어는 단일 기본 곱하기 명령어로이 기능을 제공합니다. C에서만 액세스 할 수있는 것은 아닙니다.

이것은 가장 우아하고 프로그래밍하기 쉬운 솔루션이 실제로 어셈블리 언어를 사용하는 드문 경우 중 하나입니다. 그러나 그것은 확실히 이식성이 없습니다 :-(


아마도이 문제를 해결하는 가장 좋은 방법은 두 개의 UInt64를 곱하여 UInt64 쌍, UInt128 결과의 위쪽 부분과 아래쪽 부분을 생성하는 함수를 갖는 것입니다. 다음은 결과를 16 진수로 표시하는 함수를 포함한 솔루션입니다. 아마도 C ++ 솔루션을 선호 할 것 같지만 문제를 관리하는 방법을 보여주는 Swift-Solution이 작동합니다.

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))

저는 요즘이 문제를 다루고 있는데, 사람들이 오버플로가 있는지 알 수있는 가장 좋은 방법은 결과를 나누는 것입니다. 그것은 완전히 비효율적이며 불필요한. 이 기능의 요점은 가능한 한 빨라야한다는 것입니다.

오버플로 감지에는 두 가지 옵션이 있습니다.

1º- 가능한 경우 곱셈기보다 두 배 큰 결과 변수를 만듭니다. 예를 들면 다음과 같습니다.

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)
}

오버플로가 있었는지 즉시 알 수 있으며 코드는 기계 코드로 작성하지 않고도 가능한 가장 빠릅니다. 컴파일러에 따라이 코드는 기계 코드에서 향상 될 수 있습니다.

2º- multipliers 변수보다 두 배 큰 결과 변수를 생성하는 것은 불가능합니다. 그러면 if 조건을 사용하여 최상의 경로를 결정해야합니다. 예를 들어 계속 :

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);
}
}

이 코드가 매우 효율적인 프로그램을 만드는 데 도움이되기를 바라며, 코드가 명확하기를 바랍니다.

친애하는.


clang 및 gcc로 쉽고 빠르게 :

unsigned long long t a, b, result;
if (__builtin_umulll_overflow(a, b, &result)) {
    // overflow!!
}

가능한 경우 오버플로 감지를 위해 하드웨어 지원을 사용합니다. 컴파일러 확장으로 C ++에서는 정의되지 않은 동작이지만 부호있는 정수 오버플로 (umul을 smul로 대체)를 처리 할 수도 있습니다.


다음은 두 개의 부호없는 정수의 곱셈이 오버플로되는지 여부를 감지하는 방법입니다.

N 비트의 이진수와 M 비트의 이진수를 곱하면 곱은 N + M 비트보다 많지 않다는 관찰을합니다.

예를 들어, 3 비트 숫자와 29 비트 숫자를 곱하라는 요청을 받으면 32 비트가 오버플로 되지 않는다는 것을 알고 있습니다.

#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 form\n");
    return EXIT_FAILURE;
  }

  a = strtoul(argv[1], &endptr, 0);

  if (*endptr != 0) {
    printf("%s is garbage\n", argv[1]);
    return EXIT_FAILURE;
  }

  b = strtoul(argv[2], &endptr, 0);

  if (*endptr != 0) {
    printf("%s is garbage\n", argv[2]);
    return EXIT_FAILURE;
  }

  if (might_be_mul_oflow(a, b))
    printf("might be multiplication overflow\n");

  {
    unsigned long c = a * b;
    printf("%lu * %lu = %lu\n", a, b, c);
    if (a != 0 && c / a != b)
      printf("confirmed multiplication overflow\n");
  }

  return 0;
}

일련의 테스트 : (64 비트 시스템에서) :

$ ./uflow 0x3 0x3FFFFFFFFFFFFFFF
3 * 4611686018427387903 = 13835058055282163709

$ ./uflow 0x7 0x3FFFFFFFFFFFFFFF
곱셈 오버플로 일 수 있습니다.
7 * 4611686018427387903 = 13835058055282163705
확인 된 곱셈 오버플로

$ ./uflow 0x4 0x3FFFFFFFFFFFFFFF
곱셈 오버플로 일 수 있습니다.
4 * 4611686018427387903 = 18446744073709551612

$ ./uflow 0x5 0x3FFFFFFFFFFFFFFF
곱셈 오버플로 일 수 있습니다.
5 * 4611686018427387903 = 4611686018427387899
확인 된 곱셈 오버플로

The steps in might_be_mul_oflow are almost certainly slower than just doing the division test, at least on mainstream processors used in desktop workstations, servers and mobile devices. On chips without good division support, it could be useful.


It occurs to me that there is another way to do this early rejection test.

  1. We start with a pair of numbers arng and brng which are initialized to 0x7FFF...FFFF and 1.

  2. If a <= arng and b <= brng we can conclude that there is no overflow.

  3. Otherwise, we shift arng to the right, and shift brng to the left, adding one bit to brng, so that they are 0x3FFF...FFFF and 3.

  4. If arng is zero, finish; otherwise repeat at 2.

The function now looks like:

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;
  }
}

If you just want to detect overflow, how about converting to double, doing the multiplication and if

|x| < 2^53, convert to int64

|x| < 2^63, make the multiplication using int64

otherwise produce whatever error you want?

This seems to work:

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;
}

참고URL : https://stackoverflow.com/questions/1815367/catch-and-compute-overflow-during-multiplication-of-two-large-integers

반응형